DM-30993: Implement reference catalog culling for astrometry in SFM

python/lsst/meas/algorithms/sourceSelector.py

@@ -562,6 +562,83 @@ def apply(self, catalog):
         return selected
 
 
+class CullFromMaskedRegion(pexConfig.Config):
+    """Deselect sources that lie in a "bad" mask plane.
+
+    This will select against objects whose image coordinates lie in a region
+    with any of the mask bits in the `badMaskNames` list set. Namely used for
+    a reference catalog for which the flag columns we would get from the
+    measurement plugins do not exist.
+
+    NOTE: In the context of reference objects, it is recommended NOT to include
+    EDGE in the `badMaskNames` list as that will remove all the reference objects
+    outside the detector but within the pixelMargin (thus nulling the pixelMargin
+    padding all together!)
+    """
+    badMaskNames = pexConfig.ListField(
+        dtype=str,
+        default=["NO_DATA", "NOT_DEBLENDED"],
+        doc="List of mask planes for which sources should be removed if a bit is set.",
+    )
+    xColName = pexConfig.Field(
+        dtype=str,
+        default="centroid_x",
+        doc="Name of column for image x coordinate."
+    )
+    yColName = pexConfig.Field(
+        dtype=str,
+        default="centroid_y",
+        doc="Name of column for image y coordinate."
+    )
+
+    def apply(self, catalog, exposure):
+        """Apply the mask plane requirements to a catalog.
+
+        Returns whether the sources were selected.
+
+        Parameters
+        ----------
+        catalog : `lsst.afw.table.SourceCatalog` or `pandas.DataFrame`
+                  or `astropy.table.Table`
+            Catalog of sources to which the requirements will be applied.
+        exposure : `lsst.afw.image.Exposure` or None
+            The exposure whose mask plane is to be respected.
+
+
+        Returns
+        -------
+        selected : `numpy.ndarray`
+            Boolean array indicating for each source whether it is selected
+            (True means selected).
+
+        Raises
+        ------
+        RuntimeError
+            Raised if exposure passed is `None`.
+        """
+        if exposure is None:
+            raise RuntimeError("Must provide an exposure to CullFromMaskedRegion selection.")
+        xRefList = catalog[self.xColName]
+        yRefList = catalog[self.yColName]
+        # Convert x, y coords to integers to map to indices in mask plane.
+        # If reference object nominally lies outside the exposure, consider
+        # it to be at the edge (and thus obeys those mask planes).
+        x0, y0 = exposure.getXY0()
+        xMax, yMax = exposure.getDimensions()
+        xRefList = [int(min(max(0, xRef - x0), xMax - 1)) for xRef in xRefList]
+        yRefList = [int(min(max(0, yRef - y0), yMax - 1)) for yRef in yRefList]
+        badMaskNames = []
+        maskPlaneDict = exposure.getMask().getMaskPlaneDict()
+        for badName in self.badMaskNames:
+            if badName in maskPlaneDict:
+                badMaskNames.append(badName)
+        bitmask = exposure.mask.getPlaneBitMask(badMaskNames)
+        toKeep = ((exposure.mask.array & bitmask) == 0)
+        selected = toKeep[yRefList, xRefList]  # x & y flipped for numpy arrays
+
+        return selected
+
+
 class ScienceSourceSelectorConfig(pexConfig.Config):
     """Configuration for selecting science sources"""
     doFluxLimit = pexConfig.Field(dtype=bool, default=False, doc="Apply flux limit?")
@@ -660,6 +737,8 @@ class ReferenceSourceSelectorConfig(pexConfig.Config):
     doMagError = pexConfig.Field(dtype=bool, default=False, doc="Apply magnitude error limit?")
     doRequireFiniteRaDec = pexConfig.Field(dtype=bool, default=True,
                                            doc="Apply finite sky coordinate check?")
+    doCullFromMaskedRegion = pexConfig.Field(dtype=bool, default=False,
+                                             doc="Apply image masked region culling?")
     magLimit = pexConfig.ConfigField(dtype=MagnitudeLimit, doc="Magnitude limit to apply")
     flags = pexConfig.ConfigField(dtype=RequireFlags, doc="Flags to require")
     unresolved = pexConfig.ConfigField(dtype=RequireUnresolved, doc="Star/galaxy separation to apply")
@@ -669,6 +748,8 @@ class ReferenceSourceSelectorConfig(pexConfig.Config):
     magError = pexConfig.ConfigField(dtype=MagnitudeErrorLimit, doc="Magnitude error limit to apply")
     colorLimits = pexConfig.ConfigDictField(keytype=str, itemtype=ColorLimit, default={},
                                             doc="Color limits to apply; key is used as a label only")
+    cullFromMaskedRegion = pexConfig.ConfigField(dtype=CullFromMaskedRegion,
+                                                 doc="Image mask plane criteria to apply")
 
 
 @pexConfig.registerConfigurable("references", sourceSelectorRegistry)
@@ -718,7 +799,8 @@ def selectSources(self, sourceCat, matches=None, exposure=None):
             selected &= self.config.requireFiniteRaDec.apply(sourceCat)
         for limit in self.config.colorLimits.values():
             selected &= limit.apply(sourceCat)
-
+        if self.config.doCullFromMaskedRegion:
+            selected &= self.config.cullFromMaskedRegion.apply(sourceCat, exposure)
         self.log.info("Selected %d/%d references", selected.sum(), len(sourceCat))
 
         return pipeBase.Struct(selected=selected)

tests/test_sourceSelector.py

@@ -26,7 +26,9 @@
 import astropy.units as u
 import warnings
 
+import lsst.afw.image
 import lsst.afw.table
+import lsst.geom
 import lsst.meas.algorithms
 import lsst.meas.base.tests
 import lsst.pipe.base
@@ -56,28 +58,35 @@ def setUp(self):
         schema.addField("nChild", np.int32, "Number of children")
         schema.addField("detect_isPrimary", "Flag", "Is primary detection?")
         schema.addField("sky_source", "Flag", "Empty sky region.")
+
+        self.xCol = "centroid_x"
+        self.yCol = "centroid_y"
+        schema.addField(self.xCol, float, "Centroid x value.")
+        schema.addField(self.yCol, float, "Centroid y value.")
+
         self.catalog = lsst.afw.table.SourceCatalog(schema)
         self.catalog.reserve(10)
         self.config = self.Task.ConfigClass()
+        self.exposure = None
 
     def tearDown(self):
         del self.catalog
 
     def check(self, expected):
         task = self.Task(config=self.config)
-        results = task.run(self.catalog)
+        results = task.run(self.catalog, exposure=self.exposure)
         self.assertListEqual(results.selected.tolist(), expected)
         self.assertListEqual([src.getId() for src in results.sourceCat],
                              [src.getId() for src, ok in zip(self.catalog, expected) if ok])
 
         # Check with pandas.DataFrame version of catalog
-        results = task.run(self.catalog.asAstropy().to_pandas())
+        results = task.run(self.catalog.asAstropy().to_pandas(), exposure=self.exposure)
         self.assertListEqual(results.selected.tolist(), expected)
         self.assertListEqual(list(results.sourceCat['id']),
                              [src.getId() for src, ok in zip(self.catalog, expected) if ok])
 
         # Check with astropy.table.Table version of catalog
-        results = task.run(self.catalog.asAstropy())
+        results = task.run(self.catalog.asAstropy(), exposure=self.exposure)
         self.assertListEqual(results.selected.tolist(), expected)
         self.assertListEqual(list(results.sourceCat['id']),
                              [src.getId() for src, ok in zip(self.catalog, expected) if ok])
@@ -369,6 +378,42 @@ def testFiniteRaDec(self):
 
         self.check([False, False, True, True, True])
 
+    def testCullFromMaskedRegion(self):
+        # Test that objects whose centroids land on specified mask(s) are
+        # culled.
+        maskNames = ["NO_DATA", "BLAH"]
+        num = 5
+        for _ in range(num):
+            self.catalog.addNew()
+
+        for x0, y0 in [[0, 0], [3, 8]]:
+            self.exposure = lsst.afw.image.ExposureF(5, 5)
+            self.exposure.setXY0(lsst.geom.Point2I(x0, y0))
+            mask = self.exposure.mask
+            for maskName in maskNames:
+                if maskName not in mask.getMaskPlaneDict():
+                    mask.addMaskPlane(maskName)
+            self.catalog[self.xCol][:] = x0 + 5.0
+            self.catalog[self.yCol][:] = y0 + 5.0
+            noDataPoints = [[0 + x0, 0 + y0], [3 + x0, 2 + y0]]
+            # Set first two entries in catalog to land in maskNames region.
+            for i, noDataPoint in enumerate(noDataPoints):
+                # Flip x & y for numpy array convention.
+                mask.array[noDataPoint[1] - y0][noDataPoint[0] - x0] = mask.getPlaneBitMask(
+                    maskNames[min(i, len(maskNames) - 1)]
+                )
+                self.catalog[self.xCol][i] = noDataPoint[0]
+                self.catalog[self.yCol][i] = noDataPoint[1]
+            self.config.doCullFromMaskedRegion = True
+            self.config.cullFromMaskedRegion.xColName = self.xCol
+            self.config.cullFromMaskedRegion.yColName = self.yCol
+            self.config.cullFromMaskedRegion.badMaskNames = maskNames
+            self.check([False, False, True, True, True])
+
+        # Reset config back to False and None for other tests.
+        self.config.doCullFromMaskedRegion = False
+        self.exposure = None
+
 
 class TestBaseSourceSelector(lsst.utils.tests.TestCase):
     """Test the API of the Abstract Base Class with a trivial example."""