DM-3766: Add Butler access to calibration data in obs_decam

config/ingestCalibs.py

@@ -0,0 +1,26 @@
+from lsst.obs.decam.ingestCalibs import DecamCalibsParseTask
+config.parse.retarget(DecamCalibsParseTask)
+config.parse.hdu = 1
+# N30 is not included becasue it is not functional.
+config.parse.extnames = ['S1', 'S2', 'S3', 'S4', 'S5', 'S6', 'S7', 'S8', 'S9', 'S10', 'S11', 'S12', 'S13',
+                       'S14', 'S15', 'S16', 'S17', 'S18', 'S19', 'S20', 'S21', 'S22', 'S23', 'S24', 'S25',
+                       'S26', 'S27', 'S28', 'S29', 'S30', 'S31', 'N1', 'N2', 'N3', 'N4', 'N5', 'N6', 'N7',
+                       'N8', 'N9', 'N10', 'N11', 'N12', 'N13', 'N14', 'N15', 'N16', 'N17', 'N18', 'N19',
+                       'N20', 'N21', 'N22', 'N23', 'N24', 'N25', 'N26', 'N27', 'N28', 'N29', 'N31',
+                       ]
+config.parse.translators = {'filter': 'translate_filter',
+                          'ccdnum': 'translate_ccdnum',
+                          'calibDate': 'translate_date',
+                          'validStart': 'translate_date',
+                          'validEnd': 'translate_date',
+                          }
+config.register.columns = {'filter': 'text',
+                         'ccdnum': 'int',
+                         'path': 'text',
+                         'calibDate': 'text',
+                         'validStart': 'text',
+                         'validEnd': 'text',
+                         }
+config.register.unique = ['filter', 'ccdnum', 'calibDate']
+config.register.tables = ['bias', 'flat', 'fringe']
+config.register.visit = ['calibDate']

policy/DecamMapper.paf

@@ -1,6 +1,6 @@
 #<?cfg paf policy ?>
 
-needCalibRegistry: false
+needCalibRegistry: true
 
 levels: {
     # Keys that are NOT relevant for a particular level
@@ -71,6 +71,51 @@ exposures: {
     }
 }
 
+calibrations: {
+    bias: {
+        # Given that MasterCal hdu number is identical to ccdnum
+        template:    "%(path)s[%(ccdnum)d]"
+        python:      "lsst.afw.image.DecoratedImageF"
+        persistable: "DecoratedImageF"
+        storage:     "FitsStorage"
+        level:       "ccd"
+        tables:      "bias"
+        columns:     "date"
+        reference:   "raw_visit"
+        refCols:     "visit"
+        filter:      false
+        validRange:  true
+        obsTimeName: date
+    }
+    flat: {
+        template:    "%(path)s[%(ccdnum)d]"
+        python:      "lsst.afw.image.DecoratedImageF"
+        persistable: "DecoratedImageF"
+        storage:     "FitsStorage"
+        level:       "ccd"
+        tables:      "flat"
+        columns:     "filter" "date"
+        reference:   "raw_visit"
+        refCols:     "visit"
+        filter:      true
+        validRange:  true
+        obsTimeName: date
+    }
+    fringe: {
+        template:    "%(path)s"
+        python:      "lsst.afw.image.ImageF"
+        persistable: "ImageF"
+        storage:     "FitsStorage"
+        level:       "ccd"
+        tables:      "fringe"
+        columns:     "filter"
+        reference:   "raw_visit"
+        refCols:     "visit"
+        filter:      true
+        validRange:  false
+    }
+}
+
 datasets: {
 
     ccdExposureId: {

python/lsst/obs/decam/decamMapper.py

@@ -19,13 +19,12 @@
 # the GNU General Public License along with this program.  If not,
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
-import os, re
+import re
 import numpy as np
 import lsst.afw.image as afwImage
 import lsst.afw.image.utils as afwImageUtils
-from lsst.daf.butlerUtils import CameraMapper, exposureFromImage
+from lsst.daf.butlerUtils import CameraMapper
 import lsst.pex.policy as pexPolicy
-import lsst.daf.base as dafBase
 
 np.seterr(divide="ignore")
 
@@ -46,8 +45,8 @@ def __init__(self, inputPolicy=None, **kwargs):
         afwImageUtils.defineFilter('y', lambdaEff=1000, alias=['Y DECam c0005 10095.0 1130.0'])
 
     def _extractDetectorName(self, dataId):
-        nameTuple = self.registry.executeQuery(['side','ccd'], ['raw',], [('ccdnum','?'), ('visit','?')], None,
-                                   (dataId['ccdnum'], dataId['visit']))
+        nameTuple = self.registry.executeQuery(['side','ccd'], ['raw',], [('ccdnum','?'), ('visit','?')],
+                                               None, (dataId['ccdnum'], dataId['visit']))
         if len(nameTuple) > 1:
             raise RuntimeError("More than one name returned")
         if len(nameTuple) == 0:
@@ -134,3 +133,40 @@ def bypass_instcal(self, datasetType, pythonType, butlerLocation, dataId):
 
         exp.setMetadata(md) # Do we need to remove WCS/calib info?
         return exp
+
+    def _standardizeMasterCal(self, datasetType, item, dataId, setFilter=False):
+        """Standardize a MasterCal image obtained from NOAO archive into Exposure
+
+        These MasterCal images are MEF files with one HDU for each detector.
+        Some WCS header, eg CTYPE1, exists only in the zeroth extensionr,
+        so info in the zeroth header need to be copied over to metadata.
+
+        @param datasetType: Dataset type ("bias" or "flat")
+        @param item: The image read by the butler
+        @param dataId: Data identifier
+        @param setFilter: Whether to set the filter in the Exposure
+        @return (lsst.afw.image.Exposure) the standardized Exposure
+        """
+        mi = afwImage.makeMaskedImage(item.getImage())
+        md = item.getMetadata()
+        masterCalMap = getattr(self, "map_" + datasetType)
+        masterCalPath = masterCalMap(dataId).getLocations()[0]
+        headerPath = re.sub(r'[\[](\d+)[\]]$', "[0]", masterCalPath)
+        md0 = afwImage.readMetadata(headerPath)
+        for kw in md0.paramNames():
+            if kw not in md.paramNames():
+                md.add(kw, md0.get(kw))
+        wcs = afwImage.makeWcs(md, True)
+        exp = afwImage.makeExposure(mi, wcs)
+        exp.setMetadata(md)
+        return self._standardizeExposure(self.calibrations[datasetType], exp, dataId, filter=setFilter)
+
+    def std_bias(self, item, dataId):
+        return self._standardizeMasterCal("bias", item, dataId, setFilter=False)
+
+    def std_flat(self, item, dataId):
+        return self._standardizeMasterCal("flat", item, dataId, setFilter=True)
+
+    def std_fringe(self, item, dataId):
+        exp = afwImage.makeExposure(afwImage.makeMaskedImage(item))
+        return self._standardizeExposure(self.calibrations["fringe"], exp, dataId)

python/lsst/obs/decam/ingest.py

@@ -37,8 +37,15 @@ def parseExtname(md):
 class DecamIngestArgumentParser(IngestArgumentParser):
     def __init__(self, *args, **kwargs):
         super(DecamIngestArgumentParser, self).__init__(*args, **kwargs)
-        self.add_argument("--filetype", default="instcal", choices=["instcal", "raw"], help="Data processing level of the files to be ingested")
-        self.description = "To ingest instcal data, the following directory structure is expected:\n    dqmask/ instcal/ wtmap/ \nThe science pixels, mask, and weight (inverse variance) are stored in separate files each with a unique name but with a common unique identifier EXPNUM in the FITS header. The 3 files of the same EXPNUM will be aggregated. For example, the user creates the registry by running \n    ingestImagesDecam.py outputRepository --mode=link instcal/*fits"
+        self.add_argument("--filetype", default="instcal", choices=["instcal", "raw"],
+                          help="Data processing level of the files to be ingested")
+        self.description = "To ingest instcal data, the following directory structure is expected:"\
+            "\n    dqmask/ instcal/ wtmap/"\
+            "\nThe science pixels, mask, and weight (inverse variance) are stored in"\
+            "\nseparate files each with a unique name but with a common unique identifier"\
+            "\nEXPNUM in the FITS header. The 3 files of the same EXPNUM will be aggregated."\
+            "\nFor example, the user creates the registry by running"\
+            "\n    ingestImagesDecam.py outputRepository --mode=link instcal/*fits"
 
 
 class DecamIngestTask(IngestTask):
@@ -49,17 +56,23 @@ def __init__(self, *args, **kwargs):
     def run(self, args):
         """Ingest all specified files and add them to the registry"""
         if args.filetype == "instcal":
-            with self.register.openRegistry(args.butler, create=args.create) if not args.dryrun else None as registry:
+            with self.register.openRegistry(args.butler, create=args.create, dryrun=args.dryrun) as registry:
                 for infile in args.files:
                     fileInfo, hduInfoList = self.parse.getInfo(infile, args.filetype)
                     if len(hduInfoList) > 0:
-                        outfileInstcal = os.path.join(args.butler, self.parse.getDestination(args.butler, hduInfoList[0], infile, "instcal"))
-                        outfileDqmask = os.path.join(args.butler, self.parse.getDestination(args.butler, hduInfoList[0], infile, "dqmask"))
-                        outfileWtmap = os.path.join(args.butler, self.parse.getDestination(args.butler, hduInfoList[0], infile, "wtmap"))
-
-                        ingestedInstcal = self.ingest(fileInfo["instcal"], outfileInstcal, mode=args.mode, dryrun=args.dryrun)
-                        ingestedDqmask = self.ingest(fileInfo["dqmask"], outfileDqmask, mode=args.mode, dryrun=args.dryrun)
-                        ingestedWtmap = self.ingest(fileInfo["wtmap"], outfileWtmap, mode=args.mode, dryrun=args.dryrun)
+                        outfileInstcal = os.path.join(args.butler, self.parse.getDestination(args.butler,
+                                                      hduInfoList[0], infile, "instcal"))
+                        outfileDqmask = os.path.join(args.butler, self.parse.getDestination(args.butler,
+                                                     hduInfoList[0], infile, "dqmask"))
+                        outfileWtmap = os.path.join(args.butler, self.parse.getDestination(args.butler,
+                                                    hduInfoList[0], infile, "wtmap"))
+
+                        ingestedInstcal = self.ingest(fileInfo["instcal"], outfileInstcal,
+                                                      mode=args.mode, dryrun=args.dryrun)
+                        ingestedDqmask = self.ingest(fileInfo["dqmask"], outfileDqmask,
+                                                     mode=args.mode, dryrun=args.dryrun)
+                        ingestedWtmap = self.ingest(fileInfo["wtmap"], outfileWtmap,
+                                                    mode=args.mode, dryrun=args.dryrun)
 
                         if not (ingestedInstcal or ingestedDqmask or ingestedWtmap):
                             continue
@@ -74,7 +87,8 @@ def run(self, args):
                 for infile in args.files:
                     fileInfo, hduInfoList = self.parse.getInfo(infile, args.filetype)
                     fileInfo['hdu'] = 0
-                    outfileRaw = super(DecamParseTask, self.parse).getDestination(args.butler, fileInfo, infile)
+                    outfileRaw = super(DecamParseTask, self.parse).getDestination(args.butler,
+                                                                                  fileInfo, infile)
                     self.ingest(infile, outfileRaw, mode=args.mode, dryrun=args.dryrun)
                     for info in hduInfoList:
                         self.register.addRow(registry, info, dryrun=args.dryrun, create=args.create)

python/lsst/obs/decam/ingestCalibs.py

@@ -0,0 +1,63 @@
+import collections
+import re
+from lsst.pipe.tasks.ingestCalibs import CalibsParseTask
+
+class DecamCalibsParseTask(CalibsParseTask):
+    def getInfo(self, filename):
+        """Get information about the image from the filename and/or its contents
+
+        @param filename: Name of file to inspect
+        @return File properties; list of file properties for each extension
+        """
+        phuInfo, infoList = CalibsParseTask.getInfo(self, filename)
+        # Single-extension fits without EXTNAME can be a valid CP calibration product
+        # Use info of primary header unit
+        if not infoList:
+            infoList.append(phuInfo)
+        return phuInfo, infoList
+
+    def translate_ccdnum(self, md):
+        """Return CCDNUM as a integer
+
+        @param md (PropertySet) FITS header metadata
+        """
+        if md.exists("CCDNUM"):
+            ccdnum = md.get("CCDNUM")
+        else:
+            self.log.warn("Unable to find value for CCDNUM")
+            ccdnum = None
+        # Some MasterCal from NOAO Archive has 2 CCDNUM keys in each HDU
+        # Make sure only one integer is returned.
+        if isinstance(ccdnum, collections.Sequence):
+            try:
+                ccdnum = ccdnum[0]
+            except IndexError:
+                ccdnum = None
+        return ccdnum
+
+    def translate_date(self, md):
+        """Extract the date as a strong in format YYYY-MM-DD from the FITS header DATE-OBS.
+        Return "unknown" if the value cannot be found or converted.
+
+        @param md (PropertySet) FITS header metadata
+        """
+        if md.exists("DATE-OBS"):
+            date = md.get("DATE-OBS")
+            found = re.search('(\d\d\d\d-\d\d-\d\d)', date)
+            if found:
+                date = found.group(1)
+            else:
+                self.log.warn("DATE-OBS does not match format YYYY-MM-DD")
+                date = "unknown"
+        else:
+            self.log.warn("Unable to find value for DATE-OBS")
+            date = "unknown"
+        return date
+
+    @staticmethod
+    def getExtensionName(md):
+        """ Get the name of the extension
+
+        @param md (PropertySet) FITS header metadata
+        """
+        return md.get('EXTNAME')

tests/getRaw.py

@@ -1,5 +1,5 @@
 #!/usr/bin/env python
-
+from __future__ import print_function
 #
 # LSST Data Management System
 # Copyright 2012, 2015 LSST Corporation.
@@ -31,10 +31,7 @@
 
 import lsst.pex.exceptions as pexExcept
 import lsst.daf.persistence as dafPersist
-from lsst.obs.decam import DecamMapper
 
-import lsst.afw.cameraGeom as cameraGeom
-import lsst.afw.cameraGeom.utils as cameraGeomUtils
 
 class GetRawTestCase(unittest.TestCase):
     """Testing butler raw image retrieval"""
@@ -46,10 +43,11 @@ def setUp(self):
             message = "testdata_decam not setup. Skipping."
             warnings.warn(message)
             raise unittest.SkipTest(message)
-        self.butler = dafPersist.Butler(root=os.path.join(datadir, "rawData"))
+        self.butler = dafPersist.Butler(root=os.path.join(datadir, "rawData"),
+                                        calibRoot=os.path.join(datadir, "calib"))
         self.size = (2160, 4146)
-        self.dataId = {'visit': 237628, 'ccdnum': 10}
-        self.filter = "i"
+        self.dataId = {'visit': 229388, 'ccdnum': 13}
+        self.filter = "z"
 
     def tearDown(self):
         del self.butler
@@ -61,10 +59,10 @@ def testRaw(self):
         """Test retrieval of raw image"""
         exp = self.butler.get("raw", self.dataId)
 
-        print "dataId: ", self.dataId
-        print "width: ", exp.getWidth()
-        print "height: ", exp.getHeight()
-        print "detector id: ", exp.getDetector().getId()
+        print("dataId: %s" % self.dataId)
+        print("width: %s" % exp.getWidth())
+        print("height: %s" % exp.getHeight())
+        print("detector id: %s" % exp.getDetector().getId())
 
         self.assertEqual(exp.getWidth(), self.size[0])
         self.assertEqual(exp.getHeight(), self.size[1])
@@ -74,11 +72,36 @@ def testRaw(self):
     def testRawMetadata(self):
         """Test retrieval of metadata"""
         md = self.butler.get("raw_md", self.dataId)
-        print "EXPNUM(visit): ",md.get('EXPNUM')
-        print "ccdnum:", md.get('CCDNUM')
+        print("EXPNUM(visit): %s" % md.get('EXPNUM'))
+        print("ccdnum: %s" % md.get('CCDNUM'))
         self.assertEqual(md.get('EXPNUM'), self.dataId["visit"])
         self.assertEqual(md.get('CCDNUM'), self.dataId["ccdnum"])
 
+    def testBias(self):
+        """Test retrieval of bias image"""
+        exp = self.butler.get("bias", self.dataId)
+        print("dataId: %s" % self.dataId)
+        print("detector id: %s" % exp.getDetector().getId())
+        self.assertEqual(exp.getDetector().getId(), self.dataId["ccdnum"])
+
+    def testFlat(self):
+        """Test retrieval of flat image"""
+        exp = self.butler.get("flat", self.dataId)
+        print("dataId: %s" % self.dataId)
+        print("detector id: %s" % exp.getDetector().getId())
+        print("filter: %s" % self.filter)
+        self.assertEqual(exp.getDetector().getId(), self.dataId["ccdnum"])
+        self.assertEqual(exp.getFilter().getFilterProperty().getName(), self.filter)
+
+    def testFringe(self):
+        """Test retrieval of fringe image"""
+        exp = self.butler.get("fringe", self.dataId)
+        print("dataId: %s" % self.dataId)
+        print("detector id: %s" % exp.getDetector().getId())
+        print("filter: %s" % self.filter)
+        self.assertEqual(exp.getDetector().getId(), self.dataId["ccdnum"])
+        self.assertEqual(exp.getFilter().getFilterProperty().getName(), self.filter)
+
 #-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-
 
 def suite():