Tickets/dm 5503

.gitignore

@@ -18,6 +18,7 @@ tests/.tests
 tests/data/postISRCCD/
 tests/data/snapExp/
 version.py
+bin/extractPhosimGainSaturation.py
 bin/genDefectRegistry.py
 bin/genInputRegistry.py
 bin/ingestSimImages.py

bin.src/genInputRegistry.py

@@ -36,6 +36,7 @@
 import lsst.afw.image as afwImage
 import lsst.skypix as skypix
 
+
 def process(dirList, inputRegistry, outputRegistry="registry.sqlite3"):
     if os.path.exists(outputRegistry):
         print >>sys.stderr, "Output registry exists; will not overwrite."
@@ -92,29 +93,31 @@ def process(dirList, inputRegistry, outputRegistry="registry.sqlite3"):
         conn.execute("CREATE INDEX ix_skyTile_tile ON raw_skyTile (skyTile)")
         conn.close()
 
+
 def processVisit(visitDir, conn, done, qsp):
     print >>sys.stderr, visitDir, "... started"
     for raftDir in glob.glob(os.path.join(visitDir, "E00[01]", "R[0-4][0-4]")):
         processRaft(raftDir, conn, done, qsp)
     print >>sys.stderr, visitDir, "... completed"
 
+
 def processRaft(raftDir, conn, done, qsp):
     nProcessed = 0
     nSkipped = 0
     nUnrecognized = 0
     for fits in glob.glob(os.path.join(raftDir, "S[0-2][0-2]",
-        "imsim_*_R[0-4][0-4]_S[0-2][0-2]_C[01][0-7]_E00[01].fits*")):
+                                       "imsim_*_R[0-4][0-4]_S[0-2][0-2]_C[01][0-7]_E00[01].fits*")):
         m = re.search(r'v(\d+)-f(\w)/E00(\d)/R(\d)(\d)/S(\d)(\d)/' +
-                r'imsim_\1_R\4\5_S\6\7_C(\d)(\d)_E00\3\.fits', fits)
+                      r'imsim_\1_R\4\5_S\6\7_C(\d)(\d)_E00\3\.fits', fits)
         if not m:
             print >>sys.stderr, "Warning: Unrecognized file:", fits
             nUnrecognized += 1
             continue
 
         (visit, filter, snap, raft1, raft2, sensor1, sensor2,
-                channel1, channel2) = m.groups()
+         channel1, channel2) = m.groups()
         key = "%s_F%s_E%s_R%s,%s_S%s,%s_C%s,%s" % (visit, filter,
-                snap, raft1, raft2, sensor1, sensor2, channel1, channel2)
+                                                   snap, raft1, raft2, sensor1, sensor2, channel1, channel2)
         if key in done:
             nSkipped += 1
             continue
@@ -126,30 +129,30 @@ def processRaft(raftDir, conn, done, qsp):
                 dafBase.DateTime.TAI).toString(dafBase.DateTime.UTC)[:-1]
         conn.execute("""INSERT INTO raw VALUES
             (NULL, ?, ?, ?, ?, ?, ?, ?, ?)""",
-            (visit, filter, snap, "%s,%s" % (raft1, raft2),
-                "%s,%s" % (sensor1, sensor2),
-                "%s,%s" % (channel1, channel2), taiObs, expTime))
+                     (visit, filter, snap, "%s,%s" % (raft1, raft2),
+                      "%s,%s" % (sensor1, sensor2),
+                      "%s,%s" % (channel1, channel2), taiObs, expTime))
 
         for row in conn.execute("SELECT last_insert_rowid()"):
             id = row[0]
             break
 
         wcs = afwImage.makeWcs(md)
         poly = skypix.imageToPolygon(wcs,
-                md.get("NAXIS1"), md.get("NAXIS2"),
-                padRad=0.000075) # about 15 arcsec
+                                     md.get("NAXIS1"), md.get("NAXIS2"),
+                                     padRad=0.000075)  # about 15 arcsec
         pix = qsp.intersect(poly)
         for skyTileId in pix:
             conn.execute("INSERT INTO raw_skyTile VALUES(?, ?)",
-                    (id, skyTileId))
+                         (id, skyTileId))
 
         conn.commit()
 
         nProcessed += 1
 
     print >>sys.stderr, raftDir, \
-            "... %d processed, %d skipped, %d unrecognized" % \
-            (nProcessed, nSkipped, nUnrecognized)
+        "... %d processed, %d skipped, %d unrecognized" % \
+        (nProcessed, nSkipped, nUnrecognized)
 
 if __name__ == "__main__":
     parser = OptionParser(usage="""%prog [options] DIR ...
@@ -158,7 +161,7 @@ def processRaft(raftDir, conn, done, qsp):
 or a visit subdirectory.""")
     parser.add_option("-i", dest="inputRegistry", help="input registry")
     parser.add_option("-o", dest="outputRegistry", default="registry.sqlite3",
-            help="output registry (default=registry.sqlite3)")
+                      help="output registry (default=registry.sqlite3)")
     (options, args) = parser.parse_args()
     if len(args) < 1:
         parser.error("Missing directory argument(s)")

bin.src/ingestSimImages.py

@@ -1,3 +1,3 @@
 #!/usr/bin/env python
 from lsst.obs.lsstSim.ingest import SimIngestTask
-SimIngestTask.parseAndRun()
+SimIngestTask.parseAndRun()

bin.src/makeDefectMaps.py

@@ -36,18 +36,18 @@
 if len(sys.argv) > 5:
     ha = str(sys.argv[5])
 if ha is None:
-    mi = ai.MaskedImageF("QE_R%i%i_S%i%i.fits.gz"%(rx,ry,sx,sy))
+    mi = ai.MaskedImageF("QE_R%i%i_S%i%i.fits.gz"%(rx, ry, sx, sy))
 elif ha == 'A':
-    mi = ai.MaskedImageF("QE_R%i%i_S%i%i_C0.fits.gz"%(rx,ry,sx,sy))
+    mi = ai.MaskedImageF("QE_R%i%i_S%i%i_C0.fits.gz"%(rx, ry, sx, sy))
 elif ha == 'B':
-    mi = ai.MaskedImageF("QE_R%i%i_S%i%i_C1.fits.gz"%(rx,ry,sx,sy))
+    mi = ai.MaskedImageF("QE_R%i%i_S%i%i_C1.fits.gz"%(rx, ry, sx, sy))
 else:
     raise ValueError("passed an invalid value for ha")
 
 im = mi.getImage()
 arr = im.getArray()
 deadidx = numpy.where(arr == 0.)
-hotidx = numpy.where(arr> 1.1)
+hotidx = numpy.where(arr > 1.1)
 arr[deadidx] = 2.
 arr[hotidx] = 3.
 im2 = ai.makeImageFromArray(arr)
@@ -66,16 +66,17 @@
         height.append(bbox.getHeight())
 
 head = pyfits.Header()
-cmap = {'A':(0,0), 'B':(0,1)}
+cmap = {'A': (0, 0), 'B': (0, 1)}
 if ha is not None:
-    head.update('SERIAL',int('%i%i%i%i%i%i'%(rx,ry,sx,sy,cmap[ha][0],cmap[ha][1])),'Serial of the sensor')
-    head.update('NAME','R:%i,%i S:%i,%i,%c'%(rx,ry,sx,sy,ha),'Name of sensor for this defect map')
+    head.update('SERIAL', int('%i%i%i%i%i%i' %
+                              (rx, ry, sx, sy, cmap[ha][0], cmap[ha][1])), 'Serial of the sensor')
+    head.update('NAME', 'R:%i,%i S:%i,%i,%c'%(rx, ry, sx, sy, ha), 'Name of sensor for this defect map')
 else:
-    head.update('SERIAL',int('%i%i%i%i'%(rx,ry,sx,sy)),'Serial of the sensor')
-    head.update('NAME','R:%i,%i S:%i,%i'%(rx,ry,sx,sy),'Name of sensor for this defect map')
-head.update('CDATE',time.asctime(time.gmtime()),'UTC of creation')
+    head.update('SERIAL', int('%i%i%i%i'%(rx, ry, sx, sy)), 'Serial of the sensor')
+    head.update('NAME', 'R:%i,%i S:%i,%i'%(rx, ry, sx, sy), 'Name of sensor for this defect map')
+head.update('CDATE', time.asctime(time.gmtime()), 'UTC of creation')
 
-#Need to transpose from the phosim on disk orientation
+# Need to transpose from the phosim on disk orientation
 col1 = pyfits.Column(name='y0', format='I', array=numpy.array(y0))
 col2 = pyfits.Column(name='x0', format='I', array=numpy.array(x0))
 col3 = pyfits.Column(name='width', format='I', array=numpy.array(width))
@@ -85,6 +86,6 @@
 hdu = pyfits.PrimaryHDU()
 thdulist = pyfits.HDUList([hdu, tbhdu])
 if ha is None:
-    thdulist.writeto("defects%i%i%i%i.fits"%(rx,ry,sx,sy))
+    thdulist.writeto("defects%i%i%i%i.fits"%(rx, ry, sx, sy))
 else:
-    thdulist.writeto("defects%i%i%i%i%c.fits"%(rx,ry,sx,sy,ha))
+    thdulist.writeto("defects%i%i%i%i%c.fits"%(rx, ry, sx, sy, ha))

bin.src/processCalibLsstSim.py

@@ -31,8 +31,8 @@
     parser = ArgumentParser(name="argumentParser")
     namespace = parser.parse_args(config=TaskClass.ConfigClass())
     sensorDataRefLists = {}
-    typeMap = {'0':'bias', '1':'dark', '2':'flat_u', '3':'flat_g', '4':'flat_r',
-               '5':'flat_i', '6':'flat_z', '7':'flat_y'}
+    typeMap = {'0': 'bias', '1': 'dark', '2': 'flat_u', '3': 'flat_g', '4': 'flat_r',
+               '5': 'flat_i', '6': 'flat_z', '7': 'flat_y'}
     types = []
     for dr in namespace.dataRefList:
         tdict = eval(dr.dataId.__repr__())
@@ -53,7 +53,7 @@
     task = TaskClass()
     if type in ('flat_u', 'flat_g', 'flat_r', 'flat_i', 'flat_z', 'flat_y'):
         type = 'flat'
-    
+
     for k in sensorDataRefLists.keys():
         try:
             task.run(sensorDataRefLists[k], type)

python/lsst/obs/lsstSim/__init__.py

@@ -24,3 +24,4 @@
 from .lsstSimIsrTask import *
 from .utils import *
 from .sensorLib import *
+from .makeLsstSimRawVisitInfo import *

python/lsst/obs/lsstSim/lsstSimIsrTask.py

@@ -29,6 +29,7 @@
 
 __all__ = ["LsstSimIsrTask"]
 
+
 class LsstSimIsrConfig(IsrTask.ConfigClass):
     doWriteSnaps = pexConfig.Field(
         dtype = bool,
@@ -47,13 +48,14 @@ class LsstSimIsrConfig(IsrTask.ConfigClass):
 
     def setDefaults(self):
         IsrTask.ConfigClass.setDefaults(self)
-        self.doDark = False # LSSTSims do not include darks at this time
+        self.doDark = False  # LSSTSims do not include darks at this time
         self.snapCombine.averageKeys = ("TAI", "MJD-OBS", "AIRMASS", "AZIMUTH", "ZENITH",
-            "ROTANG", "SPIDANG", "ROTRATE")
+                                        "ROTANG", "SPIDANG", "ROTRATE")
         self.snapCombine.sumKeys = ("EXPTIME", "CREXPTM", "DARKTIME")
 
 
 class LsstSimIsrTask(IsrTask):
+
     """
     \section obs_lsstSim_isr_Debug Debug variables
 
@@ -88,9 +90,9 @@ def unmaskSatHotPixels(self, exposure):
         maskarr = mask.getArray()
         badBitmask = numpy.array(mask.getPlaneBitMask("BAD"), dtype=maskarr.dtype)
         satBitmask = numpy.array(mask.getPlaneBitMask("SAT"), dtype=maskarr.dtype)
-        orBitmask = badBitmask|satBitmask
+        orBitmask = badBitmask | satBitmask
         andMask = ~satBitmask
-        idx = numpy.where((maskarr&orBitmask)==orBitmask)
+        idx = numpy.where((maskarr & orBitmask) == orBitmask)
         maskarr[idx] &= andMask
 
     def saturationInterpolation(self, ccdExposure):
@@ -108,10 +110,10 @@ def saturationInterpolation(self, ccdExposure):
     @pipeBase.timeMethod
     def runDataRef(self, sensorRef):
         """Do instrument signature removal on an exposure
-        
+
         Correct for saturation, bias, overscan, dark, flat..., perform CCD assembly,
         optionally combine snaps, and interpolate over defects and saturated pixels.
-        
+
         If config.doSnapCombine true then combine the two ISR-corrected snaps to produce the final exposure.
         If config.doSnapCombine false then uses ISR-corrected snap 0 as the final exposure.
         In either case, the final exposure is persisted as "postISRCCD" if config.doWriteSpans is True,
@@ -133,14 +135,14 @@ def runDataRef(self, sensorRef):
             isrData = self.readIsrData(snapRef, ccdExposure)
             ccdExposure = self.run(ccdExposure, **isrData.getDict()).exposure
             snapDict[snapId] = ccdExposure
-    
+
             if self.config.doWriteSnaps:
                 sensorRef.put(ccdExposure, "snapExp", snap=snapId)
 
             frame = getDebugFrame(self._display, "snapExp%d" % (snapId,))
             if frame:
                 getDisplay(frame).mtv(ccdExposure)
-        
+
         if self.config.doSnapCombine:
             loadSnapDict(snapDict, snapIdList=(0, 1), sensorRef=sensorRef)
             postIsrExposure = self.snapCombine.run(snapDict[0], snapDict[1]).exposure
@@ -155,14 +157,15 @@ def runDataRef(self, sensorRef):
         frame = getDebugFrame(self._display, "postISRCCD")
         if frame:
             getDisplay(frame).mtv(postIsrExposure)
-                
+
         return pipeBase.Struct(
             exposure = postIsrExposure,
         )
 
+
 def loadSnapDict(snapDict, snapIdList, sensorRef):
     """Load missing snaps from disk.
-    
+
     @paramp[in,out] snapDict: a dictionary of snapId: snap exposure ("snapExp")
     @param[in] snapIdList: a list of snap IDs
     @param[in] sensorRef: sensor reference for snap, excluding the snap ID.
@@ -173,4 +176,3 @@ def loadSnapDict(snapDict, snapIdList, sensorRef):
             if snapExposure is None:
                 raise RuntimeError("Could not find snapExp for snap=%s; id=%s" % (snapId, sensorRef.dataId))
             snapDict[snapId] = snapExposure
-

python/lsst/obs/lsstSim/lsstSimMapper.py

@@ -20,6 +20,7 @@
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
 
+import math
 import re
 
 import lsst.daf.base as dafBase
@@ -28,6 +29,7 @@
 import lsst.afw.coord as afwCoord
 import lsst.afw.geom as afwGeom
 import lsst.pex.policy as pexPolicy
+from .makeLsstSimRawVisitInfo import MakeLsstSimRawVisitInfo
 
 from lsst.daf.butlerUtils import CameraMapper
 
@@ -37,6 +39,8 @@
 class LsstSimMapper(CameraMapper):
     packageName = 'obs_lsstSim'
 
+    MakeRawVisitInfoClass = MakeLsstSimRawVisitInfo
+
     _CcdNameRe = re.compile(r"R:(\d,\d) S:(\d,\d(?:,[AB])?)$")
 
     def __init__(self, inputPolicy=None, **kwargs):
@@ -219,8 +223,9 @@ def _computeCoaddExposureId(self, dataId, singleFilter):
         if singleFilter:
             return id * 8 + self.filterIdMap[dataId['filter']]
         return id
-    
+
     _nbit_id = 30
+
     def bypass_deepMergedCoaddId_bits(self, *args, **kwargs):
         """The number of bits used up for patch ID bits"""
         return 64 - self._nbit_id
@@ -248,12 +253,12 @@ def _setCcdExposureId(self, propertyList, dataId):
 
     def std_raw(self, item, dataId):
         exposure = super(LsstSimMapper, self).std_raw(item, dataId)
-
         md = exposure.getMetadata()
         if md.exists("VERSION") and md.getInt("VERSION") < 16952:
-            # Precess WCS based on actual observation date
-            epoch = dafBase.DateTime(md.get("MJD-OBS"), dafBase.DateTime.MJD,
-                                     dafBase.DateTime.TAI).get(dafBase.DateTime.EPOCH)
+            # Precess crval of WCS from date of observation to J2000
+            epoch = exposure.getInfo().getVisitInfo().getDate().get(dafBase.DateTime.EPOCH)
+            if math.isnan(epoch):
+                raise RuntimeError("Date not found in raw exposure %s" % (dataId,))
             wcs = exposure.getWcs()
             origin = wcs.getSkyOrigin()
             refCoord = afwCoord.Fk5Coord(