DM-10765: Replace existing WCS classes with SkyWcs

bin.src/genCoaddRegistry.py

@@ -28,13 +28,10 @@
 import os
 import re
 import shutil
-try:
-    import sqlite3
-except ImportError:
-    # try external pysqlite package; deprecated
-    import sqlite as sqlite3
+import sqlite3
 import sys
-import lsst.afw.image as afwImage
+from lsst.afw.fits import readMetadata
+from lsst.afw.geom import makeSkyWcs
 import lsst.skypix as skypix
 
 
@@ -102,15 +99,15 @@ def processBand(filterDir, conn, done, qsp):
             nSkipped += 1
             continue
 
-        md = afwImage.readMetadata(fits)
+        md = readMetadata(fits)
         conn.execute("""INSERT INTO raw VALUES
             (NULL, ?, ?, ?, ?)""", (run, filter, camcol, field))
 
         for row in conn.execute("SELECT last_insert_rowid()"):
             id = row[0]
             break
 
-        wcs = afwImage.makeWcs(md)
+        wcs = makeSkyWcs(md)
         poly = skypix.imageToPolygon(wcs,
                                      md.get("NAXIS1"), md.get("NAXIS2"),
                                      padRad=0.000075)  # about 15 arcsec
@@ -124,9 +121,10 @@ def processBand(filterDir, conn, done, qsp):
             conn.commit()
 
     conn.commit()
-    print(filterDir, \
-        "... %d processed, %d skipped, %d unrecognized" % \
-        (nProcessed, nSkipped, nUnrecognized), file=sys.stderr)
+    print(filterDir,
+          "... %d processed, %d skipped, %d unrecognized" %
+          (nProcessed, nSkipped, nUnrecognized), file=sys.stderr)
+
 
 if __name__ == "__main__":
     parser = OptionParser(usage="""%prog [options] DIR ...

bin.src/genInputRegistry.py

@@ -28,14 +28,11 @@
 import os
 import re
 import shutil
-try:
-    import sqlite3
-except ImportError:
-    # try external pysqlite package; deprecated
-    import sqlite as sqlite3
+import sqlite3
 import sys
 import lsst.daf.base as dafBase
-import lsst.afw.image as afwImage
+from lsst.afw.fits import readMetadata
+from lsst.afw.geom import makeSkyWcs
 import lsst.skypix as skypix
 
 
@@ -110,7 +107,7 @@ def processRun(runDir, conn, done, qsp):
             nSkipped += 1
             continue
 
-        md = afwImage.readMetadata(fits)
+        md = readMetadata(fits)
         date = md.get("DATE-OBS")
         if date.find("-") != -1:
             (year, month, day) = md.get("DATE-OBS").split("-")
@@ -134,7 +131,7 @@ def processRun(runDir, conn, done, qsp):
             id = row[0]
             break
 
-        wcs = afwImage.makeWcs(md)
+        wcs = makeSkyWcs(md)
         poly = skypix.imageToPolygon(wcs,
                                      md.get("NAXIS1"), md.get("NAXIS2"),
                                      padRad=0.000075)  # about 15 arcsec
@@ -148,9 +145,10 @@ def processRun(runDir, conn, done, qsp):
             conn.commit()
 
     conn.commit()
-    print(runDir, \
-        "... %d processed, %d skipped, %d unrecognized" % \
-        (nProcessed, nSkipped, nUnrecognized), file=sys.stderr)
+    print(runDir,
+          "... %d processed, %d skipped, %d unrecognized" %
+          (nProcessed, nSkipped, nUnrecognized), file=sys.stderr)
+
 
 if __name__ == "__main__":
     parser = OptionParser(usage="""%prog [options] DIR ...

policy/SdssMapper.yaml

@@ -82,7 +82,7 @@ datasets:
     template: '%(run)d/%(rerun)d/objcs/%(camcol)d/psField-%(run)06d-%(camcol)d-%(field)04d.fit'
   asTrans:
     persistable: ignored
-    python: lsst.afw.image.Wcs
+    python: lsst.afw.geom.skyWcs.SkyWcs
     storage: FitsStorage
     tables: raw
     template: '%(run)d/%(rerun)d/astrom/asTrans-%(run)06d.fit'

python/lsst/obs/sdss/convertasTrans.py

@@ -1,5 +1,4 @@
 from __future__ import absolute_import, division, print_function
-from builtins import map
 from builtins import range
 from builtins import object
 #
@@ -30,6 +29,7 @@
 
 import lsst.afw.image as afwImage
 import lsst.afw.geom as afwGeom
+from lsst.afw.geom import makeSkyWcs
 import lsst.afw.coord as afwCoord
 import lsst.afw.table as afwTable
 import lsst.meas.astrom.sip as sip
@@ -84,9 +84,9 @@ def __init__(self, node_rad, incl_rad, dRow0, dRow1, dRow2, dRow3, dCol0, dCol1,
         self.cOff = cOffset
 
     def xyToMuNu(self, x, y):
-        rowm     = (y+self.cOff) + self.dRow0 + self.dRow1*(x+self.cOff) + self.dRow2*((x+self.cOff)**2) + \
+        rowm = (y+self.cOff) + self.dRow0 + self.dRow1*(x+self.cOff) + self.dRow2*((x+self.cOff)**2) + \
             self.dRow3*((x+self.cOff)**3)
-        colm     = (x+self.cOff) + self.dCol0 + self.dCol1*(x+self.cOff) + self.dCol2*((x+self.cOff)**2) + \
+        colm = (x+self.cOff) + self.dCol0 + self.dCol1*(x+self.cOff) + self.dCol2*((x+self.cOff)**2) + \
             self.dCol3*((x+self.cOff)**3)
 
         mu_deg = self.a + self.b * rowm + self.c * colm
@@ -152,7 +152,7 @@ def createWcs(x, y, mapper, order=4, cOffset=1.0):
 
     # CRVAL1  = RA at Reference Pixel
     # CRVAL2  = DEC at Reference Pixel
-    crval = afwCoord.Coord(afwGeom.Point2D(ra_rad[0], dec_rad[0]), afwGeom.radians)
+    crval = afwCoord.IcrsCoord(afwGeom.Point2D(ra_rad[0], dec_rad[0]), afwGeom.radians)
 
     # CD1_1   = RA  degrees per column pixel
     # CD1_2   = RA  degrees per row pixel
@@ -162,16 +162,19 @@ def createWcs(x, y, mapper, order=4, cOffset=1.0):
     ULl = mapper.xyToRaDec(0., 1.)
     LRl = mapper.xyToRaDec(1., 0.)
 
-    LLc = afwCoord.Coord(afwGeom.Point2D(LLl[0], LLl[1]), afwGeom.radians)
-    ULc = afwCoord.Coord(afwGeom.Point2D(ULl[0], ULl[1]), afwGeom.radians)
-    LRc = afwCoord.Coord(afwGeom.Point2D(LRl[0], LRl[1]), afwGeom.radians)
+    LLc = afwCoord.IcrsCoord(afwGeom.Point2D(LLl[0], LLl[1]), afwGeom.radians)
+    ULc = afwCoord.IcrsCoord(afwGeom.Point2D(ULl[0], ULl[1]), afwGeom.radians)
+    LRc = afwCoord.IcrsCoord(afwGeom.Point2D(LRl[0], LRl[1]), afwGeom.radians)
 
     cdN_1 = LLc.getTangentPlaneOffset(LRc)
     cdN_2 = LLc.getTangentPlaneOffset(ULc)
     cd1_1, cd2_1 = cdN_1[0].asDegrees(), cdN_1[1].asDegrees()
     cd1_2, cd2_2 = cdN_2[0].asDegrees(), cdN_2[1].asDegrees()
 
-    linearWcs = afwImage.makeWcs(crval, crpix, cd1_1, cd2_1, cd1_2, cd2_2)
+    cdMatrix = np.array([cd1_1, cd2_1, cd1_2, cd2_2], dtype=float)
+    cdMatrix.shape = (2, 2)
+
+    linearWcs = makeSkyWcs(crpix=crpix, crval=crval, cdMatrix=cdMatrix)
     wcs = sip.makeCreateWcsWithSip(matches, linearWcs, order).getNewWcs()
 
     return wcs
@@ -181,7 +184,7 @@ def validate(xs, ys, mapper, wcs):
     dists = []
     for i in range(len(xs)):
         tuple1 = mapper.xyToRaDec(xs[i], ys[i])
-        coord1 = afwCoord.Coord(afwGeom.Point2D(tuple1[0], tuple1[1]), afwGeom.radians)
+        coord1 = afwCoord.IcrsCoord(afwGeom.Point2D(tuple1[0], tuple1[1]), afwGeom.radians)
         coord2 = wcs.pixelToSky(xs[i], ys[i])
         dist = coord1.angularSeparation(coord2).asArcseconds()
         dists.append(dist)
@@ -264,6 +267,7 @@ def convertasTrans(infile, filt, camcol, field, stepSize=50, doValidate=False):
 
     return wcs
 
+
 if __name__ == '__main__':
     infile = sys.argv[1]
     filt = sys.argv[2]

python/lsst/obs/sdss/convertfpM.py

@@ -27,7 +27,7 @@
 import os
 import re
 import pyfits
-import numpy as num
+
 import lsst.afw.image as afwImage
 import lsst.afw.geom as afwGeom
 
@@ -171,7 +171,7 @@ def convertfpM(infile, allPlanes=False):
 
     convertfpM(infile).writeFits(outfile)
 
-    comparison = """
+    comparison = """  # noqa ignore the really long line
 import lsst.afw.image as afwImage
 import numpy as num
 import os

python/lsst/obs/sdss/converttsField.py

@@ -79,6 +79,7 @@ def converttsField(infile, filt, exptime=53.907456):
         airmass=airmass,
     )
 
+
 if __name__ == '__main__':
     infile = sys.argv[1]
     filt = sys.argv[2]

python/lsst/obs/sdss/makeCamera.py

@@ -186,7 +186,6 @@ def makeCamera(name="SDSS", outputDir=None):
             ampInfoCatDict[k].writeFits(os.path.join(outputDir, "%s.fits"%(k)))
     return makeCameraFromCatalogs(camConfig, ampInfoCatDict)
 
-#************************************************************************************************************
 #
 # Print a Ccd
 #
@@ -233,12 +232,10 @@ def printCamera(title, camera):
     print(title, "Camera:", camera.getName())
 
     for det in camera:
-        print("%s %dx%d centre (mm): %s" % \
-            (det.getName(),
-             det.getBBox().getWidth(), det.getBBox().getHeight(),
-             det.getCenter(FOCAL_PLANE).getPoint()))
-
-#************************************************************************************************************
+        print("%s %dx%d centre (mm): %s" %
+              (det.getName(),
+               det.getBBox().getWidth(), det.getBBox().getHeight(),
+               det.getCenter(FOCAL_PLANE).getPoint()))
 
 
 def main():
@@ -254,5 +251,6 @@ def main():
     print()
     printCcd("Trimmed ", ccd, trimmed=True)
 
+
 if __name__ == "__main__":
     main()

python/lsst/obs/sdss/scaleSdssZeroPoint.py

@@ -2,7 +2,6 @@
 from builtins import zip
 from builtins import range
 from builtins import object
-#!/usr/bin/env python
 #
 # LSST Data Management System
 # Copyright 2008, 2009, 2010 LSST Corporation.

python/lsst/obs/sdss/selectFluxMag0.py

@@ -1,7 +1,6 @@
 from __future__ import absolute_import, division, print_function
 from builtins import str
 from builtins import range
-#!/usr/bin/env python
 #
 # LSST Data Management System
 # Copyright 2008, 2009, 2010 LSST Corporation.

python/lsst/obs/sdss/selectSdssImages.py

@@ -1,6 +1,5 @@
 from __future__ import absolute_import, division, print_function
 from builtins import range
-#!/usr/bin/env python
 #
 # LSST Data Management System
 # Copyright 2008, 2009, 2010 LSST Corporation.
@@ -180,6 +179,7 @@ def getColumnNames():
             "run rerun camcol field filter ra1 dec1 ra2 dec2 ra3 dec3 ra4 dec4".split() +
             "strip psfWidth sky airmass quality isblacklisted".split()
         )
+
     def __hash__(self):
         """
         Objects used in sets are required to be hashable, which requires (among other things) a
@@ -280,7 +280,7 @@ def run(self, coordList, filter, strip=None):
         exposureInfoList = [ExposureInfo(result) for result in cursor]
 
         runExpInfoSetDict = dict()
-        
+
         for expInfo in exposureInfoList:
             run = expInfo.dataId["run"]
             expInfoSet = runExpInfoSetDict.get(run)

python/lsst/obs/sdss/yanny.py

@@ -511,8 +511,8 @@ def isarray(self, structure, variable):
             typ = self.type(structure, variable)
             character_array = re.compile(r'char[[<]\d*[]>][[<]\d*[]>]')
             if ((character_array.search(typ) is not None) or
-                (typ.find('char') < 0 and (typ.find('[') >= 0
-                                           or typ.find('<') >= 0))):
+                (typ.find('char') < 0 and (typ.find('[') >= 0 or
+                                           typ.find('<') >= 0))):
                 cache[variable] = True
             else:
                 cache[variable] = False
@@ -1118,7 +1118,7 @@ def _parse(self):
                 #
                 line = line.strip()
                 line = self.trailing_comment(line)
-                #line = trailing_comments.sub('',line)
+                # line = trailing_comments.sub('',line)
                 line = double_braces.sub('""', line)
                 #
                 # Now if the first word on the line does not match a

tests/test_sdss.py

@@ -30,6 +30,9 @@
 import lsst.daf.persistence as dafPersist
 from lsst.daf.base import DateTime
 import lsst.afw.image
+from lsst.afw.coord import IcrsCoord
+from lsst.afw.geom import SkyWcs
+from lsst.afw.geom import degrees
 import lsst.afw.detection
 
 
@@ -68,10 +71,12 @@ def testGetDR7(self):
             self.assertEqual(h, 31)
 
             wcs = ref.get("asTrans")
-            self.assertEqual(wcs.__class__, lsst.afw.image.TanWcs)
-            self.assertFalse(wcs.isFlipped())
-            self.assertAlmostEqual(wcs.getFitsMetadata().get("CRPIX1"), 1.0, 5)
-            self.assertAlmostEqual(wcs.getFitsMetadata().get("CRPIX2"), 1.0, 5)
+            self.assertIsInstance(wcs, SkyWcs)
+            self.assertFalse(wcs.isFlipped)
+            # comparison is to results from lsst.afw.image.TanWcs class
+            self.assertCoordsAlmostEqual(wcs.pixelToSky(700, 1000),
+                                         IcrsCoord(343.6507738304687 * degrees,
+                                                   -0.3509870420713227 * degrees))
 
             tsField = ref.get("tsField")
             self.assertAlmostEqual(tsField.gain, 4.72, 2)
@@ -90,6 +95,7 @@ class TestMemory(lsst.utils.tests.MemoryTestCase):
 def setup_module(module):
     lsst.utils.tests.init()
 
+
 if __name__ == "__main__":
     lsst.utils.tests.init()
     unittest.main()

tests/test_selectFluxMag0.py

@@ -89,7 +89,7 @@ def makeTestExposure(self, xNumPix=2060, yNumPix=1967):
         metadata.set("LTV1", -341970)
         metadata.set("LTV2", -11412)
         # exposure needs a wcs and a bbox
-        wcs = afwImage.makeWcs(metadata)
+        wcs = afwGeom.makeSkyWcs(metadata)
         bbox = afwGeom.Box2I(afwGeom.Point2I(341970, 11412), afwGeom.Extent2I(xNumPix, yNumPix))
         exposure = afwImage.ExposureF(bbox, wcs)
         mi = exposure.getMaskedImage()
@@ -165,6 +165,7 @@ class TestMemory(lsst.utils.tests.MemoryTestCase):
 def setup_module(module):
     lsst.utils.tests.init()
 
+
 if __name__ == "__main__":
     lsst.utils.tests.init()
     unittest.main()

tests/test_selectSdssImages.py

@@ -365,6 +365,7 @@ class TestMemory(lsst.utils.tests.MemoryTestCase):
 def setup_module(module):
     lsst.utils.tests.init()
 
+
 if __name__ == "__main__":
     lsst.utils.tests.init()
     unittest.main()