DM-16490: Update computeCcdOffsets() to use makeSkyWcs

cookbook/fgcmFitCycleHscCookbook_cycle00_config.py

@@ -29,8 +29,6 @@
 config.epochMjds = (56700., 57420., 57606.)
 # Latitude of the observatory
 config.latitude = 19.8256
-# Pixel scale (arcseconds)
-config.pixelScale = 0.17
 # Amount of gray extinction to be considered "photometric".  This will
 # get ratcheded down in further cycles.
 config.expGrayPhotometricCut = (-0.05, -0.05, -0.05, -0.05, -0.05)

python/lsst/fgcmcal/fgcmBuildStars.py

@@ -431,6 +431,7 @@ def _fillVisitCatalog(self, visitCat, groupedDataRefs):
             rec['telra'] = radec.getRa().asDegrees()
             rec['teldec'] = radec.getDec().asDegrees()
             rec['telha'] = visitInfo.getBoresightHourAngle().asDegrees()
+            rec['telrot'] = visitInfo.getBoresightRotAngle().asDegrees()
             rec['mjd'] = visitInfo.getDate().get(system=DateTime.MJD)
             rec['exptime'] = visitInfo.getExposureTime()
             # convert from Pa to millibar
@@ -849,6 +850,7 @@ def _makeFgcmVisitSchema(self, nCcd):
         schema.addField('telra', type=np.float64, doc="Pointing RA (deg)")
         schema.addField('teldec', type=np.float64, doc="Pointing Dec (deg)")
         schema.addField('telha', type=np.float64, doc="Pointing Hour Angle (deg)")
+        schema.addField('telrot', type=np.float64, doc="Camera rotation (deg)")
         schema.addField('mjd', type=np.float64, doc="MJD of visit")
         schema.addField('exptime', type=np.float32, doc="Exposure time")
         schema.addField('pmb', type=np.float32, doc="Pressure (millibar)")

python/lsst/fgcmcal/fgcmCalibrateTract.py

@@ -260,7 +260,10 @@ def runDataRef(self, butler, dataRefs):
                                     True, False, tract=tract)
         # Turn off plotting in tract mode
         configDict['doPlots'] = False
-        ccdOffsets = computeCcdOffsets(camera, self.config.fgcmFitCycle.pixelScale)
+
+        # Use the first orientation.
+        # TODO: DM-21215 will generalize to arbitrary camera orientations
+        ccdOffsets = computeCcdOffsets(camera, fgcmExpInfo['TELROT'][0])
         del camera
 
         # Set up the fit cycle task

python/lsst/fgcmcal/fgcmFitCycle.py

@@ -214,11 +214,12 @@ class FgcmFitCycleConfig(pexConfig.Config):
         dtype=float,
         default=None,
     )
-    # TODO: DM-16490 will make this unneccessary
     pixelScale = pexConfig.Field(
         doc="Pixel scale (arcsec/pixel) (temporary)",
         dtype=float,
-        default=None,
+        deprecated=("This field is no longer used, and has been deprecated by DM-16490. "
+                    "It will be removed after v19."),
+        optional=True,
     )
     brightObsGrayMax = pexConfig.Field(
         doc="Maximum gray extinction to be considered bright observation",
@@ -612,7 +613,9 @@ def _fgcmFitCycle(self, butler):
         fgcmExpInfo = translateVisitCatalog(visitCat)
         del visitCat
 
-        ccdOffsets = computeCcdOffsets(camera, self.config.pixelScale)
+        # Use the first orientation.
+        # TODO: DM-21215 will generalize to arbitrary camera orientations
+        ccdOffsets = computeCcdOffsets(camera, fgcmExpInfo['TELROT'][0])
 
         noFitsDict = {'lutIndex': lutIndexVals,
                       'lutStd': lutStd,

python/lsst/fgcmcal/utilities.py

@@ -30,7 +30,9 @@
 
 import lsst.afw.cameraGeom as afwCameraGeom
 import lsst.afw.table as afwTable
+import lsst.afw.image as afwImage
 import lsst.geom as geom
+from lsst.obs.base import createInitialSkyWcs
 
 import fgcm
 
@@ -348,6 +350,7 @@ def translateVisitCatalog(visitCat):
                                                  ('TELHA', 'f8'),
                                                  ('TELRA', 'f8'),
                                                  ('TELDEC', 'f8'),
+                                                 ('TELROT', 'f8'),
                                                  ('PMB', 'f8'),
                                                  ('FILTERNAME', 'a2')])
     fgcmExpInfo['VISIT'][:] = visitCat['visit']
@@ -357,6 +360,7 @@ def translateVisitCatalog(visitCat):
     fgcmExpInfo['TELHA'][:] = visitCat['telha']
     fgcmExpInfo['TELRA'][:] = visitCat['telra']
     fgcmExpInfo['TELDEC'][:] = visitCat['teldec']
+    fgcmExpInfo['TELROT'][:] = visitCat['telrot']
     fgcmExpInfo['PMB'][:] = visitCat['pmb']
     fgcmExpInfo['PSFSIGMA'][:] = visitCat['psfSigma']
     fgcmExpInfo['DELTA_APER'][:] = visitCat['deltaAper']
@@ -368,22 +372,23 @@ def translateVisitCatalog(visitCat):
     return fgcmExpInfo
 
 
-def computeCcdOffsets(camera, pixelScale):
+def computeCcdOffsets(camera, defaultOrientation):
     """
     Compute the CCD offsets in ra/dec and x/y space
 
     Parameters
     ----------
-    camera: ``
-    pixelScale: `float`
-       Pixel scale of the camera.  Required because of HSC.
+    camera: `lsst.afw.cameraGeom.Camera`
+    defaultOrientation: `float`
+       Default camera orientation (degrees)
 
     Returns
     -------
     ccdOffsets: `numpy.ndarray`
-       Numpy array with ccd offset information for input to FGCM
+       Numpy array with ccd offset information for input to FGCM.
+       Angular units are degrees, and x/y units are pixels.
     """
-    # TODO: DM-16490 will simplify and generalize the math.
+    # TODO: DM-21215 will fully generalize to arbitrary camera orientations
 
     # and we need to know the ccd offsets from the camera geometry
     ccdOffsets = np.zeros(len(camera), dtype=[('CCDNUM', 'i4'),
@@ -394,30 +399,40 @@ def computeCcdOffsets(camera, pixelScale):
                                               ('X_SIZE', 'i4'),
                                               ('Y_SIZE', 'i4')])
 
+    # Generate fake WCSs centered at 180/0 to avoid the RA=0/360 problem,
+    # since we are looking for relative positions
+    boresight = geom.SpherePoint(180.0*geom.degrees, 0.0*geom.degrees)
+
+    # TODO: DM-17597 will update testdata_jointcal so that the test data
+    # does not have nan as the boresight angle for HSC data.  For the
+    # time being, there is this ungainly hack.
+    if camera.getName() == 'HSC' and np.isnan(defaultOrientation):
+        orientation = 270*geom.degrees
+    else:
+        orientation = defaultOrientation*geom.degrees
+    flipX = False
+
+    # Create a temporary visitInfo for input to createInitialSkyWcs
+    visitInfo = afwImage.VisitInfo(boresightRaDec=boresight,
+                                   boresightRotAngle=orientation,
+                                   rotType=afwImage.visitInfo.RotType.SKY)
+
     for i, detector in enumerate(camera):
-        # new version, using proper rotations
-        #  but I worry this only works with HSC, as there's a unit inconsistency
-        camPoint = detector.getCenter(afwCameraGeom.FOCAL_PLANE)
-        pixelSize = detector.getPixelSize()[0]  # Assumes x and y pixel sizes in arcsec are the same
-        camPoint.scale(1.0/pixelSize)
+        ccdOffsets['CCDNUM'][i] = detector.getId()
+
+        wcs = createInitialSkyWcs(visitInfo, detector, flipX)
+
+        detCenter = wcs.pixelToSky(detector.getCenter(afwCameraGeom.PIXELS))
+        ccdOffsets['DELTA_RA'][i] = (detCenter.getRa() - boresight.getRa()).asDegrees()
+        ccdOffsets['DELTA_DEC'][i] = (detCenter.getDec() - boresight.getDec()).asDegrees()
+
         bbox = detector.getBBox()
 
-        ccdOffsets['CCDNUM'][i] = detector.getId()
-        # this requires a pixelScale
-        # Note that this now works properly with HSC, but I need to work on
-        # generalizing this properly.  I expect the updates in DM-16490 will
-        # generalize these computations, and appropriate tests can be added
-        # on that ticket.
-        ccdOffsets['DELTA_RA'][i] = -camPoint.getY() * pixelScale / 3600.0
-        ccdOffsets['DELTA_DEC'][i] = -camPoint.getX() * pixelScale / 3600.0
-        boxXform = detector.transform(geom.Point2D(bbox.getMaxX(), bbox.getMaxY()),
-                                      afwCameraGeom.PIXELS, afwCameraGeom.FOCAL_PLANE)
-        boxXform.scale(1.0/pixelSize)
-        # but this does not (for the delta)
-        ccdOffsets['RA_SIZE'][i] = 2. * np.abs(boxXform.getY() -
-                                               camPoint.getY()) / 3600.0
-        ccdOffsets['DEC_SIZE'][i] = 2. * np.abs(boxXform.getX() -
-                                                camPoint.getX()) / 3600.0
+        detCorner1 = wcs.pixelToSky(geom.Point2D(bbox.getMin()))
+        detCorner2 = wcs.pixelToSky(geom.Point2D(bbox.getMax()))
+
+        ccdOffsets['RA_SIZE'][i] = np.abs((detCorner2.getRa() - detCorner1.getRa()).asDegrees())
+        ccdOffsets['DEC_SIZE'][i] = np.abs((detCorner2.getDec() - detCorner1.getDec()).asDegrees())
 
         ccdOffsets['X_SIZE'][i] = bbox.getMaxX()
         ccdOffsets['Y_SIZE'][i] = bbox.getMaxY()

tests/test_ccdoffsets.py

@@ -0,0 +1,111 @@
+# See COPYRIGHT file at the top of the source tree.
+#
+# This file is part of fgcmcal.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# 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 GNU General Public License
+# along with this program.  If not, see <https://www.gnu.org/licenses/>.
+"""Test the fgcmcal computeCcdOffsets code with testdata_jointcal.
+"""
+
+import unittest
+import os
+
+import lsst.utils
+import lsst.daf.persistence as dafPersist
+
+import lsst.fgcmcal as fgcmcal
+
+
+class FgcmCcdOffsetsTest(lsst.utils.tests.TestCase):
+    @classmethod
+    def setUpClass(cls):
+        try:
+            cls.dataDir = lsst.utils.getPackageDir('testdata_jointcal')
+        except lsst.pex.exceptions.NotFoundError:
+            raise unittest.SkipTest("testdata_jointcal not setup")
+
+    def test_fgcmCcdOffsetsHsc(self):
+        """
+        Test computation of ccd offsets for HSC.
+        """
+
+        lsst.log.setLevel("HscMapper", lsst.log.FATAL)
+
+        butler = dafPersist.Butler(os.path.join(self.dataDir, 'hsc'))
+
+        visit = 903986
+        ccd = 16
+
+        visitInfo = butler.get('calexp_visitInfo', visit=visit, ccd=ccd)
+        rotAngle = visitInfo.getBoresightRotAngle().asDegrees()
+
+        camera = butler.get('camera')
+
+        ccdOffsets = fgcmcal.utilities.computeCcdOffsets(camera, rotAngle)
+
+        # Spot check relative orientations of some of the ccds
+        # This is based on
+        # https://subarutelescope.org/Observing/Instruments/HSC/CCDPosition_20170212.png
+        # The goal here is to check that North is Up, South is Down,
+        # East is Left, and West is Right.
+        self.assertLess(ccdOffsets['DELTA_RA'][15], ccdOffsets['DELTA_RA'][10])
+        self.assertLess(ccdOffsets['DELTA_RA'][95], ccdOffsets['DELTA_RA'][90])
+        self.assertGreater(ccdOffsets['DELTA_RA'][46], 0.0)
+
+        self.assertLess(ccdOffsets['DELTA_DEC'][15], ccdOffsets['DELTA_DEC'][95])
+        self.assertLess(ccdOffsets['DELTA_DEC'][10], ccdOffsets['DELTA_DEC'][90])
+        self.assertGreater(ccdOffsets['DELTA_DEC'][97], 0.0)
+
+        # Check the sizes
+        # The projected size of the ccds varies with radius over the large
+        # HSC field-of-view. Empirically, the x size is between 0.07 and 0.10 deg
+        # and the y size is between 0.17 and 0.20 deg for the non-rotated CCDs.
+        # This test checks that the orientations of the CCDs are as expected for
+        # rotated/non-rotated CCDs (relative size of RA/DEC), and that the size
+        # is roughly correct.  Because these values are only used for visualization
+        # in the fgcm code, this does not need to be perfect.
+
+        rotatedCcds = [100, 101, 102, 103]
+
+        for i, ccd in enumerate(ccdOffsets['CCDNUM']):
+            if ccd in rotatedCcds:
+                self.assertLess(ccdOffsets['RA_SIZE'][i], ccdOffsets['DEC_SIZE'][i])
+                self.assertGreater(ccdOffsets['DEC_SIZE'][i], 0.17)
+                self.assertLess(ccdOffsets['DEC_SIZE'][i], 0.20)
+                self.assertGreater(ccdOffsets['RA_SIZE'][i], 0.07)
+                self.assertLess(ccdOffsets['RA_SIZE'][i], 0.10)
+            else:
+                self.assertGreater(ccdOffsets['RA_SIZE'][i], ccdOffsets['DEC_SIZE'][i])
+                self.assertGreater(ccdOffsets['RA_SIZE'][i], 0.17)
+                self.assertLess(ccdOffsets['RA_SIZE'][i], 0.20)
+                self.assertGreater(ccdOffsets['DEC_SIZE'][i], 0.07)
+                self.assertLess(ccdOffsets['DEC_SIZE'][i], 0.10)
+
+
+class TestMemory(lsst.utils.tests.MemoryTestCase):
+    pass
+
+
+def setup_module(module):
+    lsst.utils.tests.init()
+
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()

tests/test_fgcmcal_hsc.py

@@ -267,7 +267,6 @@ def fillDefaultFitCycleConfig(self, config):
         config.washMjds = (0.0, )
         config.epochMjds = (0.0, 100000.0)
         config.latitude = 19.8256
-        config.pixelScale = 0.17
         config.expGrayPhotometricCut = (-0.05, -0.05)
         config.expGrayHighCut = (0.2, 0.2)
         config.aperCorrFitNBins = 0