Add methods to write explicit skymap without building the BTreeMap

Author: fxpineau

src/nested/iter/sortedhash2hashcount.rs

@@ -255,10 +255,10 @@ where
 
 #[cfg(test)]
 mod tests {
-  use super::*;
-  use crate::depth;
   use std::iter;
 
+  use super::*;
+
   #[test]
   fn test_count_iterators_ok() {
     // Test coutn from empty iterators

src/nested/map/fits/write.rs

@@ -94,6 +94,37 @@ where
   }
   write_final_padding(writer, n_elems * naxis1)
 }
+// Made for python wrapper
+pub fn write_explicit_skymap_fits_from_parts<R, K, V>(
+  mut writer: R,
+  depth: u8,
+  n_elems: usize,
+  it_keys: K,
+  it_vals: V,
+) -> Result<(), FitsError>
+where
+  R: Write,
+  K: Iterator<Item = u64>,
+  V: Iterator,
+  V::Item: SkyMapValue,
+{
+  let naxis1 = (<u64 as HHash>::FITS_NAXIS1 + V::Item::FITS_NAXIS1) as usize;
+  write_explicit_skymap_fits_header_from_parts::<_, u64, V::Item>(&mut writer, depth, n_elems)?;
+  let mut n = 0;
+  for (h, v) in it_keys.zip(it_vals) {
+    writer
+      .write_all(h.to_be_bytes().as_ref())
+      .and_then(|()| writer.write_all(v.to_be_bytes().as_ref()))?;
+    n += 1;
+  }
+  if n != n_elems {
+    return Err(FitsError::new_custom(format!(
+      "Wrong number of HEALPix cells writen at depth {}. Expected: {}. Actual: {}.",
+      depth, n_elems, n
+    )));
+  }
+  write_final_padding(writer, n_elems * naxis1)
+}
 
 /// Possible add blanks at the end of the FITS file to complete the last
 /// 2880 bytes block.
@@ -216,15 +247,12 @@ fn write_implicit_skymap_fits_header<R: Write, T: SkyMapValue>(
   writer.write_all(&header_block[..]).map_err(FitsError::Io)
 }
 
-fn write_explicit_skymap_fits_header<'a, R, S>(
-  mut writer: R,
-  skymap: &'a S,
-) -> Result<(), FitsError>
+fn write_explicit_skymap_fits_header<'a, R, S>(writer: R, skymap: &'a S) -> Result<(), FitsError>
 where
   R: Write,
   S: SkyMap<'a>,
 {
-  let nside = nside(skymap.depth());
+  /*let nside = nside(skymap.depth());
   let n_cells = skymap.len() as u64;
 
   write_primary_hdu(&mut writer)?;
@@ -267,5 +295,65 @@ where
   );
   it.next().unwrap()[0..3].copy_from_slice(b"END");
   // Do write the header
+  writer.write_all(&header_block[..]).map_err(FitsError::Io)*/
+  write_explicit_skymap_fits_header_from_parts::<_, S::HashType, S::ValueType>(
+    writer,
+    skymap.depth(),
+    skymap.len(),
+  )
+}
+fn write_explicit_skymap_fits_header_from_parts<R, K, V>(
+  mut writer: R,
+  depth: u8,
+  n_elems: usize,
+) -> Result<(), FitsError>
+where
+  R: Write,
+  K: HHash,
+  V: SkyMapValue,
+{
+  let nside = nside(depth);
+  let n_cells = n_elems as u64;
+
+  write_primary_hdu(&mut writer)?;
+  let mut header_block = [b' '; 2880];
+  let mut it = header_block.chunks_mut(80);
+  // Write BINTABLE specific keywords in the buffer
+  it.next().unwrap()[0..20].copy_from_slice(b"XTENSION= 'BINTABLE'");
+  it.next().unwrap()[0..30].copy_from_slice(b"BITPIX  =                    8");
+  it.next().unwrap()[0..30].copy_from_slice(b"NAXIS   =                    2");
+  write_uint_mandatory_keyword_record(
+    it.next().unwrap(),
+    b"NAXIS1  ",
+    (K::FITS_NAXIS1 + V::FITS_NAXIS1) as u64,
+  );
+  write_uint_mandatory_keyword_record(it.next().unwrap(), b"NAXIS2  ", n_cells);
+  it.next().unwrap()[0..30].copy_from_slice(b"PCOUNT  =                    0");
+  it.next().unwrap()[0..30].copy_from_slice(b"GCOUNT  =                    1");
+  it.next().unwrap()[0..30].copy_from_slice(b"TFIELDS =                    2");
+  it.next().unwrap()[0..20].copy_from_slice(b"TTYPE1  = 'PIXEL   '");
+  it.next().unwrap()[0..20].copy_from_slice(b"TTYPE2  = 'VALUE   '");
+  write_str_mandatory_keyword_record(it.next().unwrap(), b"TFORM1  ", K::FITS_TFORM);
+  write_str_mandatory_keyword_record(it.next().unwrap(), b"TFORM2  ", V::FITS_TFORM);
+  it.next().unwrap()[0..20].copy_from_slice(b"PIXTYPE = 'HEALPIX '");
+  it.next().unwrap()[0..20].copy_from_slice(b"INDXSCHM= 'EXPLICIT'");
+  it.next().unwrap()[0..20].copy_from_slice(b"ORDERING= 'NESTED  '");
+  it.next().unwrap()[0..20].copy_from_slice(b"COORDSYS= 'C       '");
+  it.next().unwrap()[0..20].copy_from_slice(b"EXTNAME = 'xtension'");
+  write_uint_mandatory_keyword_record(it.next().unwrap(), b"NSIDE   ", nside as u64);
+  write_uint_mandatory_keyword_record(it.next().unwrap(), b"OBS_NPIX", n_cells);
+  // it.next().unwrap()[0..28].copy_from_slice(b"CREATOR = 'CDS HEALPix Rust'");
+  write_keyword_record(
+    it.next().unwrap(),
+    b"CREATOR ",
+    format!(
+      "'Rust crate {} {}'",
+      env!("CARGO_PKG_NAME"),
+      env!("CARGO_PKG_VERSION")
+    )
+    .as_str(),
+  );
+  it.next().unwrap()[0..3].copy_from_slice(b"END");
+  // Do write the header
   writer.write_all(&header_block[..]).map_err(FitsError::Io)
 }

src/nested/map/mod.rs

@@ -148,6 +148,6 @@ mod tests {
   fn test_read_skymap() {
     init_logger();
     let path = "test/resources/skymap/gaiadr3.nside64.densmap.fits";
-    let skymap = SkyMapEnum::from_fits_file(path).unwrap();
+    SkyMapEnum::from_fits_file(path).unwrap();
   }
 }

src/nested/mod.rs

@@ -5188,8 +5188,7 @@ mod tests {
   #[test]
   fn testok_cone_approx_bmoc_2() {
     // From https://github.com/cds-astro/cds-healpix-rust/issues/21
-    let actual_res = cone_coverage_approx(1, 0.0, -0.3956593324046537, 0.00986851403158301);
-
+    //let actual_res = cone_coverage_approx(1, 0.0, -0.3956593324046537, 0.00986851403158301);
     let actual_res = cone_coverage_approx(10, 0.0, 0.0, (5f64 / 3600.0).to_radians());
     println!("@@@@@ FLAT VIEW");
     for cell in actual_res.flat_iter() {