DM-44261: Add support for parsing IVOA POS string

.pre-commit-config.yaml

@@ -1,6 +1,6 @@
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.4.0
+    rev: v4.6.0
     hooks:
       - id: check-yaml
         args:
@@ -9,20 +9,25 @@ repos:
       - id: trailing-whitespace
       - id: check-toml
   - repo: https://github.com/psf/black
-    rev: 23.1.0
+    rev: 24.4.2
     hooks:
       - id: black
         # It is recommended to specify the latest version of Python
         # supported by your project here, or alternatively use
         # pre-commit's default_language_version, see
         # https://pre-commit.com/#top_level-default_language_version
-        language_version: python3.10
+        language_version: python3.11
   - repo: https://github.com/pycqa/isort
-    rev: 5.12.0
+    rev: 5.13.2
     hooks:
       - id: isort
         name: isort (python)
   - repo: https://github.com/PyCQA/flake8
-    rev: 6.0.0
+    rev: 7.0.0
     hooks:
       - id: flake8
+  - repo: https://github.com/astral-sh/ruff-pre-commit
+    # Ruff version.
+    rev: v0.4.3
+    hooks:
+      - id: ruff

python/lsst/sphgeom/__init__.py

@@ -29,6 +29,7 @@
 """lsst.sphgeom
 """
 
+from ._continue_class import *
 from ._healpixPixelization import *
 from ._sphgeom import *
 from ._sphgeom import Pixelization

python/lsst/sphgeom/_continue_class.py

@@ -0,0 +1,217 @@
+# This file is part of sphgeom.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (http://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This software is dual licensed under the GNU General Public License and also
+# under a 3-clause BSD license. Recipients may choose which of these licenses
+# to use; please see the files gpl-3.0.txt and/or bsd_license.txt,
+# respectively.  If you choose the GPL option then the following text applies
+# (but note that there is still no warranty even if you opt for BSD instead):
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+
+"""Extend any of the C++ Python classes by adding additional methods."""
+
+# Nothing to export.
+__all__ = []
+
+import math
+import sys
+import typing
+
+from ._sphgeom import Angle, Box, Circle, ConvexPolygon, LonLat, Region, UnitVector3d
+
+# Copy and paste from lsst.utils.wrappers:
+# * INTRINSIC_SPECIAL_ATTRIBUTES
+# * isAttributeSafeToTransfer
+# * continueClass
+_INTRINSIC_SPECIAL_ATTRIBUTES = frozenset(
+    (
+        "__qualname__",
+        "__module__",
+        "__metaclass__",
+        "__dict__",
+        "__weakref__",
+        "__class__",
+        "__subclasshook__",
+        "__name__",
+        "__doc__",
+    )
+)
+
+
+def _isAttributeSafeToTransfer(name: str, value: typing.Any) -> bool:
+    if name.startswith("__") and (
+        value is getattr(object, name, None) or name in _INTRINSIC_SPECIAL_ATTRIBUTES
+    ):
+        return False
+    return True
+
+
+def _continueClass(cls):
+    orig = getattr(sys.modules[cls.__module__], cls.__name__)
+    for name in dir(cls):
+        # Common descriptors like classmethod and staticmethod can only be
+        # accessed without invoking their magic if we use __dict__; if we use
+        # getattr on those we'll get e.g. a bound method instance on the dummy
+        # class rather than a classmethod instance we can put on the target
+        # class.
+        attr = cls.__dict__.get(name, None) or getattr(cls, name)
+        if _isAttributeSafeToTransfer(name, attr):
+            setattr(orig, name, attr)
+    return orig
+
+
+def _inf_to_limit(value: float, min: float, max: float) -> float:
+    """Map a value to a fixed range if infinite."""
+    if not math.isinf(value):
+        return value
+    if value > 0.0:
+        return max
+    return min
+
+
+def _inf_to_lat(lat: float) -> float:
+    """Map latitude +Inf to +90 and -Inf to -90 degrees."""
+    return _inf_to_limit(lat, -90.0, 90.0)
+
+
+def _inf_to_lon(lat: float) -> float:
+    """Map longitude +Inf to +360 and -Inf to 0 degrees."""
+    return _inf_to_limit(lat, 0.0, 360.0)
+
+
+@_continueClass
+class Region:
+    """A minimal interface for 2-dimensional regions on the unit sphere."""
+
+    @classmethod
+    def from_ivoa_pos(cls, pos: str) -> Region:
+        """Create a Region from an IVOA POS string.
+
+        Parameters
+        ----------
+        pos : `str`
+            A string using the IVOA SIAv2 POS syntax.
+
+        Returns
+        -------
+        region : `Region`
+            A region equivalent to the POS string.
+
+        Notes
+        -----
+        See
+        https://ivoa.net/documents/SIA/20151223/REC-SIA-2.0-20151223.html#toc12
+        for a description of the POS parameter but in summary the options are:
+
+        * ``CIRCLE <longitude> <latitude> <radius>``
+        * ``RANGE <longitude1> <longitude2> <latitude1> <latitude2>``
+        * ``POLYGON <longitude1> <latitude1> ... (at least 3 pairs)``
+
+        Units are degrees in all coordinates.
+        """
+        shape, *coordinates = pos.split()
+        coordinates = tuple(float(c) for c in coordinates)
+        n_floats = len(coordinates)
+        if shape == "CIRCLE":
+            if n_floats != 3:
+                raise ValueError(f"CIRCLE requires 3 numbers but got {n_floats} in '{pos}'.")
+            center = LonLat.fromDegrees(coordinates[0], coordinates[1])
+            radius = Angle.fromDegrees(coordinates[2])
+            return Circle(UnitVector3d(center), radius)
+
+        if shape == "RANGE":
+            if n_floats != 4:
+                raise ValueError(f"RANGE requires 4 numbers but got {n_floats} in '{pos}'.")
+            # POS allows +Inf and -Inf in ranges. These are not allowed by
+            # Box and so must be converted.
+            return Box(
+                LonLat.fromDegrees(_inf_to_lon(coordinates[0]), _inf_to_lat(coordinates[2])),
+                LonLat.fromDegrees(_inf_to_lon(coordinates[1]), _inf_to_lat(coordinates[3])),
+            )
+
+        if shape == "POLYGON":
+            if n_floats % 2 != 0:
+                raise ValueError(f"POLYGON requires even number of floats but got {n_floats} in '{pos}'.")
+            if n_floats < 6:
+                raise ValueError(
+                    f"POLYGON specification requires at least 3 coordinates, got {n_floats // 2} in '{pos}'"
+                )
+            # Coordinates are x1, y1, x2, y2, x3, y3...
+            # Get pairs by skipping every other value.
+            pairs = list(zip(coordinates[0::2], coordinates[1::2]))
+            vertices = [LonLat.fromDegrees(lon, lat) for lon, lat in pairs]
+            return ConvexPolygon([UnitVector3d(c) for c in vertices])
+
+        raise ValueError(f"Unrecognized shape in POS string '{pos}'")
+
+    def to_ivoa_pos(self) -> str:
+        """Represent the region as an IVOA POS string.
+
+        Returns
+        -------
+        pos : `str`
+            The region in ``POS`` format.
+        """
+        raise NotImplementedError("This region can not be converted to an IVOA POS string.")
+
+
+@_continueClass
+class Circle:  # noqa: F811
+    """A circular region on the unit sphere that contains its boundary."""
+
+    def to_ivoa_pos(self) -> str:
+        # Docstring inherited.
+        center = LonLat(self.getCenter())
+        lon = center.getLon().asDegrees()
+        lat = center.getLat().asDegrees()
+        rad = self.getOpeningAngle().asDegrees()
+        return f"CIRCLE {lon} {lat} {rad}"
+
+
+@_continueClass
+class Box:  # noqa: F811
+    """A rectangle in spherical coordinate space that contains its boundary."""
+
+    def to_ivoa_pos(self) -> str:
+        # Docstring inherited.
+        lon_range = self.getLon()
+        lat_range = self.getLat()
+
+        lon1 = lon_range.getA().asDegrees()
+        lon2 = lon_range.getB().asDegrees()
+        lat1 = lat_range.getA().asDegrees()
+        lat2 = lat_range.getB().asDegrees()
+
+        # Do not attempt to map to +/- Inf -- there is no way to know if
+        # that is any better than 0. -> 360.
+        return f"RANGE {lon1} {lon2} {lat1} {lat2}"
+
+
+@_continueClass
+class ConvexPolygon:  # noqa: F811
+    """A rectangle in spherical coordinate space that contains its boundary."""
+
+    def to_ivoa_pos(self) -> str:
+        # Docstring inherited.
+        coords = (LonLat(v) for v in self.getVertices())
+        coord_strings = [f"{c.getLon().asDegrees()} {c.getLat().asDegrees()}" for c in coords]
+
+        return f"POLYGON {' '.join(coord_strings)}"

tests/test_ivoa.py

@@ -0,0 +1,132 @@
+# This file is part of sphgeom.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (http://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This software is dual licensed under the GNU General Public License and also
+# under a 3-clause BSD license. Recipients may choose which of these licenses
+# to use; please see the files gpl-3.0.txt and/or bsd_license.txt,
+# respectively.  If you choose the GPL option then the following text applies
+# (but note that there is still no warranty even if you opt for BSD instead):
+#
+# 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 <http://www.gnu.org/licenses/>.
+
+import unittest
+
+from lsst.sphgeom import Box, Circle, ConvexPolygon, LonLat, Region, UnitVector3d
+
+
+class IvoaTestCase(unittest.TestCase):
+    """Test Box."""
+
+    def test_construction(self):
+        # Example POS strings found in IVOA documentation.
+        example_pos = (
+            "CIRCLE 12.0 34.0 0.5",
+            "RANGE 12.0 12.5 34.0 36.0",
+            "POLYGON 12.0 34.0 14.0 35.0 14. 36.0 12.0 35.0",
+            "RANGE 0 360.0 -2.0 2.0",
+            "RANGE 0 360.0 89.0 +Inf",
+            "RANGE -Inf +Inf -Inf +Inf",
+            "POLYGON 12 34 14 34 14 36 12 36",
+            "RANGE 0 360 89 90",
+        )
+        for pos in example_pos:
+            region = Region.from_ivoa_pos(pos)
+            self.assertIsInstance(region, Region)
+
+        # Badly formed strings raising ValueError.
+        bad_pos = (
+            "circle 12 34 0.5",
+            "CIRCLE 12 34 1 1",
+            "RANGE 0 360",
+            "POLYGON 0 1 2 3",
+            "POLYGON 0 1 2 3 4 5 6",
+            "CONVEXPOLYGON 0 1 2 3 4 5",
+        )
+        for pos in bad_pos:
+            with self.assertRaises(ValueError):
+                Region.from_ivoa_pos(pos)
+
+    def _split_pos(self, pos: str):
+        """Split POS into type and floats."""
+        region_type, *coordstr = pos.split()
+        coordinates = [float(c) for c in coordstr]
+        return region_type, coordinates
+
+    def assert_pos_equal(self, pos1: str, pos2: str):
+        """Compare two POS strings and check for equality taking into account
+        floating point differences.
+        """
+        region_type1, coords1 = self._split_pos(pos1)
+        region_type2, coords2 = self._split_pos(pos2)
+        self.assertEqual(region_type1, region_type2)
+        self.assertEqual(len(coords1), len(coords2))
+        for c1, c2 in zip(coords1, coords2):
+            self.assertAlmostEqual(c1, c2)
+
+    def test_circle(self):
+        """Test circle construction."""
+        pos = "CIRCLE 12.0 34.0 5"
+        circle = Region.from_ivoa_pos(pos)
+        self.assertIsInstance(circle, Circle)
+        self.assertTrue(circle.contains(UnitVector3d(LonLat.fromDegrees(13.0, 33.0))))
+        self.assertFalse(circle.contains(UnitVector3d(LonLat.fromDegrees(12.0, 40.0))))
+
+        self.assert_pos_equal(circle.to_ivoa_pos(), pos)
+
+    def test_range(self):
+        """Test range construction."""
+        box = Region.from_ivoa_pos("RANGE 1 2 5 6")
+        self.assertIsInstance(box, Box)
+        self.assertTrue(box.contains(UnitVector3d(LonLat.fromDegrees(1.5, 5.4))))
+        self.assertFalse(box.contains(UnitVector3d(LonLat.fromDegrees(4, 10))))
+        self.assert_pos_equal(box.to_ivoa_pos(), "RANGE 1 2 5 6")
+
+        box = Region.from_ivoa_pos("RANGE 1 2 20 +Inf")
+        self.assertTrue(box.contains(UnitVector3d(LonLat.fromDegrees(1.7, 80))))
+        self.assertFalse(box.contains(UnitVector3d(LonLat.fromDegrees(1.7, 10))))
+        self.assert_pos_equal(box.to_ivoa_pos(), "RANGE 1 2 20 90")
+
+        box = Region.from_ivoa_pos("RANGE 50 +Inf 20 30")
+        self.assertTrue(box.contains(UnitVector3d(LonLat.fromDegrees(60, 25))))
+        self.assertFalse(box.contains(UnitVector3d(LonLat.fromDegrees(49, 21))))
+        self.assert_pos_equal(box.to_ivoa_pos(), "RANGE 50 360 20 30")
+
+        box = Region.from_ivoa_pos("RANGE -Inf +50 20 30")
+        self.assertTrue(box.contains(UnitVector3d(LonLat.fromDegrees(40, 25))))
+        self.assertFalse(box.contains(UnitVector3d(LonLat.fromDegrees(60, 21))))
+        self.assert_pos_equal(box.to_ivoa_pos(), "RANGE 0 50 20 30")
+
+        box = Region.from_ivoa_pos("RANGE -Inf +Inf 20 30")
+        self.assertTrue(box.contains(UnitVector3d(LonLat.fromDegrees(10, 25))))
+        self.assertTrue(box.contains(UnitVector3d(LonLat.fromDegrees(359, 25))))
+        self.assertFalse(box.contains(UnitVector3d(LonLat.fromDegrees(49, 19))))
+        self.assert_pos_equal(box.to_ivoa_pos(), "RANGE 0 360. 20 30")
+
+    def test_polygon(self):
+        """Test polygon construction."""
+        pos = "POLYGON 12.0 34.0 14.0 35.0 14. 36.0 12.0 35.0"
+        poly = Region.from_ivoa_pos(pos)
+        self.assertIsInstance(poly, ConvexPolygon)
+        self.assertTrue(poly.contains(UnitVector3d(LonLat.fromDegrees(13, 35))))
+        self.assertFalse(poly.contains(UnitVector3d(LonLat.fromDegrees(14, 34))))
+        self.assert_pos_equal(poly.to_ivoa_pos(), pos)
+
+
+if __name__ == "__main__":
+    unittest.main()