DM-48509: Move tract and patch corner calculation in analysis_tools to utility file

python/lsst/analysis/tools/__init__.py

@@ -36,5 +36,6 @@
 """
 
 from .math import *
+from .utils import *
 from .version import *  # Generated by sconsUtils
 from .warning_control import *

python/lsst/analysis/tools/tasks/makeMetricTable.py

@@ -33,6 +33,8 @@
 from lsst.pipe.base import connectionTypes as ct
 from lsst.skymap import BaseSkyMap
 
+from ..utils import getTractCorners
+
 
 class MakeMetricTableConnections(
     pipeBase.PipelineTaskConnections,
@@ -181,7 +183,7 @@ def run(self, dataIdInfo, metricBundles, skymap):
 
         # Add tract corners if inputs are at the tract-level
         if "tract" in self.config.inputDataDimensions:
-            corners = self.getTractCorners(skymap, dataIdInfo[0]["tract"])
+            corners = getTractCorners(skymap, dataIdInfo[0]["tract"])
             metricsDict["corners"] = [corners]
 
         # Add the metrics and units from the first bundle to the dicts
@@ -219,7 +221,7 @@ def run(self, dataIdInfo, metricBundles, skymap):
                 metricsDict[key].append(value)
 
             if "tract" in self.config.inputDataDimensions:
-                corners = self.getTractCorners(skymap, dataIdInfo[0]["tract"])
+                corners = getTractCorners(skymap, dataIdInfo[0]["tract"])
                 metricsDict["corners"].append(corners)
 
             metricNames = list(metricsDict)
@@ -248,35 +250,3 @@ def run(self, dataIdInfo, metricBundles, skymap):
 
         metricTableStruct = pipeBase.Struct(metricTable=Table(metricsDict, units=metricUnits))
         return metricTableStruct
-
-    def getTractCorners(self, skymap, tract):
-        """Calculate the corners of a tract, given  skymap.
-
-        Parameters
-        ----------
-        skymap : `lsst.skymap`
-        tract : `int`
-
-        Returns
-        -------
-        corners : `list` of `tuples` of `float`
-
-        Notes
-        -----
-        Corners are returned in degrees and wrapped in ra.
-        """
-        # Find the tract corners
-        tractCorners = skymap[tract].getVertexList()
-        corners = [(corner.getRa().asDegrees(), corner.getDec().asDegrees()) for corner in tractCorners]
-        minRa = np.min([corner[0] for corner in corners])
-        maxRa = np.max([corner[0] for corner in corners])
-        # If the tract needs wrapping in ra, wrap it
-        if maxRa - minRa > 10:
-            x = maxRa
-            maxRa = 360 + minRa
-            minRa = x
-            minDec = np.min([corner[1] for corner in corners])
-            maxDec = np.max([corner[1] for corner in corners])
-            corners = [(minRa, minDec), (maxRa, minDec), (maxRa, maxDec), (minRa, maxDec)]
-
-        return corners

python/lsst/analysis/tools/utils.py

@@ -0,0 +1,111 @@
+# This file is part of analysis_tools.
+#
+# 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/>.
+
+__all__ = (
+    "getTractCorners",
+    "getPatchCorners",
+)
+
+import numpy as np
+from lsst.geom import Box2D
+
+
+def getTractCorners(skymap, tractId):
+    """Calculate the corners of a tract, given  skymap.
+
+     Parameters
+    ----------
+    skymap : `lsst.skymap`
+    tractId : `int`
+        Identification number of the tract whose corner coordinates
+        are returned.
+
+     Returns
+    -------
+    corners : `list` of `tuples` of `float`
+
+     Notes
+    -----
+    Corners are returned in degrees and wrapped in ra.
+    """
+    tractCorners = skymap[tractId].getVertexList()
+    corners = _wrapRa([(corner.getRa().asDegrees(), corner.getDec().asDegrees()) for corner in tractCorners])
+
+    return corners
+
+
+def getPatchCorners(tractInfo, patchId):
+    """Calculate the corners of a patch, given tractInfo.
+
+     Parameters
+    ----------
+    tractInfo : `lsst.skymap.tractInfo.ExplicitTractInfo`
+        Tract info object of the tract containing the patch whose
+        corner coordinates are returned.
+    patchId : `int`
+        Identification number of the patch whose corner coordinates
+        are returned.
+
+     Returns
+    -------
+    corners : `list` of `tuples` of `float`
+
+     Notes
+    -----
+    Corners are returned in degrees and are wrapped in ra.
+    """
+    patchInfo = tractInfo.getPatchInfo(patchId)
+    patchCorners = Box2D(patchInfo.getInnerBBox()).getCorners()
+
+    tractWcs = tractInfo.getWcs()
+    patchCorners = tractWcs.pixelToSky(patchCorners)
+    corners = _wrapRa([(corner.getRa().asDegrees(), corner.getDec().asDegrees()) for corner in patchCorners])
+
+    return corners
+
+
+def _wrapRa(corners):
+    """Wrap in right ascension if the corners span RA=0
+
+    Parameters
+    ----------
+    corners :  `list` of `tuples` of `float`
+        Pairs of coordinates representing tract or patch corners.
+
+    Returns
+    -------
+    corners : `list` of `tuples` of `float`
+        Pairs of coordinates representing tract or patch corners,
+        wrapped in RA.
+    """
+
+    minRa = np.min([corner[0] for corner in corners])
+    maxRa = np.max([corner[0] for corner in corners])
+    # If the tract needs wrapping in ra, wrap it
+    if maxRa - minRa > 10:
+        x = maxRa
+        maxRa = 360 + minRa
+        minRa = x
+        minDec = np.min([corner[1] for corner in corners])
+        maxDec = np.max([corner[1] for corner in corners])
+        corners = [(minRa, minDec), (maxRa, minDec), (maxRa, maxDec), (minRa, maxDec)]
+
+    return corners

tests/test_utils.py

@@ -0,0 +1,91 @@
+#
+# 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/>.
+
+import unittest
+
+import lsst.skymap as skyMap
+import lsst.utils.tests
+from lsst.analysis.tools.utils import getPatchCorners, getTractCorners
+
+
+class TestTractPatchUtils(lsst.utils.tests.TestCase):
+    """Test to see if the tract and patch corner calculations are working.
+    Includes a test case in which the tract spans the RA=360 line to ensure
+    that RA wrapping is working as expected.
+    """
+
+    def setUp(self):
+
+        # Tract 0, Patch 8 of the following SkyMap spans RA=360.
+        self.tractIds = [0, 1]
+        self.patchIds = [8, 0]
+        skyMapConfig = skyMap.discreteSkyMap.DiscreteSkyMapConfig()
+        skyMapConfig.raList = [0, 180]
+        skyMapConfig.decList = [-1, 1]
+        skyMapConfig.radiusList = [0.1, 0.1]
+        self.skyMap = skyMap.DiscreteSkyMap(skyMapConfig)
+
+        self.tractCorners = [
+            [
+                (358.8900906668926, -2.1096270745156804),
+                (361.10981685029367, -2.1096270745156804),
+                (361.10981685029367, 0.11009493220816949),
+                (358.8900906668926, 0.11009493220816949),
+            ],
+            [
+                (181.10981686420706, -0.11000243565052142),
+                (178.89009065297807, -0.110002449565178),
+                (178.88933977564471, 2.109719508502047),
+                (181.1105676789934, 2.1097195571483582),
+            ],
+        ]
+        self.patchCorners = [
+            [
+                (359.9999537372586, -1.7399434656613333),
+                (360.37005438077574, -1.7399434656613333),
+                (360.37005438077574, -1.369927758180603),
+                (359.9999537372586, -1.369927758180603),
+            ],
+            [
+                (181.10981695667596, -0.11000252814707863),
+                (180.73987541079327, -0.1099561458277576),
+                (180.73992023207617, 0.2600039998656897),
+                (181.10988418299306, 0.25993833301062175),
+            ],
+        ]
+
+    def testTractCorners(self):
+
+        for i, tractId in enumerate(self.tractIds):
+            self.assertEqual(
+                getTractCorners(self.skyMap, tractId),
+                self.tractCorners[i],
+            )
+
+    def testPatchCorners(self):
+
+        for i, (tractId, patchId) in enumerate(zip(self.tractIds, self.patchIds)):
+            tractInfo = self.skyMap.generateTract(tractId)
+            self.assertEqual(getPatchCorners(tractInfo, patchId), self.patchCorners[i])
+
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()