Merge pull request #37 from lsst/tickets/DM-19887

DM-19887: Use geom rather than afwGeom

examples/lsstSkyMap.py

@@ -22,7 +22,7 @@
 #
 import math
 
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 from lsst.skymap import DodecaSkyMap
 
 print("Warning: this does not take into account the extra space required by patch borders")
@@ -41,21 +41,21 @@
         numPix = dimensions[0] * dimensions[1]
         totNumPix += numPix
         wcs = tractInfo.getWcs()
-        posBBox = afwGeom.Box2D(bbox)
+        posBBox = geom.Box2D(bbox)
         ctrPixPos = posBBox.getCenter()
         ctrCoord = wcs.pixelToSky(ctrPixPos)
-        ctrSkyPosDeg = ctrCoord.getPosition(afwGeom.degrees)
-        leftCoord = wcs.pixelToSky(afwGeom.Point2D(posBBox.getMinX(), ctrPixPos[1]))
-        rightCoord = wcs.pixelToSky(afwGeom.Point2D(posBBox.getMaxX(), ctrPixPos[1]))
-        topCoord = wcs.pixelToSky(afwGeom.Point2D(ctrPixPos[0], posBBox.getMinY()))
-        bottomCoord = wcs.pixelToSky(afwGeom.Point2D(ctrPixPos[0], posBBox.getMaxY()))
+        ctrSkyPosDeg = ctrCoord.getPosition(geom.degrees)
+        leftCoord = wcs.pixelToSky(geom.Point2D(posBBox.getMinX(), ctrPixPos[1]))
+        rightCoord = wcs.pixelToSky(geom.Point2D(posBBox.getMaxX(), ctrPixPos[1]))
+        topCoord = wcs.pixelToSky(geom.Point2D(ctrPixPos[0], posBBox.getMinY()))
+        bottomCoord = wcs.pixelToSky(geom.Point2D(ctrPixPos[0], posBBox.getMaxY()))
         xSpan = leftCoord.separation(rightCoord).asDegrees()
         ySpan = bottomCoord.separation(topCoord).asDegrees()
         print("%3d   %7.1f %7.1f %10.2e  %10.2e %10.1e %6.1f %6.1f" %
               (tractInfo.getId(), ctrSkyPosDeg[0], ctrSkyPosDeg[1],
                dimensions[0], dimensions[1], numPix, xSpan, ySpan))
 
-    pixelScaleRad = afwGeom.Angle(skyMap.config.pixelScale, afwGeom.arcseconds).asRadians()
+    pixelScaleRad = geom.Angle(skyMap.config.pixelScale, geom.arcseconds).asRadians()
     nomPixelAreaRad2 = pixelScaleRad**2  # nominal area of a pixel in rad^2
     numPixToTileSphere = 4 * math.pi / nomPixelAreaRad2
     print("total pixels = %.1e" % (totNumPix,))

examples/plotDodecaSkyMap.py

@@ -28,7 +28,7 @@
 from mpl_toolkits.mplot3d import Axes3D  # noqa F401 used by fig.gca
 import matplotlib.pyplot as plt
 
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 from lsst.skymap import DodecaSkyMap
 
 skyMap = DodecaSkyMap()
@@ -54,9 +54,9 @@
     bbox = tractInfo.getBBox()
     outerPixPos = [
         bbox.getMin(),
-        afwGeom.Point2I(bbox.getMaxX(), bbox.getMinY()),
+        geom.Point2I(bbox.getMaxX(), bbox.getMinY()),
         bbox.getMax(),
-        afwGeom.Point2I(bbox.getMinX(), bbox.getMaxY()),
+        geom.Point2I(bbox.getMinX(), bbox.getMaxY()),
         bbox.getMin(),
     ]
     outerPoints = [numpy.array(wcs.pixelToSky(p[0], p[1]).getVector()) for p in outerPixPos]

examples/plotSkyMap.py

@@ -31,7 +31,7 @@
 from mpl_toolkits.mplot3d import Axes3D  # noqa F401 used by fig.gca
 import matplotlib.pyplot as plt
 
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 
 
 def reportSkyMapInfo(skyMap):
@@ -45,14 +45,14 @@ def reportSkyMapInfo(skyMap):
     print("\nSkyMap has %d tracts:" % (len(skyMap)))
     for tractInfo in skyMap:
         wcs = tractInfo.getWcs()
-        posBox = afwGeom.Box2D(tractInfo.getBBox())
+        posBox = geom.Box2D(tractInfo.getBBox())
         pixelPosList = (
             posBox.getMin(),
-            afwGeom.Point2D(posBox.getMaxX(), posBox.getMinY()),
+            geom.Point2D(posBox.getMaxX(), posBox.getMinY()),
             posBox.getMax(),
-            afwGeom.Point2D(posBox.getMinX(), posBox.getMaxY()),
+            geom.Point2D(posBox.getMinX(), posBox.getMaxY()),
         )
-        skyPosList = [wcs.pixelToSky(pos).getPosition(afwGeom.degrees) for pos in pixelPosList]
+        skyPosList = [wcs.pixelToSky(pos).getPosition(geom.degrees) for pos in pixelPosList]
         posStrList = ["(%0.3f, %0.3f)" % tuple(skyPos) for skyPos in skyPosList]
         print("tract %s has corners %s (RA, Dec deg) and %s x %s patches" %
               (tractInfo.getId(), ", ".join(posStrList),
@@ -73,7 +73,7 @@ def plotSkyMap3d(skyMap):
     for tractInfo in skyMap:
         # display outer edge; scale to be approximately in the same plane as the inner region
         wcs = tractInfo.getWcs()
-        posBox = afwGeom.Box2D(tractInfo.getBBox())
+        posBox = geom.Box2D(tractInfo.getBBox())
         xRange = posBox.getMinX(), posBox.getMaxX()
         yRange = posBox.getMinY(), posBox.getMaxY()
 
@@ -181,7 +181,7 @@ def plotSkyMap2d(skyMap):
                            (xList, yMax*numpy.ones(num)),
                            (xMin*numpy.ones(num), yList),
                            ):
-                coords = [wcs.pixelToSky(afwGeom.Point2D(x1, y1)) for x1, y1 in zip(xs, ys)]
+                coords = [wcs.pixelToSky(geom.Point2D(x1, y1)) for x1, y1 in zip(xs, ys)]
                 bounds = [proj.projectWithRecenter(c) for c in coords]
                 axes.plot([b[0] for b in bounds], [b[1] for b in bounds], color + '-')
 

examples/showVisitSkyMap.py

@@ -26,6 +26,7 @@
 import matplotlib.pyplot as pyplot
 
 import lsst.afw.cameraGeom as cameraGeom
+import lsst.geom as geom
 import lsst.afw.geom as afwGeom
 import lsst.daf.persistence as dafPersist
 from lsst.pipe.base.argumentParser import IdValueAction, DataIdContainer
@@ -35,7 +36,7 @@ def bboxToRaDec(bbox, wcs):
     """Get the corners of a BBox and convert them to lists of RA and Dec."""
     corners = []
     for corner in bbox.getCorners():
-        p = afwGeom.Point2D(corner.getX(), corner.getY())
+        p = geom.Point2D(corner.getX(), corner.getY())
         coord = wcs.pixelToSky(p)
         corners.append([coord.getRa().asDegrees(), coord.getDec().asDegrees()])
     ra, dec = zip(*corners)
@@ -91,10 +92,10 @@ def main(rootDir, tracts, visits, ccds=None, ccdKey='ccd', showPatch=False, save
 
                     # add CCD serial numbers
                     if showCcds:
-                        minPoint = afwGeom.Point2D(min(ra), min(dec))
-                        maxPoint = afwGeom.Point2D(max(ra), max(dec))
+                        minPoint = geom.Point2D(min(ra), min(dec))
+                        maxPoint = geom.Point2D(max(ra), max(dec))
                         # Use doubles in Box2D to check overlap
-                        bboxDouble = afwGeom.Box2D(minPoint, maxPoint)
+                        bboxDouble = geom.Box2D(minPoint, maxPoint)
                         overlaps = [not bboxDouble.overlaps(otherBbox) for otherBbox in bboxesPlotted]
                         if all(overlaps):
                             pyplot.text(percent(ra), percent(dec), str(ccdId), fontsize=6,

python/lsst/skymap/detail/utils.py

@@ -25,7 +25,7 @@
 import numpy
 
 import lsst.sphgeom
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 
 
 _TinyFloat = numpy.finfo(float).tiny
@@ -54,5 +54,5 @@ def coordFromVec(vec, defRA=None):
             decDeg = 90.0
         else:
             decDeg = -90.0
-        return afwGeom.SpherePoint(defRA, decDeg*afwGeom.degrees)
-    return afwGeom.SpherePoint(lsst.sphgeom.Vector3d(*vec))
+        return geom.SpherePoint(defRA, decDeg*geom.degrees)
+    return geom.SpherePoint(lsst.sphgeom.Vector3d(*vec))

python/lsst/skymap/detail/wcsFactory.py

@@ -22,6 +22,7 @@
 
 __all__ = ['WcsFactory']
 
+import lsst.geom as geom
 import lsst.afw.geom as afwGeom
 
 
@@ -43,7 +44,7 @@ class WcsFactory:
         Flip the X axis?
     """
 
-    def __init__(self, pixelScale, projection, rotation=0*afwGeom.radians, flipX=False):
+    def __init__(self, pixelScale, projection, rotation=0*geom.radians, flipX=False):
         if len(projection) != 3:
             raise RuntimeError("projection=%r; must have length 3" % (projection,))
         self._projection = projection

python/lsst/skymap/discreteSkyMap.py

@@ -25,7 +25,7 @@
 import struct
 
 from lsst.pex.config import ListField
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 from .cachingSkyMap import CachingSkyMap
 from .tractInfo import ExplicitTractInfo
 
@@ -69,11 +69,11 @@ def __init__(self, config, version=0):
 
     def generateTract(self, index):
         """Generate TractInfo for the specified tract index."""
-        center = afwGeom.SpherePoint(self.config.raList[index], self.config.decList[index], afwGeom.degrees)
+        center = geom.SpherePoint(self.config.raList[index], self.config.decList[index], geom.degrees)
         radius = self.config.radiusList[index]
-        wcs = self._wcsFactory.makeWcs(crPixPos=afwGeom.Point2D(0, 0), crValCoord=center)
+        wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=center)
         return ExplicitTractInfo(index, self.config.patchInnerDimensions, self.config.patchBorder, center,
-                                 radius * afwGeom.degrees, self.config.tractOverlap * afwGeom.degrees, wcs)
+                                 radius*geom.degrees, self.config.tractOverlap*geom.degrees, wcs)
 
     def updateSha1(self, sha1):
         """Add subclass-specific state or configuration options to the SHA1."""

python/lsst/skymap/dodecaSkyMap.py

@@ -30,7 +30,7 @@
 import struct
 
 import lsst.pex.config as pexConfig
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 from . import detail
 from .baseSkyMap import BaseSkyMap
 from .tractInfo import TractInfo
@@ -69,16 +69,16 @@ def __init__(self, config=None):
         BaseSkyMap.__init__(self, config)
         self._dodecahedron = detail.Dodecahedron(withFacesOnPoles=self.config.withTractsOnPoles)
 
-        tractOverlap = afwGeom.Angle(self.config.tractOverlap, afwGeom.degrees)
+        tractOverlap = geom.Angle(self.config.tractOverlap, geom.degrees)
 
         for id in range(12):
             tractVec = self._dodecahedron.getFaceCtr(id)
-            tractCoord = detail.coordFromVec(tractVec, defRA=afwGeom.Angle(0))
+            tractCoord = detail.coordFromVec(tractVec, defRA=geom.Angle(0))
             tractRA = tractCoord.getLongitude()
             vertexVecList = self._dodecahedron.getVertices(id)
 
             # make initial WCS; don't worry about crPixPos because TractInfo will shift it as required
-            wcs = self._wcsFactory.makeWcs(crPixPos=afwGeom.Point2D(0, 0), crValCoord=tractCoord)
+            wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=tractCoord)
 
             self._tractInfoList.append(
                 TractInfo(

python/lsst/skymap/equatSkyMap.py

@@ -25,7 +25,7 @@
 import struct
 
 import lsst.pex.config as pexConfig
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 from .baseSkyMap import BaseSkyMap
 from .tractInfo import TractInfo
 
@@ -64,29 +64,29 @@ class EquatSkyMap(BaseSkyMap):
     def __init__(self, config=None):
         BaseSkyMap.__init__(self, config)
 
-        decRange = tuple(afwGeom.Angle(dr, afwGeom.degrees) for dr in self.config.decRange)
+        decRange = tuple(geom.Angle(dr, geom.degrees) for dr in self.config.decRange)
         midDec = (decRange[0] + decRange[1]) / 2.0
-        tractWidthRA = afwGeom.Angle(360.0 / self.config.numTracts, afwGeom.degrees)
-        tractOverlap = afwGeom.Angle(self.config.tractOverlap, afwGeom.degrees)
+        tractWidthRA = geom.Angle(360.0 / self.config.numTracts, geom.degrees)
+        tractOverlap = geom.Angle(self.config.tractOverlap, geom.degrees)
 
         for id in range(self.config.numTracts):
             begRA = tractWidthRA * id
             endRA = begRA + tractWidthRA
             vertexCoordList = (
-                afwGeom.SpherePoint(begRA, decRange[0]),
-                afwGeom.SpherePoint(endRA, decRange[0]),
-                afwGeom.SpherePoint(endRA, decRange[1]),
-                afwGeom.SpherePoint(begRA, decRange[1]),
+                geom.SpherePoint(begRA, decRange[0]),
+                geom.SpherePoint(endRA, decRange[0]),
+                geom.SpherePoint(endRA, decRange[1]),
+                geom.SpherePoint(begRA, decRange[1]),
             )
 
             midRA = begRA + tractWidthRA / 2.0
-            ctrCoord = afwGeom.SpherePoint(midRA, midDec)
+            ctrCoord = geom.SpherePoint(midRA, midDec)
 
             # CRVal must have Dec=0 for symmetry about the equator
-            crValCoord = afwGeom.SpherePoint(midRA, afwGeom.Angle(0.0))
+            crValCoord = geom.SpherePoint(midRA, geom.Angle(0.0))
 
             # make initial WCS; don't worry about crPixPos because TractInfo will shift it as required
-            wcs = self._wcsFactory.makeWcs(crPixPos=afwGeom.Point2D(0, 0), crValCoord=crValCoord)
+            wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=crValCoord)
 
             self._tractInfoList.append(TractInfo(
                 id=id,

python/lsst/skymap/healpixSkyMap.py

@@ -40,7 +40,7 @@ def __getattr__(self, name, e=e):
     healpy = DummyHealpy()
 
 from lsst.pex.config import Field
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 from .cachingSkyMap import CachingSkyMap
 from .tractInfo import TractInfo
 
@@ -52,7 +52,7 @@ def angToCoord(thetaphi):
     of healpy functions can be directed to this function without
     additional translation.
     """
-    return afwGeom.SpherePoint(float(thetaphi[1]), float(thetaphi[0] - 0.5*numpy.pi), afwGeom.radians)
+    return geom.SpherePoint(float(thetaphi[1]), float(thetaphi[0] - 0.5*numpy.pi), geom.radians)
 
 
 def coordToAng(coord):
@@ -126,9 +126,9 @@ def findTract(self, coord):
     def generateTract(self, index):
         """Generate TractInfo for the specified tract index."""
         center = angToCoord(healpy.pix2ang(self._nside, index, nest=self.config.nest))
-        wcs = self._wcsFactory.makeWcs(crPixPos=afwGeom.Point2D(0, 0), crValCoord=center)
+        wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=center)
         return HealpixTractInfo(self._nside, index, self.config.nest, self.config.patchInnerDimensions,
-                                self.config.patchBorder, center, self.config.tractOverlap*afwGeom.degrees,
+                                self.config.patchBorder, center, self.config.tractOverlap*geom.degrees,
                                 wcs)
 
     def updateSha1(self, sha1):

python/lsst/skymap/patchInfo.py

@@ -23,7 +23,7 @@
 __all__ = ["PatchInfo", "makeSkyPolygonFromBBox"]
 
 from lsst.sphgeom import ConvexPolygon
-from lsst.afw.geom import Box2D
+from lsst.geom import Box2D
 
 
 def makeSkyPolygonFromBBox(bbox, wcs):

python/lsst/skymap/ringsSkyMap.py

@@ -26,7 +26,7 @@
 import math
 
 from lsst.pex.config import Field
-import lsst.afw.geom as afwGeom
+import lsst.geom as geom
 from .cachingSkyMap import CachingSkyMap
 from .tractInfo import ExplicitTractInfo
 
@@ -87,7 +87,7 @@ def __init__(self, config, version=1):
             self._ringNums.append(int(2*math.pi*math.cos(dec)/self._ringSize) + 1)
         numTracts = sum(self._ringNums) + 2
         super(RingsSkyMap, self).__init__(numTracts, config, version)
-        self._raStart = self.config.raStart*afwGeom.degrees
+        self._raStart = self.config.raStart*geom.degrees
 
     def getRingIndices(self, index):
         """Calculate ring indices given a numerical index of a tract.
@@ -129,13 +129,13 @@ def generateTract(self, index):
             ra, dec = 0, 0.5*math.pi
         else:
             dec = self._ringSize*(ringNum + 1) - 0.5*math.pi
-            ra = ((2*math.pi*tractNum/self._ringNums[ringNum])*afwGeom.radians +
+            ra = ((2*math.pi*tractNum/self._ringNums[ringNum])*geom.radians +
                   self._raStart).wrap().asRadians()
 
-        center = afwGeom.SpherePoint(ra, dec, afwGeom.radians)
-        wcs = self._wcsFactory.makeWcs(crPixPos=afwGeom.Point2D(0, 0), crValCoord=center)
+        center = geom.SpherePoint(ra, dec, geom.radians)
+        wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=center)
         return ExplicitTractInfo(index, self.config.patchInnerDimensions, self.config.patchBorder, center,
-                                 0.5*self._ringSize*afwGeom.radians, self.config.tractOverlap*afwGeom.degrees,
+                                 0.5*self._ringSize*geom.radians, self.config.tractOverlap*geom.degrees,
                                  wcs)
 
     def _decToRingNum(self, dec):