DM-4692: Refactor ProcessCcdTask and sub-tasks

.gitignore

@@ -12,6 +12,7 @@ config.log
 doc/*.tag
 doc/html
 doc/doxygen.conf
+doc/xml
 python/lsst/pipe/tasks/version.py
 tests/.tests
 ups/*.cfgc

examples/calibrateTask.py

@@ -27,18 +27,24 @@
 import numpy as np
 
 import lsst.utils
-import lsst.pipe.base              as pipeBase
-import lsst.daf.base               as dafBase
-import lsst.afw.coord              as afwCoord
-import lsst.afw.geom               as afwGeom
-import lsst.afw.table              as afwTable
-import lsst.afw.image              as afwImage
-import lsst.afw.display.ds9        as ds9
-from lsst.meas.astrom import AstrometryTask
+import lsst.pipe.base as pipeBase
+from lsst.daf.butlerUtils import ExposureIdInfo
+import lsst.afw.display as afwDisplay
+import lsst.afw.table as afwTable
+import lsst.afw.image as afwImage
+from lsst.pex.config import Config
+from lsst.afw.coord import Coord
+from lsst.afw.geom import PointD
+from lsst.meas.algorithms import LoadReferenceObjectsTask
+from lsst.meas.astrom import createMatchMetadata
+
+from lsst.pipe.tasks.characterizeImage import CharacterizeImageTask
 from lsst.pipe.tasks.calibrate import CalibrateTask
 
 np.random.seed(1)
 
+FilterName = "r"
+
 def loadData(pixelScale=1.0):
     """Prepare the data we need to run the example"""
 
@@ -52,81 +58,101 @@ def loadData(pixelScale=1.0):
     calib.setExptime(1.0)
     exposure.setCalib(calib)
     # add a filter
-    filterName = "r"
-    afwImage.Filter.define(afwImage.FilterProperty(filterName, 600, True))
-    exposure.setFilter(afwImage.Filter(filterName))
+    afwImage.Filter.define(afwImage.FilterProperty(FilterName, 600, True))
+    exposure.setFilter(afwImage.Filter(FilterName))
     # and a trivial WCS (needed by MyAstrometryTask)
-    pixelScale /= 3600.0                # degrees per pixel
-    wcs = afwImage.makeWcs(afwCoord.Coord(afwGeom.PointD(15, 1)), afwGeom.PointD(0, 0),
-                           pixelScale, 0.0, 0.0, pixelScale)
+    pixelScale /= 3600.0  # degrees per pixel
+    wcs = afwImage.makeWcs(Coord(PointD(15, 1)), PointD(0, 0), pixelScale, 0.0, 0.0, pixelScale)
     exposure.setWcs(wcs)
 
     return exposure
 
-class MyAstrometryTask(AstrometryTask):
-    """An override for CalibrateTask's astrometry task"""
-    def __init__(self, *args, **kwargs):
-        super(MyAstrometryTask, self).__init__(*args, **kwargs)
+class MyAstrometryTask(pipeBase.Task):
+    """An override for CalibrateTask's astrometry task that fakes a solution"""
+    ConfigClass = Config
+    _defaultName = "astrometry"
+
+    def __init__(self, schema=None, **kwargs):
+        pipeBase.Task(**kwargs)
 
     def run(self, exposure, sourceCat):
-        """My run method that totally fakes the astrometric solution"""
+        """Fake an astrometric solution
 
-        filterName = exposure.getFilter().getName()
-        wcs = exposure.getWcs()
-        #
-        # Fake a reference catalogue by copying fluxes from the list of Sources
-        #
-        schema = afwTable.SimpleTable.makeMinimalSchema()
-        schema.addField(afwTable.Field[float]("{}_flux".format(filterName), "Reference flux"))
-        schema.addField(afwTable.Field[float]("photometric", "I am a reference star"))
-        refCat = afwTable.SimpleCatalog(schema)
+        Pretend the current solution is perfect
+        and make a reference catalog that matches the source catalog
+        """
+        return self.loadAndMatch(exposure=exposure, sourceCat=sourceCat)
 
-        for s in sourceCat:
-            m = refCat.addNew()
-            flux = 1e-3*s.getPsfFlux()*np.random.normal(1.0, 2e-2)
-            m.set("{}_flux".format(filterName), flux)
-            m.setCoord(wcs.pixelToSky(s.getCentroid()))
-
-        refCat.get("photometric")[:] = True
-        #
-        # Perform the "match"
-        #
-        matches = []
-        md = dafBase.PropertyList()
-        for m, s in zip(refCat, sourceCat):
-            matches.append(afwTable.ReferenceMatch(m, s, 0.0))
+    def loadAndMatch(self, exposure, sourceCat):
+        """!Fake loading and matching
+
+        Copy the source catalog to a reference catalog and producing a match list
+        """
+        wcs = exposure.getWcs()
+        refSchema = LoadReferenceObjectsTask.makeMinimalSchema(
+            filterNameList = [FilterName],
+            addIsPhotometric = True,
+        )
+        refCat = afwTable.SimpleCatalog(refSchema)
+        refFluxKey = refSchema[FilterName + "_flux"].asKey()
+        refIsPhotoKey = refSchema["photometric"].asKey()
+
+        matches = lsst.afw.table.ReferenceMatchVector()
+        for src in sourceCat:
+            flux = 1e-3*src.getPsfFlux()*np.random.normal(1.0, 2e-2)
+            refObj = refCat.addNew()
+            refObj.set(refFluxKey, flux)
+            refObj.setCoord(wcs.pixelToSky(src.getCentroid()))
+            refObj.set(refIsPhotoKey, True)
+            match = lsst.afw.table.ReferenceMatch(refObj, src, 0)
+            matches.append(match)
 
         return pipeBase.Struct(
-            matches=matches,
-            matchMeta=md
-            )
+            refCat = refCat,
+            matches = matches,
+            matchMeta = createMatchMetadata(
+                bbox = exposure.getBBox(),
+                wcs = exposure.getWcs(),
+                filterName = FilterName,
+            ),
+        )
 
 def run(display=False):
     #
-    # Create the task
+    # Create the tasks
     #
+    charImageConfig = CharacterizeImageTask.ConfigClass()
+    charImageTask = CharacterizeImageTask(config=charImageConfig)
+
     config = CalibrateTask.ConfigClass()
-    config.initialPsf.pixelScale = 0.185 # arcsec per pixel
-    config.initialPsf.fwhm = 1.0
     config.astrometry.retarget(MyAstrometryTask)
     calibrateTask = CalibrateTask(config=config)
-    #
-    # Process the data
-    #
-    exposure = loadData(config.initialPsf.pixelScale)
-    result = calibrateTask.run(exposure)
 
-    exposure0, exposure = exposure, result.exposure
-    sources = result.sources
+    # load the data
+    # Exposure ID and the number of bits required for exposure IDs are usually obtained from a data repo,
+    # but here we pick reasonable values (there are 64 bits to share between exposure IDs and source IDs).
+    exposure = loadData()
+    exposureIdInfo = ExposureIdInfo(expId=1, expBits=5)
 
-    if display:                         # display on ds9 (see also --debug argparse option)
-        frame = 1
-        ds9.mtv(exposure, frame=frame)
+    # characterize the exposure to repair cosmic rays and fit a PSF model
+    # display now because CalibrateTask modifies the exposure in place
+    charRes = charImageTask.characterize(exposure=exposure, exposureIdInfo=exposureIdInfo)
+    if display:
+        displayFunc(charRes.exposure, charRes.sourceCat, frame=1)
 
-        with ds9.Buffering():
-            for s in sources:
-                xy = s.getCentroid()
-                ds9.dot('+', *xy, ctype=ds9.CYAN if s.get("flags_negative") else ds9.GREEN, frame=frame)
+    # calibrate the exposure
+    calRes = calibrateTask.calibrate(exposure=charRes.exposure, exposureIdInfo=exposureIdInfo)
+    if display:
+        displayFunc(calRes.exposure, calRes.sourceCat, frame=2)
+
+def displayFunc(exposure, sourceCat, frame):
+    display = afwDisplay.getDisplay(frame)
+    display.mtv(exposure)
+
+    with display.Buffering():
+        for s in sourceCat:
+            xy = s.getCentroid()
+            display.dot('+', *xy, ctype=afwDisplay.CYAN if s.get("flags_negative") else afwDisplay.GREEN)
 
 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 
@@ -135,7 +161,8 @@ def run(display=False):
     parser = argparse.ArgumentParser(description="Demonstrate the use of CalibrateTask")
 
     parser.add_argument('--debug', '-d', action="store_true", help="Load debug.py?", default=False)
-    parser.add_argument('--ds9', action="store_true", help="Display sources on ds9", default=False)
+    parser.add_argument('--display', action="store_true",
+        help="Display images in this example task (not using debug.py)", default=False)
 
     args = parser.parse_args()
 
@@ -145,4 +172,4 @@ def run(display=False):
         except ImportError as e:
             print >> sys.stderr, e
 
-    run(display=args.ds9)
+    run(display=args.display)

examples/detectAndMeasureExample.py

@@ -0,0 +1,108 @@
+#!/usr/bin/env python
+
+#
+# LSST Data Management System
+# Copyright 2008-2015 AURA/LSST.
+#
+# This product includes software developed by the
+# LSST Project (http://www.lsst.org/).
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the LSST License Statement and
+# the GNU General Public License along with this program.  If not,
+# see <http://www.lsstcorp.org/LegalNotices/>.
+#
+"""Example use of DetectAndMeasureTask
+"""
+
+import os
+import numpy as np
+
+import lsst.utils
+from lsst.daf.butlerUtils import ExposureIdInfo
+from lsst.afw.detection import GaussianPsf    
+import lsst.afw.image as afwImage
+from lsst.meas.astrom import displayAstrometry
+from lsst.meas.algorithms import estimateBackground
+from lsst.pipe.tasks.calibrate import DetectAndMeasureTask
+from lsst.pipe.tasks.repair import RepairTask
+
+np.random.seed(1)
+
+FilterName = "r"
+
+def loadData(psfSigma=1.5):
+    """Prepare the data we need to run the example"""
+
+    # Load sample input from disk
+    mypath = lsst.utils.getPackageDir('afwdata')
+    imFile = os.path.join(mypath, "CFHT", "D4", "cal-53535-i-797722_small_1.fits")
+
+    exposure = afwImage.ExposureF(imFile)
+    # add a filter
+    afwImage.Filter.define(afwImage.FilterProperty(FilterName, 600, True))
+    exposure.setFilter(afwImage.Filter(FilterName))
+    # add a simple Gaussian PSF model
+    psfModel = GaussianPsf(11, 11, psfSigma)
+    exposure.setPsf(psfModel)
+
+    return exposure
+
+
+def run(display=False):
+    """Subtract background, mask cosmic rays, then detect and measure
+    """
+    # Create the tasks; note that background estimation is performed by a function,
+    # not a task, though it has a config
+    repairConfig = RepairTask.ConfigClass()
+    repairTask = RepairTask(config=repairConfig)
+
+    backgroundConfig = estimateBackground.ConfigClass()
+
+    damConfig = DetectAndMeasureTask.ConfigClass()
+    damConfig.detection.thresholdValue = 5.0
+    damConfig.detection.includeThresholdMultiplier = 1.0
+    damConfig.measurement.doApplyApCorr = "yes"
+    detectAndMeasureTask = DetectAndMeasureTask(config=damConfig)
+
+    # load the data
+    # Exposure ID and the number of bits required for exposure IDs are usually obtained from a data repo,
+    # but here we pick reasonable values (there are 64 bits to share between exposure IDs and source IDs).
+    exposure = loadData()
+    exposureIdInfo = ExposureIdInfo(expId=1, expBits=5)
+
+    # repair cosmic rays
+    repairTask.run(exposure=exposure)
+
+    # subtract an initial estimate of background level
+    estBg, exposure = estimateBackground(
+        exposure = exposure,
+        backgroundConfig = backgroundConfig,
+        subtract = True,
+    )
+
+    # detect and measure
+    damRes = detectAndMeasureTask.run(exposure=exposure, exposureIdInfo=exposureIdInfo)
+    if display:
+        displayAstrometry(frame=2, exposure=damRes.exposure, sourceCat=damRes.sourceCat, pause=False)
+
+
+if __name__ == "__main__":
+    import argparse
+    parser = argparse.ArgumentParser(description="Demonstrate the use of DetectAndMeasureTask")
+
+    parser.add_argument('--display', action="store_true",
+        help="Display the output image and source catalog", default=False)
+
+    args = parser.parse_args()
+
+    run(display=args.display)