DM-32402: Add TractBuilders and sub-patch cell capability.

python/lsst/skymap/__init__.py

@@ -1,3 +1,5 @@
+from .cellInfo import *
+from .tractBuilder import *
 from .patchInfo import *
 from .tractInfo import *
 from .baseSkyMap import *
@@ -6,3 +8,4 @@
 from .discreteSkyMap import *
 from .skyMapRegistry import *
 from .version import *
+from .detail import makeSkyPolygonFromBBox, Index2D

python/lsst/skymap/baseSkyMap.py

@@ -27,27 +27,21 @@
 __all__ = ["BaseSkyMapConfig", "BaseSkyMap"]
 
 import hashlib
-import struct
 import numpy as np
 
 import lsst.geom as geom
 import lsst.pex.config as pexConfig
 from lsst.geom import SpherePoint, Angle, arcseconds, degrees
 from . import detail
+from .tractBuilder import tractBuilderRegistry
 
 
 class BaseSkyMapConfig(pexConfig.Config):
-    patchInnerDimensions = pexConfig.ListField(
-        doc="dimensions of inner region of patches (x,y pixels)",
-        dtype=int,
-        length=2,
-        default=(4000, 4000),
-    )
-    patchBorder = pexConfig.Field(
-        doc="border between patch inner and outer bbox (pixels)",
-        dtype=int,
-        default=100,
+    tractBuilder = tractBuilderRegistry.makeField(
+        doc="Algorithm for creating patches within the tract.",
+        default="legacy"
     )
+
     tractOverlap = pexConfig.Field(
         doc="minimum overlap between adjacent sky tracts, on the sky (deg)",
         dtype=float,
@@ -72,6 +66,64 @@ class BaseSkyMapConfig(pexConfig.Config):
         default=0,
     )
 
+    # Backwards compatibility
+    # We can't use the @property decorator because it makes pexConfig sad.
+    def getPatchInnerDimensions(self):
+        """Get the patch inner dimensions, for backwards compatibility.
+
+        This value is only used with the ``legacy`` tract builder,
+        and is ignored otherwise.  In general, the config should be
+        accessed directly with config.tractBuilder["legacy"].patchInnerDimensions.
+
+        Returns
+        -------
+        innerDimensions : `list` [`int`, `int`]
+        """
+        return self.tractBuilder["legacy"].patchInnerDimensions
+
+    def setPatchInnerDimensions(self, value):
+        """Set the patch inner dimensions, for backwards compatibility.
+
+        This value is only used with the ``legacy`` tract builder,
+        and is ignored otherwise.  In general, the config should be
+        accessed directly with config.tractBuilder["legacy"].patchInnerDimensions.
+
+        Parameters
+        ----------
+        value : `list` [`int`, `int`]
+        """
+        self.tractBuilder["legacy"].patchInnerDimensions = value
+
+    patchInnerDimensions = property(getPatchInnerDimensions, setPatchInnerDimensions)
+
+    def getPatchBorder(self):
+        """Get the patch border, for backwards compatibility.
+
+        This value is only used with the ``legacy`` tract builder,
+        and is ignored otherwise.  In general, the config should be
+        accessed directly with config.tractBuilder["legacy"].patchBorder.
+
+        Returns
+        -------
+        border: `int`
+        """
+        return self.tractBuilder["legacy"].patchBorder
+
+    def setPatchBorder(self, value):
+        """Set the patch border, for backwards compatibility.
+
+        This value is only used with the ``legacy`` tract builder,
+        and is ignored otherwise.  In general, the config should be
+        accessed directly with config.tractBuilder["legacy"].patchBorder.
+
+        Parameters
+        -------
+        border: `int`
+        """
+        self.tractBuilder["legacy"].patchBorder = value
+
+    patchBorder = property(getPatchBorder, setPatchBorder)
+
 
 class BaseSkyMap:
     """A collection of overlapping Tracts that map part or all of the sky.
@@ -113,6 +165,7 @@ def __init__(self, config=None):
             rotation=Angle(self.config.rotation, degrees),
         )
         self._sha1 = None
+        self._tractBuilder = config.tractBuilder.apply()
 
     def findTract(self, coord):
         """Find the tract whose center is nearest the specified coord.
@@ -298,16 +351,7 @@ def getSha1(self):
         if self._sha1 is None:
             sha1 = hashlib.sha1()
             sha1.update(type(self).__name__.encode('utf-8'))
-            configPacked = struct.pack(
-                "<iiidd3sd",
-                self.config.patchInnerDimensions[0],
-                self.config.patchInnerDimensions[1],
-                self.config.patchBorder,
-                self.config.tractOverlap,
-                self.config.pixelScale,
-                self.config.projection.encode('ascii'),
-                self.config.rotation
-            )
+            configPacked = self._tractBuilder.getPackedConfig(self.config)
             sha1.update(configPacked)
             self.updateSha1(sha1)
             self._sha1 = sha1.digest()

python/lsst/skymap/cellInfo.py

@@ -0,0 +1,168 @@
+#
+# LSST Data Management System
+# Copyright 2008, 2009, 2010 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/>.
+#
+
+__all__ = ["CellInfo"]
+
+from .detail import makeSkyPolygonFromBBox
+
+
+class CellInfo:
+    """Information about a cell within a patch of a tract of a sky map.
+
+    See `PatchInfo` and `TractInfo` for more information.
+
+    Parameters
+    ----------
+    index : `lsst.skymap.Index2D`
+        x,y index of a cell (a pair of ints)
+    innerBBox : `lsst.geom.Box2I`
+        Inner bounding box.
+    outerBBox : `lsst.geom.Box2I`
+        Outer bounding box.
+    sequentialIndex : `int`
+        Cell sequential index.
+    tractWcs : `lsst.afw.geom.SkyWcs`
+        Tract WCS object.
+    """
+    def __init__(self, index, innerBBox, outerBBox, sequentialIndex, tractWcs):
+        self._index = index
+        self._sequentialIndex = sequentialIndex
+        self._innerBBox = innerBBox
+        self._outerBBox = outerBBox
+        self._wcs = tractWcs
+        if not outerBBox.contains(innerBBox):
+            raise RuntimeError("outerBBox=%s does not contain innerBBox=%s" % (outerBBox, innerBBox))
+
+    def getIndex(self):
+        """Return cell index: a tuple of (x, y)
+
+        Returns
+        -------
+        result : `lsst.skymap.Index2D`
+            Patch index (x, y).
+        """
+        return self._index
+
+    index = property(getIndex)
+
+    def getSequentialIndex(self):
+        """Return cell sequential index.
+
+        Returns
+        -------
+        result : `int`
+            Sequential cell index.
+        """
+        return self._sequentialIndex
+
+    sequential_index = property(getSequentialIndex)
+
+    def getWcs(self):
+        """Return the associated tract wcs
+
+        Returns
+        -------
+        wcs : `lsst.afw.geom.SkyWcs`
+            Tract WCS.
+        """
+        return self._wcs
+
+    wcs = property(getWcs)
+
+    def getInnerBBox(self):
+        """Get inner bounding box.
+
+        Returns
+        -------
+        bbox : `lsst.geom.Box2I`
+            The inner bounding Box.
+        """
+        return self._innerBBox
+
+    inner_bbox = property(getInnerBBox)
+
+    def getOuterBBox(self):
+        """Get outer bounding box.
+
+        Returns
+        -------
+        bbox : `lsst.geom.Box2I`
+            The outer bounding Box.
+        """
+        return self._outerBBox
+
+    outer_bbox = property(getOuterBBox)
+
+    def getInnerSkyPolygon(self, tractWcs=None):
+        """Get the inner on-sky region.
+
+        Parameters
+        ----------
+        tractWcs : `lsst.afw.image.SkyWcs`, optional
+            WCS for the associated tract.
+
+        Returns
+        -------
+        result : `lsst.sphgeom.ConvexPolygon`
+            The inner sky region.
+        """
+        _tractWcs = tractWcs if tractWcs is not None else self._wcs
+        return makeSkyPolygonFromBBox(bbox=self.getInnerBBox(), wcs=_tractWcs)
+
+    @property
+    def inner_sky_polygon(self):
+        return self.getInnerSkyPolygon()
+
+    def getOuterSkyPolygon(self, tractWcs=None):
+        """Get the outer on-sky region.
+
+        Parameters
+        ----------
+        tractWcs : `lsst.afw.image.SkyWcs`, optional
+            WCS for the associated tract.
+
+        Returns
+        -------
+        result : `lsst.sphgeom.ConvexPolygon`
+            The outer sky region.
+        """
+        _tractWcs = tractWcs if tractWcs is not None else self._wcs
+        return makeSkyPolygonFromBBox(bbox=self.getOuterBBox(), wcs=_tractWcs)
+
+    @property
+    def outer_sky_polygon(self):
+        return self.getOuterSkyPolygon()
+
+    def __eq__(self, rhs):
+        return (self.getIndex() == rhs.getIndex()) \
+            and (self.getInnerBBox() == rhs.getInnerBBox()) \
+            and (self.getOuterBBox() == rhs.getOuterBBox())
+
+    def __ne__(self, rhs):
+        return not self.__eq__(rhs)
+
+    def __str__(self):
+        return "CellInfo(index=%s)" % (self.getIndex(),)
+
+    def __repr__(self):
+        return "CellInfo(index=%s, innerBBox=%s, outerBBox=%s)" % \
+            (self.getIndex(), self.getInnerBBox(), self.getOuterBBox())

python/lsst/skymap/detail/utils.py

@@ -20,14 +20,14 @@
 # see <http://www.lsstcorp.org/LegalNotices/>.
 #
 
-__all__ = ["coordFromVec"]
+__all__ = ["coordFromVec", "makeSkyPolygonFromBBox", "Index2D"]
 
+from typing import NamedTuple
 import numpy
 
 import lsst.sphgeom
 import lsst.geom as geom
 
-
 _TinyFloat = numpy.finfo(float).tiny
 
 
@@ -56,3 +56,38 @@ def coordFromVec(vec, defRA=None):
             decDeg = -90.0
         return geom.SpherePoint(defRA, decDeg*geom.degrees)
     return geom.SpherePoint(lsst.sphgeom.Vector3d(*vec))
+
+
+def makeSkyPolygonFromBBox(bbox, wcs):
+    """Make an on-sky polygon from a bbox and a SkyWcs
+
+    Parameters
+    ----------
+    bbox : `lsst.geom.Box2I` or `lsst.geom.Box2D`
+        Bounding box of region, in pixel coordinates
+    wcs : `lsst.afw.geom.SkyWcs`
+        Celestial WCS
+
+    Returns
+    -------
+    polygon : `lsst.sphgeom.ConvexPolygon`
+        On-sky region
+    """
+    pixelPoints = geom.Box2D(bbox).getCorners()
+    skyPoints = wcs.pixelToSky(pixelPoints)
+    return lsst.sphgeom.ConvexPolygon([sp.getVector() for sp in skyPoints])
+
+
+class Index2D(NamedTuple):
+    """Two dimensional index for patches and cells.
+
+    This class contains the x and y values of the location of a patch
+    within a tract, or a cell within a patch.
+
+    Parameters
+    ----------
+    x : `int`
+    y : `int`
+    """
+    x: int
+    y: int

python/lsst/skymap/discreteSkyMap.py

@@ -72,7 +72,7 @@ def generateTract(self, index):
         center = geom.SpherePoint(self.config.raList[index], self.config.decList[index], geom.degrees)
         radius = self.config.radiusList[index]
         wcs = self._wcsFactory.makeWcs(crPixPos=geom.Point2D(0, 0), crValCoord=center)
-        return ExplicitTractInfo(index, self.config.patchInnerDimensions, self.config.patchBorder, center,
+        return ExplicitTractInfo(index, self._tractBuilder, center,
                                  radius*geom.degrees, self.config.tractOverlap*geom.degrees, wcs)
 
     def updateSha1(self, sha1):

python/lsst/skymap/dodecaSkyMap.py

@@ -83,8 +83,7 @@ def __init__(self, config=None):
             self._tractInfoList.append(
                 TractInfo(
                     id=id,
-                    patchInnerDimensions=self.config.patchInnerDimensions,
-                    patchBorder=self.config.patchBorder,
+                    tractBuilder=self._tractBuilder,
                     ctrCoord=tractCoord,
                     vertexCoordList=[detail.coordFromVec(vec, defRA=tractRA) for vec in vertexVecList],
                     tractOverlap=tractOverlap,

python/lsst/skymap/equatSkyMap.py

@@ -90,8 +90,7 @@ def __init__(self, config=None):
 
             self._tractInfoList.append(TractInfo(
                 id=id,
-                patchInnerDimensions=self.config.patchInnerDimensions,
-                patchBorder=self.config.patchBorder,
+                tractBuilder=self._tractBuilder,
                 ctrCoord=ctrCoord,
                 vertexCoordList=vertexCoordList,
                 tractOverlap=tractOverlap,