From 7dac39bee744ff7e33fc4e3dee05c43f86130fff Mon Sep 17 00:00:00 2001 From: Tim Jenness Date: Wed, 8 May 2024 16:08:08 -0700 Subject: [PATCH] Add support for parsing IVOA POS string --- python/lsst/sphgeom/__init__.py | 128 ++++++++++++++++++++++++++++++++ tests/test_ivoa.py | 66 ++++++++++++++++ 2 files changed, 194 insertions(+) create mode 100644 tests/test_ivoa.py diff --git a/python/lsst/sphgeom/__init__.py b/python/lsst/sphgeom/__init__.py index 7ffa077..949f3f1 100644 --- a/python/lsst/sphgeom/__init__.py +++ b/python/lsst/sphgeom/__init__.py @@ -28,6 +28,7 @@ """lsst.sphgeom """ +import typing from ._healpixPixelization import * from ._sphgeom import * @@ -37,3 +38,130 @@ from .version import * PixelizationABC.register(Pixelization) + +# Copy and paste from lsst.utils.wrappers +_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): + import sys + + 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_lat(lat: float) -> float: + """Map +Inf to +90 and -Inf to -90 degrees.""" + import math + + if not math.isinf(lat): + return lat + if lat > 0.0: + return 90.0 + return -90.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 `` + * ``RANGE `` + * ``POLYGON ... (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) + elif shape == "RANGE": + import math + + 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. + # If either of the longitude values are infinite, all longitudes + # should be included. + if math.isinf(coordinates[0]) or math.isinf(coordinates[1]): + return Box( + NormalizedAngleInterval.fromDegrees(0.0, 360.0), + AngleInterval.fromDegrees(_inf_to_lat(coordinates[2]), _inf_to_lat(coordinates[3])), + ) + else: + return Box( + LonLat.fromDegrees(coordinates[0], _inf_to_lat(coordinates[2])), + LonLat.fromDegrees(coordinates[1], _inf_to_lat(coordinates[3])), + ) + elif 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}'") + + +del typing diff --git a/tests/test_ivoa.py b/tests/test_ivoa.py new file mode 100644 index 0000000..4634b70 --- /dev/null +++ b/tests/test_ivoa.py @@ -0,0 +1,66 @@ +# 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 . + +import unittest + +from lsst.sphgeom import Region + + +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", + ) + for pos in bad_pos: + with self.assertRaises(ValueError): + Region.from_ivoa_pos(pos) + + +if __name__ == "__main__": + unittest.main()