Merge pull request #18 from lsst/tickets/DM-6349

Passed some tests with twinkles data. Merging before I rewrite the Readme to use the PhoSim Globus data. Current version is in sync with PhoSim 3.4.2.

Readme.md

@@ -0,0 +1,28 @@
+Description
+===========
+
+This obs_lsstSim package provides an interface to the phosim output for the LSST Data Management software.
+
+Updating camera description
+---------------------------
+
+The camera description FITS files are built from the phosim lsst data text files. To update the camera description to match a new version of phosim, you will need two files from `phosim_release/data/lsst`, and the gain and saturation data from the header of the phosim generated files for the entire focalplane.
+
+1. Checkout the latest version of phosim:
+    * `git clone https://bitbucket.org/phosim/phosim_release.git`
+2. Copy these files from your phosim checkout into this repo's `description/` directory:
+    * `description/lsst/focalplanelayout.txt`
+    * `description/lsst/segmentation.txt`
+3. Generate the gain/saturation file (with obs_lsstSim setup via eups):
+    1. Build phosim (see [Using Phosim](https://bitbucket.org/phosim/phosim_release/wiki/Using%20PhoSim) for build instructions). You likely will need to create an SEDs directory (even if its empty):
+    `mkdir data/SEDs`
+    2. Run phosim to generate a simple no-background simulation (give this several hours to complete):
+    `./phosim $OBS_LSSTSIM_DIR/description/nostars_allchips -c examples/nobackground`
+    3. Extract the gain and saturation values from the headers produced by the above commands:
+    `extractPhosimGainSaturation.py -v output/`
+3. Update the phosim verison in this repo's `description/phosim_version.txt` to match the version associated with the above files.
+4. Commit your above changes and push.
+
+The files you have updated above will be converted into the camera description files when this product is built by scons. You can check that this completes successfully by setting up this product and running "scons" at the root level.
+
+Note that this process does not build the wavefront sensor files.

SConstruct

@@ -1,3 +1,13 @@
 # -*- python -*-
+import lsst.sconsUtils
 from lsst.sconsUtils import scripts
-scripts.BasicSConstruct("obs_lsstSim")
+scripts.BasicSConstruct("obs_lsstSim",
+                        defaultTargets=scripts.DEFAULT_TARGETS + ("description",),
+                        sconscriptOrder=["python", "description", "tests"]
+                        )
+env = lsst.sconsUtils.env
+# Don't need to rebuild the camera descriptions if this code changes, so just Requires(), not Depends().
+env.Requires(lsst.sconsUtils.targets['description'], lsst.sconsUtils.targets['version'])
+env.Requires(lsst.sconsUtils.targets['description'], lsst.sconsUtils.targets['python'])
+# Do need to rerun tests anytime the camera descriptions change.
+env.Depends(lsst.sconsUtils.targets['tests'], lsst.sconsUtils.targets['description'])

bin.src/extractPhosimGainSaturation.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python2
+#
+# LSST Data Management System
+# Copyright 2014 LSST Corporation.
+#
+# 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/>.
+#
+"""
+Extract gain and saturation values from phosim headers to make a gainFile for
+makeLsstCameraRepository.py.
+
+You first have to run phosim to generate images for each amp. These can be the
+simplest no-background, no-source images, so long as they have a real gain and
+saturation threshold.
+"""
+from __future__ import absolute_import, division, print_function
+
+import os
+import glob
+
+from astropy.io import fits
+
+
+def read_files(path, verbose=False):
+    """Return a dictionary of amp name: (gain, saturation) for all amp files in path."""
+    amps = {}
+    files = glob.glob(os.path.join(path, 'lsst_*_R??_S??_C??*.fits.gz'))
+    for file in files:
+        if verbose:
+            print(file)
+        header = fits.getheader(file)
+        ampid = '_'.join((header['CCDID'], header['AMPID']))
+        amps[ampid] = header['GAIN'], int(header['SATURATE'])
+    return amps
+
+
+def write_gain_file(amps, filename):
+    """Write a gainFile for use by makeLsstCameraRepository.py."""
+    with open(filename, 'w') as outfile:
+        for amp in sorted(amps):
+            line = '{} {} {}\n'.format(amp, amps[amp][0], amps[amp][1])
+            outfile.write(line)
+
+
+def main():
+    import argparse
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("phosim_output_path",
+                        help="Path to phosim output directory containing amp files.")
+    parser.add_argument("--verbose", "-v", action="store_true", help="Show file processing prgress.")
+    args = parser.parse_args()
+
+    amps = read_files(args.phosim_output_path, verbose=args.verbose)
+    filename = os.path.join(os.environ['OBS_LSSTSIM_DIR'], 'description', 'gain_saturation.txt')
+    write_gain_file(amps, filename)
+    print("Wrote: ", filename)
+
+
+if __name__ == "__main__":
+    main()

bin.src/makeLsstCameraRepository.py

@@ -1,7 +1,7 @@
 #!/usr/bin/env python2
 #
 # LSST Data Management System
-# Copyright 2014 LSST Corporation.
+# Copyright 2014-2016 LSST Corporation.
 #
 # This product includes software developed by the
 # LSST Project (http://www.lsst.org/).
@@ -20,6 +20,12 @@
 # the GNU General Public License along with this program.  If not,
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
+"""
+Produce the camera description FITS files from phosim text files.
+
+Scons should have automatically run this when building obs_lsstSim. To produce
+the same files that scons would have, run with no arguments.
+"""
 from __future__ import absolute_import, division
 import argparse
 import os
@@ -32,6 +38,7 @@
 from lsst.afw.cameraGeom import DetectorConfig, CameraConfig, PUPIL, FOCAL_PLANE, PIXELS, NullLinearityType
 from lsst.obs.lsstSim import LsstSimMapper
 
+
 def expandDetectorName(abbrevName):
     """Convert a detector name of the form Rxy_Sxy[_Ci] to canonical form: R:x,y S:x,y[,c]
 
@@ -43,9 +50,10 @@ def expandDetectorName(abbrevName):
     fullName = "R:%s,%s S:%s,%s" % tuple(m.groups()[0:4])
     subSensor = m.groups()[4]
     if subSensor is not None:
-        fullName  = fullName + "," + {"0": "A", "1": "B"}[subSensor]
+        fullName = fullName + "," + {"0": "A", "1": "B"}[subSensor]
     return fullName
 
+
 def detectorIdFromAbbrevName(abbrevName):
     """Compute detector ID from an abbreviated detector name of the form Rxy_Sxy_Ci
 
@@ -59,24 +67,29 @@ def detectorIdFromAbbrevName(abbrevName):
         detectorId += 10000 * (1 + int(m.group(5)))
     return detectorId
 
+
 def makeAmpTables(segmentsFile, gainFile):
     """
     Read the segments file from a PhoSim release and produce the appropriate AmpInfo
-    @param segmentsFile -- String indicating where the file is located
+
+    @param segmentsFile (str) full path to the segmentation file.
+    @param gainFile (str) full path to the gain/saturation file.
+
+    @return (dict) per amp dictionary of ampCatalogs
     """
     gainDict = {}
     with open(gainFile) as fh:
         for l in fh:
             els = l.rstrip().split()
-            gainDict[els[0]] = {'gain':float(els[1]), 'saturation':int(els[2])}
+            gainDict[els[0]] = {'gain': float(els[1]), 'saturation': int(els[2])}
     returnDict = {}
-    #TODO currently there is no linearity provided, but we should identify
-    #how to get this information.
-    linearityCoeffs = (0.,1.,0.,0.)
+    # TODO currently there is no linearity provided, but we should identify
+    # how to get this information.
+    linearityCoeffs = (0., 1., 0., 0.)
     linearityType = NullLinearityType
-    readoutMap = {'LL':afwTable.LL, 'LR':afwTable.LR, 'UR':afwTable.UR, 'UL':afwTable.UL}
+    readoutMap = {'LL': afwTable.LL, 'LR': afwTable.LR, 'UR': afwTable.UR, 'UL': afwTable.UL}
     ampCatalog = None
-    detectorName = [] # set to a value that is an invalid dict key, to catch bugs
+    detectorName = []  # set to a value that is an invalid dict key, to catch bugs
     correctY0 = False
     with open(segmentsFile) as fh:
         for l in fh:
@@ -91,22 +104,22 @@ def makeAmpTables(segmentsFile, gainFile):
                 numy = int(els[2])
                 schema = afwTable.AmpInfoTable.makeMinimalSchema()
                 ampCatalog = afwTable.AmpInfoCatalog(schema)
-                if len(els[0].split('_')) == 3:   #wavefront sensor
+                if len(els[0].split('_')) == 3:  # wavefront sensor
                     correctY0 = True
                 else:
                     correctY0 = False
                 continue
             record = ampCatalog.addNew()
             name = els[0].split("_")[-1]
             name = '%s,%s'%(name[1], name[2])
-            #Because of the camera coordinate system, we choose an
-            #image coordinate system that requires a -90 rotation to get
-            #the correct pixel positions from the
-            #phosim segments file
+            # Because of the camera coordinate system, we choose an
+            # image coordinate system that requires a -90 rotation to get
+            # the correct pixel positions from the
+            # phosim segments file
             y0 = numy - 1 - int(els[2])
             y1 = numy - 1 - int(els[1])
-            #Another quirk of the phosim file is that one of the wavefront sensor
-            #chips has an offset of 2000 pix in y.  It's always the 'C1' chip.
+            # Another quirk of the phosim file is that one of the wavefront sensor
+            # chips has an offset of 2000 pix in y.  It's always the 'C1' chip.
             if correctY0:
                 if y0 > 0:
                     y1 -= y0
@@ -132,8 +145,8 @@ def makeAmpTables(segmentsFile, gainFile):
             else:
                 flipy = True
 
-            #Since the amps are stored in amp coordinates, the readout is the same
-            #for all amps
+            # Since the amps are stored in amp coordinates, the readout is the same
+            # for all amps
             readCorner = readoutMap['LL']
 
             ndatax = x1 - x0 + 1
@@ -144,7 +157,7 @@ def makeAmpTables(segmentsFile, gainFile):
             hoverscan = 0
             extended = 4
             voverscan = 0
-            rawBBox = afwGeom.Box2I(afwGeom.Point2I(0,0),
+            rawBBox = afwGeom.Box2I(afwGeom.Point2I(0, 0),
                                     afwGeom.Extent2I(extended+ndatax+hoverscan, prescan+ndatay+voverscan))
             rawDataBBox = afwGeom.Box2I(afwGeom.Point2I(extended, prescan), afwGeom.Extent2I(ndatax, ndatay))
             rawHorizontalOverscanBBox = afwGeom.Box2I(afwGeom.Point2I(0, prescan),
@@ -157,7 +170,7 @@ def makeAmpTables(segmentsFile, gainFile):
             extraRawY = prescan + voverscan
             rawx0 = x0 + extraRawX*(x0//ndatax)
             rawy0 = y0 + extraRawY*(y0//ndatay)
-            #Set the elements of the record for this amp
+            # Set the elements of the record for this amp
             record.setBBox(bbox)
             record.setName(name)
             record.setReadoutCorner(readCorner)
@@ -179,6 +192,7 @@ def makeAmpTables(segmentsFile, gainFile):
     returnDict[detectorName] = ampCatalog
     return returnDict
 
+
 def makeLongName(shortName):
     """
     Make the long name from the PhoSim short name
@@ -188,13 +202,14 @@ def makeLongName(shortName):
     if len(parts) == 2:
         return " ".join(["%s:%s"%(el[0], ",".join(el[1:])) for el in parts])
     elif len(parts) == 3:
-        #This must be a wavefront sensor
-        wsPartMap = {'S':{'C0':'A', 'C1':'B'},
-                     'R':{'C0':'', 'C1':''}}
+        # This must be a wavefront sensor
+        wsPartMap = {'S': {'C0': 'A', 'C1': 'B'},
+                     'R': {'C0': '', 'C1': ''}}
         return " ".join(["%s:%s"%(el[0], ",".join(el[1:]+wsPartMap[el[0]][parts[-1]])) for el in parts[:-1]])
     else:
         raise ValueError("Could not parse %s: has %i parts"%(shortName, len(parts)))
 
+
 def makeDetectorConfigs(detectorLayoutFile, phosimVersion):
     """
     Create the detector configs to use in building the Camera
@@ -205,8 +220,8 @@ def makeDetectorConfigs(detectorLayoutFile, phosimVersion):
     * deal with the extra orientation angles (not that they really matter)
     """
     detectorConfigs = []
-    detTypeMap = {"Group2":2, "Group1":3, "Group0":0}
-    #We know we need to rotate 3 times and also apply the yaw perturbation
+    detTypeMap = {"Group2": 2, "Group1": 3, "Group0": 0}
+    # We know we need to rotate 3 times and also apply the yaw perturbation
     nQuarter = 1
     with open(detectorLayoutFile) as fh:
         for l in fh:
@@ -244,60 +259,80 @@ def makeDetectorConfigs(detectorLayoutFile, phosimVersion):
             detectorConfigs.append(detConfig)
     return detectorConfigs
 
+
+def getPhosimVersion(defaultDataDir):
+    """Return the phosim version from data/phosim_version.txt"""
+    with open(os.path.join(defaultDataDir, 'phosim_version.txt')) as infile:
+        return infile.read().strip()
+
+
 if __name__ == "__main__":
     """
     Create the configs for building a camera.  This runs on the files distributed with PhoSim.
 
     Currently gain and saturation need to be supplied as well. The file should have three columns:
         on disk amp id (R22_S11_C00), gain, saturation.
-    For example:
-    - DetectorLayoutFile
-        https://dev.lsstcorp.org/cgit/LSST/sims/phosim.git/plain/data/lsst/focalplanelayout.txt?h=dev
-    - SegmentsFile
-        https://dev.lsstcorp.org/cgit/LSST/sims/phosim.git/plain/data/lsst/segmentation.txt?h=dev
     """
     baseDir = lsst.utils.getPackageDir('obs_lsstsim')
+    defaultDataDir = os.path.join(os.path.normpath(baseDir), "description")
     defaultOutDir = os.path.join(os.path.normpath(baseDir), "description", "camera")
 
-    parser = argparse.ArgumentParser()
-    parser.add_argument("DetectorLayoutFile", help="Path to detector layout file")
-    parser.add_argument("SegmentsFile", help="Path to amp segments file")
-    parser.add_argument("GainFile", help="Path to gain and saturation file")
+    parser = argparse.ArgumentParser(description=__doc__,
+                                     formatter_class=argparse.ArgumentDefaultsHelpFormatter)
+    parser.add_argument("DetectorLayoutFile",
+                        default=os.path.join(defaultDataDir, 'focalplanelayout.txt'),
+                        nargs='?',
+                        help="Path to detector layout file")
+    parser.add_argument("SegmentsFile",
+                        default=os.path.join(defaultDataDir, 'segmentation.txt'),
+                        nargs='?',
+                        help="Path to amp segments file")
+    parser.add_argument("GainFile",
+                        default=os.path.join(defaultDataDir, 'gain_saturation.txt'),
+                        nargs='?',
+                        help="Path to gain and saturation file")
     parser.add_argument("phosimVersion",
-                        help="String id of the version of phosim used to construct this camera repository")
+                        default=None,
+                        nargs='?',
+                        help="String id of the version of phosim used to construct this camera repository."
+                        "If None, use the value in data/phosim_version.txt.")
     parser.add_argument("OutputDir",
-        help = "Path to dump configs and AmpInfo Tables; defaults to %r" % (defaultOutDir,),
-        nargs = "?",
-        default = defaultOutDir,
-    )
+                        help = "Path to dump configs and AmpInfo Tables; defaults to %r" % (defaultOutDir,),
+                        nargs = "?",
+                        default = defaultOutDir,
+                        )
     parser.add_argument("--clobber", action="store_true", dest="clobber", default=False,
-        help=("remove and re-create the output directory if it already exists?"))
+                        help=("remove and re-create the output directory if it already exists?"))
     args = parser.parse_args()
     ampTableDict = makeAmpTables(args.SegmentsFile, args.GainFile)
-    detectorConfigList = makeDetectorConfigs(args.DetectorLayoutFile, args.phosimVersion)
+    if args.phosimVersion is None:
+        phosimVersion = getPhosimVersion(defaultDataDir)
+    else:
+        phosimVersion = args.phosimVersion
+    detectorConfigList = makeDetectorConfigs(args.DetectorLayoutFile, phosimVersion)
 
-    #Build the camera config.
+    # Build the camera config.
     camConfig = CameraConfig()
-    camConfig.detectorList = dict([(i,detectorConfigList[i]) for i in xrange(len(detectorConfigList))])
+    camConfig.detectorList = dict([(i, detectorConfigList[i]) for i in xrange(len(detectorConfigList))])
     camConfig.name = 'LSST'
     camConfig.plateScale = 20.0
     pScaleRad = afwGeom.arcsecToRad(camConfig.plateScale)
     pincushion = 0.925
     # Don't have this yet ticket/3155
-    #camConfig.boresiteOffset_x = 0.
-    #camConfig.boresiteOffset_y = 0.
+    # camConfig.boresiteOffset_x = 0.
+    # camConfig.boresiteOffset_y = 0.
     tConfig = afwGeom.TransformConfig()
     tConfig.transform.name = 'inverted'
     radialClass = afwGeom.xyTransformRegistry['radial']
     tConfig.transform.active.transform.retarget(radialClass)
     # According to Dave M. the simulated LSST transform is well approximated (1/3 pix)
     # by a scale and a pincusion.
     tConfig.transform.active.transform.coeffs = [0., 1./pScaleRad, 0., pincushion/pScaleRad]
-    #tConfig.transform.active.boresiteOffset_x = camConfig.boresiteOffset_x
-    #tConfig.transform.active.boresiteOffset_y = camConfig.boresiteOffset_y
+    # tConfig.transform.active.boresiteOffset_x = camConfig.boresiteOffset_x
+    # tConfig.transform.active.boresiteOffset_y = camConfig.boresiteOffset_y
     tmc = afwGeom.TransformMapConfig()
     tmc.nativeSys = FOCAL_PLANE.getSysName()
-    tmc.transforms = {PUPIL.getSysName():tConfig}
+    tmc.transforms = {PUPIL.getSysName(): tConfig}
     camConfig.transformDict = tmc
 
     def makeDir(dirPath, doClobber=False):