Merge pull request #16 from erykoff/lsstsims2

Add ability to mask bad LSSTCam amps.

decasu/configuration.py

@@ -56,6 +56,7 @@ class Configuration(object):
     border: int = 15
     amp_boundary: int = 1024
     use_two_amps: bool = True
+    mask_lsstcam_bad_amps: bool = False
     use_wcs: bool = True
     ra_corner_fields: List[str] = field(default_factory=_default_ra_corner_fields)
     dec_corner_fields: List[str] = field(default_factory=_default_dec_corner_fields)
@@ -90,6 +91,9 @@ def _validate(self):
             except ImportError:
                 raise RuntimeError("Cannot use lsst db without Rubin Science Pipelines setup.")
 
+        if self.use_two_amps and self.mask_lsstcam_bad_amps:
+            raise RuntimeError("Cannot set both use_two_amps and mask_lsstcam_bad_amps.")
+
     @classmethod
     def load_yaml(cls, configfile):
         """

decasu/lsst_wcs_db.py

@@ -172,7 +172,6 @@ def __call__(self, row):
 
         # Link to global table.
         self.table = decasu_globals.table
-        # self.camera = decasu_globals.lsst_camera
 
         detector = decasu_globals.lsst_camera[self.table['detector'][row]]
 

decasu/region_mapper.py

@@ -4,6 +4,7 @@
 import hpgeom as hpg
 import esutil
 import time
+from functools import lru_cache
 
 import coord
 
@@ -555,6 +556,7 @@ def build_region_input_map(self, indices, tilename=None, hpix=None):
                     pixels = poly.get_pixels(nside=self.config.nside)
 
             # Check for bad amps -- only in two amp mode
+            mask_pixel_arrays = []
             if self.config.use_two_amps:
                 if int(dg.table[self.config.ccd_field][ind]) in list(self.config.bad_amps):
                     ba = self.config.bad_amps[int(dg.table[self.config.ccd_field][ind])]
@@ -563,6 +565,21 @@ def build_region_input_map(self, indices, tilename=None, hpix=None):
                             pixels_a = np.array([], dtype=np.int64)
                         elif b.lower() == 'b':
                             pixels_b = np.array([], dtype=np.int64)
+            elif self.config.mask_lsstcam_bad_amps:
+                det = int(dg.table[self.config.ccd_field][ind])
+                if det in list(self.config.bad_amps):
+                    ba = self.config.bad_amps[det]
+                    for b in ba:
+                        x_coords, y_coords = _get_numpy_corners_from_detector(
+                            dg.lsst_camera[det],
+                            b,
+                        )
+                        ra, dec = wcs.pixelToSkyArray(x_coords.astype(np.float64),
+                                                      y_coords.astype(np.float64),
+                                                      degrees=True)
+                        maskpoly = healsparse.Polygon(ra=ra, dec=dec, value=1)
+
+                        mask_pixel_arrays.append(maskpoly.get_pixels(nside=self.config.nside))
 
             if int(dg.table[self.config.ccd_field][ind]) in self.config.bad_ccds:
                 if self.config.use_two_amps:
@@ -607,6 +624,10 @@ def build_region_input_map(self, indices, tilename=None, hpix=None):
                     ok = ((pixels >= pixel_min) & (pixels <= pixel_max))
                     region_input_map.set_bits_pix(pixels[ok], [bit])
 
+            if len(mask_pixel_arrays) > 0:
+                mask_pixels = np.concatenate(mask_pixel_arrays)
+                region_input_map.clear_bits_pix(mask_pixels, [bit])
+
         region_input_map.metadata = metadata
 
         return region_input_map
@@ -762,3 +783,19 @@ def _get_maskcircle_from_row(self, table_row):
                                        radius=table_row['radius']/3600.,
                                        value=1)
         return maskcircle
+
+
+@lru_cache(maxsize=200)
+def _get_numpy_corners_from_detector(detector, amp_name):
+    bbox = detector[amp_name].getBBox()
+
+    x_coords = np.array([bbox.getBeginX(),
+                         bbox.getBeginX(),
+                         bbox.getEndX(),
+                         bbox.getEndX()])
+    y_coords = np.array([bbox.getBeginY(),
+                         bbox.getEndY(),
+                         bbox.getEndY(),
+                         bbox.getBeginY()])
+
+    return x_coords, y_coords

tests/configs/config_hpix_lsst.yaml

@@ -6,13 +6,14 @@ map_types:
   exptime: ['sum']
   airmass: ['wmean', 'min', 'max']
 use_two_amps: False
+mask_lsstcam_bad_amps: True
 arcsec_per_pix: 0.20
 zp_global: 27.0
 exp_field: "observationID"
 ccd_field: "detector"
 band_field: "filter"
 mjd_field: "mjd"
-bad_amps: {}
+bad_amps: {110: ["C10", "C11"]}
 bad_ccds: []
 extra_fields: {}
 latitude: -30.244633333333333