Merge pull request #609 from lsst/tickets/DM-31990

DM-31990: Add configuration option to set coverage nside for survey-wide HealSparsePropertyMaps

python/lsst/pipe/tasks/healSparseMapping.py

@@ -20,6 +20,7 @@
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
 from collections import defaultdict
+import warnings
 import numbers
 import numpy as np
 import healpy as hp
@@ -152,10 +153,16 @@ def build_ccd_input_map(self, bbox, wcs, ccds):
         ccds : `lsst.afw.table.ExposureCatalog`
             Exposure catalog with ccd data from coadd inputs.
         """
-        self.ccd_input_map = hsp.HealSparseMap.make_empty(nside_coverage=self.config.nside_coverage,
-                                                          nside_sparse=self.config.nside,
-                                                          dtype=hsp.WIDE_MASK,
-                                                          wide_mask_maxbits=len(ccds))
+        with warnings.catch_warnings():
+            # Healsparse will emit a warning if nside coverage is greater than
+            # 128.  In the case of generating patch input maps, and not global
+            # maps, high nside coverage works fine, so we can suppress this
+            # warning.
+            warnings.simplefilter("ignore")
+            self.ccd_input_map = hsp.HealSparseMap.make_empty(nside_coverage=self.config.nside_coverage,
+                                                              nside_sparse=self.config.nside,
+                                                              dtype=hsp.WIDE_MASK,
+                                                              wide_mask_maxbits=len(ccds))
         self._wcs = wcs
         self._bbox = bbox
         self._ccds = ccds
@@ -203,12 +210,18 @@ def build_ccd_input_map(self, bbox, wcs, ccds):
 
         dtype = [(f"v{visit}", "i4") for visit in self._bits_per_visit.keys()]
 
-        cov = self.config.nside_coverage
-        ns = self.config.nside
-        self._ccd_input_bad_count_map = hsp.HealSparseMap.make_empty(nside_coverage=cov,
-                                                                     nside_sparse=ns,
-                                                                     dtype=dtype,
-                                                                     primary=dtype[0][0])
+        with warnings.catch_warnings():
+            # Healsparse will emit a warning if nside coverage is greater than
+            # 128.  In the case of generating patch input maps, and not global
+            # maps, high nside coverage works fine, so we can suppress this
+            # warning.
+            warnings.simplefilter("ignore")
+            self._ccd_input_bad_count_map = hsp.HealSparseMap.make_empty(
+                nside_coverage=self.config.nside_coverage,
+                nside_sparse=self.config.nside,
+                dtype=dtype,
+                primary=dtype[0][0])
+
         # Don't set input bad map if there are no ccds which overlap the bbox.
         if len(self._ccd_input_pixels) > 0:
             self._ccd_input_bad_count_map[self._ccd_input_pixels] = np.zeros(1, dtype=dtype)
@@ -694,8 +707,9 @@ def _compute_nside_coverage_tract(self, tract_info):
         nside_coverage_tract = 32
         while hp.nside2pixarea(nside_coverage_tract, degrees=True) > tract_area:
             nside_coverage_tract = 2*nside_coverage_tract
-        # Step back one, but don't go bigger pixels than nside=32
-        nside_coverage_tract = int(np.clip(nside_coverage_tract/2, 32, None))
+        # Step back one, but don't go bigger pixels than nside=32 or smaller
+        # than 128 (recommended by healsparse).
+        nside_coverage_tract = int(np.clip(nside_coverage_tract/2, 32, 128))
 
         return nside_coverage_tract
 
@@ -838,6 +852,12 @@ class ConsolidateHealSparsePropertyMapConfig(pipeBase.PipelineTaskConfig,
                  "dcr_e2"],
         doc="Property map computation objects",
     )
+    nside_coverage = pexConfig.Field(
+        doc="Consolidated HealSparse coverage map nside.  Must be power of 2.",
+        dtype=int,
+        default=32,
+        check=_is_power_of_two,
+    )
 
     def setDefaults(self):
         self.property_maps["exposure_time"].do_sum = True
@@ -904,29 +924,44 @@ def consolidate_map(self, sky_map, input_refs):
         # First, we read in the coverage maps to know how much memory
         # to allocate
         cov_mask = None
+        nside_coverage_inputs = None
         for tract_id in input_refs:
             cov = input_refs[tract_id].get(component='coverage')
             if cov_mask is None:
                 cov_mask = cov.coverage_mask
+                nside_coverage_inputs = cov.nside_coverage
             else:
                 cov_mask |= cov.coverage_mask
 
-        cov_pix, = np.where(cov_mask)
+        cov_pix_inputs, = np.where(cov_mask)
+
+        # Compute the coverage pixels for the desired nside_coverage
+        if nside_coverage_inputs == self.config.nside_coverage:
+            cov_pix = cov_pix_inputs
+        elif nside_coverage_inputs > self.config.nside_coverage:
+            # Converting from higher resolution coverage to lower
+            # resolution coverage.
+            bit_shift = hsp.utils._compute_bitshift(self.config.nside_coverage,
+                                                    nside_coverage_inputs)
+            cov_pix = np.right_shift(cov_pix_inputs, bit_shift)
+        else:
+            # Converting from lower resolution coverage to higher
+            # resolution coverage.
+            bit_shift = hsp.utils._compute_bitshift(nside_coverage_inputs,
+                                                    self.config.nside_coverage)
+            cov_pix = np.left_shift(cov_pix_inputs, bit_shift)
 
         # Now read in each tract map and build the consolidated map.
         consolidated_map = None
         for tract_id in input_refs:
             input_map = input_refs[tract_id].get()
             if consolidated_map is None:
-                dtype = input_map.dtype
-                sentinel = input_map._sentinel
-                nside_coverage = input_map.nside_coverage
-                nside_sparse = input_map.nside_sparse
-                consolidated_map = hsp.HealSparseMap.make_empty(nside_coverage,
-                                                                nside_sparse,
-                                                                dtype,
-                                                                sentinel=sentinel)
-                consolidated_map._reserve_cov_pix(cov_pix)
+                consolidated_map = hsp.HealSparseMap.make_empty(
+                    self.config.nside_coverage,
+                    input_map.nside_sparse,
+                    input_map.dtype,
+                    sentinel=input_map._sentinel,
+                    cov_pixels=cov_pix)
 
             # Only use pixels that are properly inside the tract.
             vpix, ra, dec = input_map.valid_pixels_pos(return_pixels=True)