Review for DM-4831: Add bright object masks to pipeline outputs

python/lsst/pipe/tasks/assembleCoadd.py

@@ -22,6 +22,8 @@
 #
 import numpy
 import lsst.pex.config as pexConfig
+import lsst.pex.exceptions as pexExceptions
+import lsst.afw.detection as afwDetect
 import lsst.afw.geom as afwGeom
 import lsst.afw.image as afwImage
 import lsst.afw.math as afwMath
@@ -118,6 +120,18 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass):
     )
     removeMaskPlanes = pexConfig.ListField(dtype=str, default=["CROSSTALK", "NOT_DEBLENDED"],\
                                  doc="Mask planes to remove before coadding")
+    #
+    # N.b. These configuration options only set the bitplane config.brightObjectMaskName
+    # To make this useful you *must* also configure the flags.pixel algorithm, for example
+    # by adding 
+    #   config.measurement.plugins["base_PixelFlags"].masksFpCenter.append("BRIGHT_OBJECT")
+    #   config.measurement.plugins["base_PixelFlags"].masksFpAnywhere.append("BRIGHT_OBJECT")
+    # to your measureCoaddSources.py config overrides
+    #
+    doMaskBrightObjects = pexConfig.Field(dtype=bool, default=False,
+                                          doc="Set mask and flag bits for bright objects?")
+    brightObjectMaskName = pexConfig.Field(dtype=str, default="BRIGHT_OBJECT",
+                                           doc="Name of mask bit used for bright objects")
 
     def setDefaults(self):
         CoaddBaseTask.ConfigClass.setDefaults(self)
@@ -136,6 +150,15 @@ def __init__(self, *args, **kwargs):
         self.makeSubtask("matchBackgrounds")
         self.makeSubtask("scaleZeroPoint")
 
+        if self.config.doMaskBrightObjects:
+            mask = afwImage.MaskU()
+            try:
+                self.brightObjectBitmask = 1 << mask.addMaskPlane(self.config.brightObjectMaskName)
+            except pexExceptions.LsstCppException:
+                raise RuntimeError("Unable to define mask plane for bright objects; planes used are %s" %
+                                   mask.getMaskPlaneDict().keys())
+            del mask
+
     @pipeBase.timeMethod
     def run(self, dataRef, selectDataList=[]):
         """Assemble a coadd from a set of coaddTempExp
@@ -200,6 +223,10 @@ def run(self, dataRef, selectDataList=[]):
             varArray = coaddExp.getMaskedImage().getVariance().getArray()
             varArray[:] = numpy.where(varArray > 0, varArray, numpy.inf)
 
+        if self.config.doMaskBrightObjects:
+            brightObjectMasks = self.readBrightObjectMasks(dataRef)
+            self.setBrightObjectMasks(coaddExp, dataRef.dataId, brightObjectMasks)
+
         if self.config.doWrite:
             self.writeCoaddOutput(dataRef, coaddExp)
 
@@ -548,6 +575,48 @@ def addBackgroundMatchingMetadata(self, coaddExposure, tempExpRefList, backgroun
                                    backgroundInfo.matchedMSE/backgroundInfo.diffImVar)
                 metadata.addDouble("CTExp_SDQA2_%d" % (ind),
                                    backgroundInfo.fitRMS)
+
+    def readBrightObjectMasks(self, dataRef):
+        """Returns None on failure"""
+        try:
+            return dataRef.get("brightObjectMask", immediate=True)
+        except Exception as e:
+            self.log.warn("Unable to read brightObjectMask for %s: %s" % (dataRef.dataId, e))
+            return None
+
+    def setBrightObjectMasks(self, exposure, dataId, brightObjectMasks):
+        """Set the bright object masks
+
+        exposure:          Exposure under consideration
+        dataId:            Data identifier dict for patch
+        brightObjectMasks: afwTable of bright objects to mask
+        """
+        #
+        # Check the metadata specifying the tract/patch/filter
+        #
+        if brightObjectMasks is None:
+            self.log.warn("Unable to apply bright object mask: none supplied")
+            return
+        self.log.info("Applying %d bright object masks to %s" % (len(brightObjectMasks), dataId))
+        md = brightObjectMasks.table.getMetadata()
+        for k in dataId:
+            if not md.exists(k):
+                self.log.warn("Expected to see %s in metadata" % k)
+            else:
+                if md.get(k) != dataId[k]:
+                    self.log.warn("Expected to see %s == %s in metadata, saw %s" % (k, md.get(k), dataId[k]))
+
+        mask = exposure.getMaskedImage().getMask()
+        wcs = exposure.getWcs()
+        plateScale = wcs.pixelScale().asArcseconds()
+
+        for rec in brightObjectMasks:
+            center = afwGeom.PointI(wcs.skyToPixel(rec.getCoord()))
+            radius = rec["radius"].asArcseconds()/plateScale   # convert to pixels
+
+            foot = afwDetect.Footprint(center, radius, exposure.getBBox())
+            afwDetect.setMaskFromFootprint(mask, foot, self.brightObjectBitmask)
+
     @classmethod
     def _makeArgumentParser(cls):
         """Create an argument parser

python/lsst/pipe/tasks/objectMasks.py

@@ -0,0 +1,151 @@
+import re
+import lsst.daf.base as dafBase
+import lsst.pex.logging as pexLog
+import lsst.afw.coord as afwCoord
+import lsst.afw.geom as afwGeom
+import lsst.afw.table as afwTable
+
+class ObjectMaskCatalog(object):
+    """Class to support bright object masks
+
+    N.b. I/O is done by providing a readFits method which fools the butler.
+    """
+    def __init__(self):
+        schema = afwTable.SimpleTable.makeMinimalSchema()
+        schema.addField("radius", "Angle", "radius of mask")
+
+        self._catalog = afwTable.SimpleCatalog(schema)
+        self._catalog.table.setMetadata(dafBase.PropertyList())
+
+        self.table  = self._catalog.table
+        self.addNew = self._catalog.addNew
+
+    def __len__(self):
+        return len(self._catalog)
+
+    def __iter__(self):
+        return iter(self._catalog)
+
+    def __getitem__(self, i):
+        return self._catalog.__getitem__(i)
+
+    def __setitem__(self, i, v):
+        return self._catalog.__setitem__(i, v)
+
+    @staticmethod
+    def readFits(fileName, hdu=0, flags=0):
+        """Read a ds9 region file, returning a ObjectMaskCatalog object
+
+        This method is called "readFits" to fool the butler. The corresponding mapper entry looks like
+        brightObjectMask: {
+            template:      "deepCoadd/BrightObjectMasks/%(tract)d/BrightObjectMask-%(tract)d-%(patch)s-%(filter)s.reg"
+            python:        "lsst.obs.subaru.objectMasks.ObjectMaskCatalog"
+            persistable:   "PurePythonClass"
+            storage:       "FitsCatalogStorage"
+        }
+        and this is the only way I know to get it to read a random file type, in this case a ds9 region file
+
+        This method expects to find files named as BrightObjectMask-%(tract)d-%(patch)s-%(filter)s.reg
+        The files should be structured as follows:
+
+        # Description of catalogue as a comment
+        # CATALOG: catalog-id-string
+        # TRACT: 0
+        # PATCH: 5,4
+        # FILTER: HSC-I
+
+        wcs; fk5
+
+        circle(RA, DEC, RADIUS) # ID: 1
+
+        The commented lines must be present, with the relevant fields such as tract patch and filter filled
+        in. The coordinate system must be listed as above. Each patch is specified as a circle, with an RA,
+        DEC, and Radius specified in decimal degrees. Only circles are supported as region definitions
+        currently.
+        """
+
+        log = pexLog.getDefaultLog().createChildLog("ObjectMaskCatalog")
+
+        brightObjects = ObjectMaskCatalog()
+        checkedWcsIsFk5 = False
+
+        with open(fileName) as fd:
+            for lineNo, line in enumerate(fd.readlines(), 1):
+                line = line.rstrip()
+
+                if re.search(r"^\s*#", line):
+                    #
+                    # Parse any line of the form "# key : value" and put them into the metadata.
+                    #
+                    # The medatdata values must be defined as outlined in the above docstring
+                    # 
+                    # The value of these three keys will be checked,
+                    # so get them right!
+                    #
+                    mat = re.search(r"^\s*#\s*([a-zA-Z][a-zA-Z0-9_]+)\s*:\s*(.*)", line)
+                    if mat:
+                        key, value = mat.group(1).lower(), mat.group(2)
+                        if key == "tract":
+                            value = int(value)
+
+                        brightObjects.table.getMetadata().set(key, value)
+
+                line = re.sub(r"^\s*#.*", "", line)
+                if not line:
+                    continue
+
+                if re.search(r"^\s*wcs\s*;\s*fk5\s*$", line, re.IGNORECASE):
+                    checkedWcsIsFk5 = True
+                    continue
+
+                # This regular expression parses the regions file for each region to be masked,
+                # with the format as specified in the above docstring.
+                mat = re.search(r"^\s*circle(?:\s+|\s*\(\s*)"
+                                "(\d+(?:\.\d*)([d]*))" "(?:\s+|\s*,\s*)"
+                                "([+-]?\d+(?:\.\d*)([d]*))" "(?:\s+|\s*,\s*)"
+                                "(\d+(?:\.\d*))([d'\"]*)" "(?:\s*|\s*\)\s*)"
+                                "\s*#\s*ID:\s*(\d+)" "\s*$"
+                                , line)
+                if mat:
+                    ra, raUnit, dec, decUnit, radius, radiusUnit, _id = mat.groups()
+
+                    _id = int(_id)
+                    ra =     convertToAngle(ra, raUnit, "ra", fileName, lineNo)
+                    dec =    convertToAngle(dec, decUnit, "dec", fileName, lineNo)
+                    radius = convertToAngle(radius, radiusUnit, "radius", fileName, lineNo)
+
+                    rec = brightObjects.addNew()
+                    # N.b. rec["coord"] = Coord is not supported, so we have to use the setter
+                    rec["id"] = _id
+                    rec.setCoord(afwCoord.Fk5Coord(ra, dec))
+                    rec["radius"] = radius
+                else:
+                    log.warn("Unexpected line \"%s\" at %s:%d" % (line, fileName, lineNo))
+
+        if not checkedWcsIsFk5:
+            raise RuntimeError("Expected to see a line specifying an fk5 wcs")
+
+        # This makes the deep copy contiguous in memory so that a ColumnView can be exposed to Numpy
+        brightObjects._catalog = brightObjects._catalog.copy(True)
+
+        return brightObjects
+
+
+def convertToAngle(var, varUnit, what, fileName, lineNo):
+    """Given a variable and its units, return an afwGeom.Angle
+
+    what, fileName, and lineNo are used to generate helpful error messages
+    """
+    var = float(var)
+
+    if varUnit in ("d", ""):
+        pass
+    elif varUnit == "'":
+        var /= 60.0
+    elif varUnit == '"':
+        var /= 3600.0
+    else:
+        raise RuntimeError("unsupported unit \"%s\" for %s at %s:%d" %
+                          (varUnit, what, fileName, lineNo))
+
+    return var*afwGeom.degrees