Add basic support for WCS pixelization.

.github/workflows/test.yml

@@ -71,8 +71,8 @@ jobs:
       fail-fast: false
       matrix:
         include:
-          - python: "3.7"
-            pyshort: "37"
+          # - python: "3.8"
+          #   pyshort: "38"
           - python: "3.9"
             pyshort: "39"
           # - python: "3.10"

setup.py

@@ -275,6 +275,7 @@ def readme():
     "pyyaml",
     "astropy",
     "healpy",
+    "pixell",
     "ephem",
 ]
 conf["extras_require"] = {

src/toast/CMakeLists.txt

@@ -110,7 +110,8 @@ install(FILES
     config.py
     job.py
     pixels.py
-    pixels_io.py
+    pixels_io_healpix.py
+    pixels_io_wcs.py
     covariance.py
     dist.py
     data.py

src/toast/_libtoast/ops_mapmaker_utils.cpp

@@ -265,7 +265,7 @@ void init_ops_mapmaker_utils(py::module & m) {
                                     use_det_flags
                                 );
                                 for (int64_t iw = 0; iw < nnz; iw++) {
-                                    #pragma omp atomic
+                                    # pragma omp atomic
                                     raw_zmap[zoff + iw] += zmap_val[iw];
                                 }
                             }

src/toast/_libtoast/ops_pixels_healpix.cpp

@@ -27,7 +27,7 @@
 // So we must duplicate this across compilation units.
 
 void pixels_healpix_qa_rotate(double const * q_in, double const * v_in,
-               double * v_out) {
+                              double * v_out) {
     // The input quaternion has already been normalized on the host.
 
     double xw =  q_in[3] * q_in[0];

src/toast/_libtoast/ops_stokes_weights.cpp

@@ -19,7 +19,7 @@
 // So we must duplicate this across compilation units.
 
 void stokes_weights_qa_rotate(double const * q_in, double const * v_in,
-               double * v_out) {
+                              double * v_out) {
     // The input quaternion has already been normalized on the host.
 
     double xw =  q_in[3] * q_in[0];

src/toast/_libtoast/template_offset.cpp

@@ -201,7 +201,7 @@ void init_template_offset(py::module & m) {
                             } else {
                                 contrib = raw_det_data[d];
                             }
-                            #pragma omp atomic
+                            # pragma omp atomic
                             raw_amplitudes[amp] += contrib;
                         }
                         offset += raw_n_amp_views[iview];

src/toast/dipole.py

@@ -2,7 +2,6 @@
 # All rights reserved.  Use of this source code is governed by
 # a BSD-style license that can be found in the LICENSE file.
 
-import healpy as hp
 import numpy as np
 import scipy.constants as constants
 from astropy import units as u

src/toast/observation.py

@@ -185,12 +185,13 @@ def __init__(
             self._uid = name_UID(self._name)
 
         if self._session is None:
-            self._session = Session(
-                name=self._name,
-                uid=self._uid,
-                start=None,
-                end=None,
-            )
+            if self._name is not None:
+                self._session = Session(
+                    name=self._name,
+                    uid=self._uid,
+                    start=None,
+                    end=None,
+                )
         elif not isinstance(self._session, Session):
             raise RuntimeError("session should be a Session instance or None")
 

src/toast/ops/CMakeLists.txt

@@ -27,6 +27,7 @@ install(FILES
     pointing_detector.py
     scan_map.py
     scan_healpix.py
+    scan_wcs.py
     mapmaker_utils.py
     mapmaker_binning.py
     mapmaker_solve.py
@@ -53,6 +54,7 @@ install(FILES
     demodulation.py
     stokes_weights.py
     pixels_healpix.py
+    pixels_wcs.py
     filterbin.py
     save_hdf5.py
     load_hdf5.py

src/toast/ops/__init__.py

@@ -36,6 +36,7 @@
 from .operator import Operator
 from .pipeline import Pipeline
 from .pixels_healpix import PixelsHealpix
+from .pixels_wcs import PixelsWCS
 from .pointing import BuildPixelDistribution
 from .pointing_detector import PointingDetectorSimple
 from .polyfilter import CommonModeFilter, PolyFilter, PolyFilter2D
@@ -45,6 +46,7 @@
 from .save_spt3g import SaveSpt3g
 from .scan_healpix import ScanHealpixMap, ScanHealpixMask
 from .scan_map import ScanMap, ScanMask, ScanScale
+from .scan_wcs import ScanWCSMap, ScanWCSMask
 from .sim_cosmic_rays import InjectCosmicRays
 from .sim_crosstalk import CrossTalk, MitigateCrossTalk
 from .sim_gaindrifts import GainDrifter

src/toast/ops/crosslinking.py

@@ -12,7 +12,7 @@
 from ..mpi import MPI
 from ..observation import default_values as defaults
 from ..pixels import PixelData, PixelDistribution
-from ..pixels_io import write_healpix_fits
+from ..pixels_io_healpix import write_healpix_fits
 from ..timing import Timer, function_timer
 from ..traits import Bool, Float, Instance, Int, Unicode, trait_docs
 from ..utils import Logger

src/toast/ops/filterbin.py

@@ -22,7 +22,7 @@
 from ..mpi import MPI
 from ..observation import default_values as defaults
 from ..pixels import PixelData, PixelDistribution
-from ..pixels_io import (
+from ..pixels_io_healpix import (
     filename_is_fits,
     filename_is_hdf5,
     read_healpix_fits,

src/toast/ops/groundfilter.py

@@ -4,7 +4,6 @@
 
 from time import time
 
-import healpy as hp
 import numpy as np
 import traitlets
 from astropy import units as u

src/toast/ops/madam_utils.py

@@ -4,7 +4,6 @@
 
 import os
 
-import healpy as hp
 import numpy as np
 
 from ..timing import function_timer

src/toast/ops/mapmaker.py

@@ -10,7 +10,8 @@
 from ..mpi import MPI
 from ..observation import default_values as defaults
 from ..pixels import PixelData, PixelDistribution
-from ..pixels_io import write_healpix_fits, write_healpix_hdf5
+from ..pixels_io_healpix import write_healpix_fits, write_healpix_hdf5
+from ..pixels_io_wcs import write_wcs_fits
 from ..timing import Timer, function_timer
 from ..traits import Bool, Float, Instance, Int, Unicode, trait_docs
 from ..utils import Logger
@@ -145,7 +146,7 @@ class MapMaker(Operator):
     write_rcond = Bool(True, help="If True, write the reciprocal condition numbers.")
 
     write_solver_products = Bool(
-        True, help="If True, write out equivalent solver products."
+        False, help="If True, write out equivalent solver products."
     )
 
     keep_solver_products = Bool(
@@ -423,9 +424,14 @@ def _exec(self, data, detectors=None, **kwargs):
                 ]
             ).apply(data)
 
-        # FIXME:  This all assumes the pointing operator is an instance of the
-        # PointingHealpix class.  We need to generalize distributed pixel data
-        # formats and associate them with the pointing operator.
+        # FIXME:  This I/O technique assumes "known" types of pixel representations.
+        # Instead, we should associate read / write functions to a particular pixel
+        # class.
+
+        is_pix_wcs = hasattr(map_binning.pixel_pointing, "wcs")
+        is_hpix_nest = None
+        if not is_pix_wcs:
+            is_hpix_nest = map_binning.pixel_pointing.nest
 
         write_del = list()
         write_del.append((self.hits_name, self.write_hits))
@@ -439,25 +445,31 @@ def _exec(self, data, detectors=None, **kwargs):
         wtimer.start()
         for prod_key, prod_write in write_del:
             if prod_write:
-                if self.write_hdf5:
-                    # Non-standard HDF5 output
-                    fname = os.path.join(self.output_dir, "{}.h5".format(prod_key))
-                    write_healpix_hdf5(
-                        data[prod_key],
-                        fname,
-                        nest=map_binning.pixel_pointing.nest,
-                        single_precision=True,
-                        force_serial=self.write_hdf5_serial,
-                    )
-                else:
-                    # Standard FITS output
+                if is_pix_wcs:
                     fname = os.path.join(self.output_dir, "{}.fits".format(prod_key))
-                    write_healpix_fits(
-                        data[prod_key],
-                        fname,
-                        nest=map_binning.pixel_pointing.nest,
-                        report_memory=self.report_memory,
-                    )
+                    write_wcs_fits(data[prod_key], fname)
+                else:
+                    if self.write_hdf5:
+                        # Non-standard HDF5 output
+                        fname = os.path.join(self.output_dir, "{}.h5".format(prod_key))
+                        write_healpix_hdf5(
+                            data[prod_key],
+                            fname,
+                            nest=is_hpix_nest,
+                            single_precision=True,
+                            force_serial=self.write_hdf5_serial,
+                        )
+                    else:
+                        # Standard FITS output
+                        fname = os.path.join(
+                            self.output_dir, "{}.fits".format(prod_key)
+                        )
+                        write_healpix_fits(
+                            data[prod_key],
+                            fname,
+                            nest=is_hpix_nest,
+                            report_memory=self.report_memory,
+                        )
                 log.info_rank(f"Wrote {fname} in", comm=comm, timer=wtimer)
             if not self.keep_final_products:
                 if prod_key in data:

src/toast/ops/mapmaker_templates.py

@@ -11,7 +11,8 @@
 from ..mpi import MPI
 from ..observation import default_values as defaults
 from ..pixels import PixelData, PixelDistribution
-from ..pixels_io import write_healpix_fits, write_healpix_hdf5
+from ..pixels_io_healpix import write_healpix_fits, write_healpix_hdf5
+from ..pixels_io_wcs import write_wcs_fits
 from ..templates import AmplitudesMap, Template
 from ..timing import Timer, function_timer
 from ..traits import Bool, Float, Instance, Int, List, Unicode, trait_docs
@@ -785,6 +786,15 @@ def _exec(self, data, detectors=None, **kwargs):
         memreport.prefix = "After solving for amplitudes"
         memreport.apply(data)
 
+        # FIXME:  This I/O technique assumes "known" types of pixel representations.
+        # Instead, we should associate read / write functions to a particular pixel
+        # class.
+
+        is_pix_wcs = hasattr(self.binning.pixel_pointing, "wcs")
+        is_hpix_nest = None
+        if not is_pix_wcs:
+            is_hpix_nest = self.binning.pixel_pointing.nest
+
         write_del = [
             self.solver_hits_name,
             self.solver_cov_name,
@@ -794,25 +804,31 @@ def _exec(self, data, detectors=None, **kwargs):
         ]
         for prod_key in write_del:
             if self.write_solver_products:
-                if self.write_hdf5:
-                    # Non-standard HDF5 output
-                    fname = os.path.join(self.output_dir, "{}.h5".format(prod_key))
-                    write_healpix_hdf5(
-                        data[prod_key],
-                        fname,
-                        nest=self.binning.pixel_pointing.nest,
-                        single_precision=True,
-                        force_serial=self.write_hdf5_serial,
-                    )
-                else:
-                    # Standard FITS output
+                if is_pix_wcs:
                     fname = os.path.join(self.output_dir, "{}.fits".format(prod_key))
-                    write_healpix_fits(
-                        data[prod_key],
-                        fname,
-                        nest=self.binning.pixel_pointing.nest,
-                        report_memory=self.report_memory,
-                    )
+                    write_wcs_fits(data[prod_key], fname)
+                else:
+                    if self.write_hdf5:
+                        # Non-standard HDF5 output
+                        fname = os.path.join(self.output_dir, "{}.h5".format(prod_key))
+                        write_healpix_hdf5(
+                            data[prod_key],
+                            fname,
+                            nest=is_hpix_nest,
+                            single_precision=True,
+                            force_serial=self.write_hdf5_serial,
+                        )
+                    else:
+                        # Standard FITS output
+                        fname = os.path.join(
+                            self.output_dir, "{}.fits".format(prod_key)
+                        )
+                        write_healpix_fits(
+                            data[prod_key],
+                            fname,
+                            nest=is_hpix_nest,
+                            report_memory=self.report_memory,
+                        )
             if not self.mc_mode and not self.keep_solver_products:
                 if prod_key in data:
                     data[prod_key].clear()

src/toast/ops/noise_estimation.py

@@ -128,7 +128,8 @@ class NoiseEstim(Operator):
     )
 
     output_dir = Unicode(
-        ".",
+        None,
+        allow_none=True,
         help="If specified, write output data products to this directory",
     )
 
@@ -1075,21 +1076,21 @@ def process_noise_estimate(
             )
             log.debug_rank("Discard outliers", timer=timer)
 
-            self.save_psds(
-                binfreq, all_psds, all_times, det1, det2, fsample, fileroot, all_cov
-            )
-
-            if nbad > 0:
+            if self.output_dir is not None:
                 self.save_psds(
-                    binfreq,
-                    good_psds,
-                    good_times,
-                    det1,
-                    det2,
-                    fsample,
-                    fileroot + "_good",
-                    good_cov,
+                    binfreq, all_psds, all_times, det1, det2, fsample, fileroot, all_cov
                 )
+                if nbad > 0:
+                    self.save_psds(
+                        binfreq,
+                        good_psds,
+                        good_times,
+                        det1,
+                        det2,
+                        fsample,
+                        fileroot + "_good",
+                        good_cov,
+                    )
 
             final_freqs = binfreq
             final_psd = np.mean(np.array(good_psds), axis=0)