diff --git a/regions/shapes/__init__.py b/regions/shapes/__init__.py index 28f24a4e..ca9c8c30 100644 --- a/regions/shapes/__init__.py +++ b/regions/shapes/__init__.py @@ -6,6 +6,7 @@ from .point import * from .polygon import * from .rectangle import * +from .range import * from .annulus import * from .line import * from .text import * diff --git a/regions/shapes/circle.py b/regions/shapes/circle.py index 2c9ded75..4e6c8967 100644 --- a/regions/shapes/circle.py +++ b/regions/shapes/circle.py @@ -11,7 +11,7 @@ from ..core.attributes import (ScalarSky, ScalarPix, QuantityLength, ScalarLength, RegionVisual, RegionMeta) -__all__ = ['CirclePixelRegion', 'CircleSkyRegion'] +__all__ = ['CirclePixelRegion', 'CircleSkyRegion', 'CircleSphericalRegion'] class CirclePixelRegion(PixelRegion): @@ -195,3 +195,15 @@ def to_pixel(self, wcs): center, scale, _ = skycoord_to_pixel_scale_angle(self.center, wcs) radius = self.radius.to('deg').value * scale return CirclePixelRegion(center, radius, self.meta, self.visual) + + +class CircleSphericalRegion(CircleSkyRegion): + """Spherical circle sky region. + + TODO: document + """ + + def contains(self, skycoord, wcs=None): + """Defined by spherical distance.""" + separation = self.center.separation(skycoord) + return separation < self.radius diff --git a/regions/shapes/range.py b/regions/shapes/range.py new file mode 100644 index 00000000..f64cd8eb --- /dev/null +++ b/regions/shapes/range.py @@ -0,0 +1,72 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from astropy.coordinates import Angle, SkyCoord +from regions.core.attributes import QuantityLength +from ..core import SkyRegion, PixelRegion, PixCoord + +__all__ = ["RangePixelRegion", "RangeSphericalRegion"] + +inf = Angle(float("inf"), "deg") + + +# TODO: needed? +class RangePixelRegion(PixelRegion): + """ + A range in pixel coordinates. + + See `RangeSphericalRegion`. + """ + + +class RangeSphericalRegion(SkyRegion): + """ + A range in longitude and latitude. + + Use case: http://www.ivoa.net/documents/SIA/20151223/REC-SIA-2.0-20151223.html#toc12 + + Parameters + ---------- + lon_min, lon_max : `~astropy.coordinates.Angle` + Range in longitude + lat_min, lat_max : `~astropy.coordinates.Angle` + Range in latitude + frame : `~astropy.coordinates.BaseCoordinateFrame` + Sky frame + meta : `~regions.RegionMeta` object, optional + A dictionary which stores the meta attributes of this region. + visual : `~regions.RegionVisual` object, optional + A dictionary which stores the visual meta attributes of this region. + """ + # TODO: should we use these properties, or SkyCoord for min / max, + # or a center-based representation? + _params = ('lon_min', 'lon_max', 'lat_min', 'lat_max') + lon_min = QuantityLength('lon_min') + lon_max = QuantityLength('lon_max') + lat_min = QuantityLength('lat_min') + lat_max = QuantityLength('lat_max') + + def __init__(self, lon_min=-inf, lon_max=inf, lat_min=-inf, lat_max=inf, frame="icrs", meta=None, visual=None): + self.lon_min = lon_min + self.lon_max = lon_max + self.lat_min = lat_min + self.lat_max = lat_max + # TODO: Is there a helper function to make a frame without a + self.frame = SkyCoord(0, 0, unit="deg", frame=frame).frame + self.meta = meta or {} + self.visual = visual or {} + + def to_pixel(self, wcs): + sky_min = SkyCoord(self.lon_min, self.lat_min, frame=self.frame) + sky_max = SkyCoord(self.lon_max, self.lat_max, frame=self.frame) + + pix_min = PixCoord.from_sky(sky_min, self.wcs) + pix_max = PixCoord.from_sky(sky_max, self.wcs) + return RangePixelRegion(pix_min.x, pix_max.x, pix_min.y, pix_max.y, self.meta, self.visual) + + def contains(self, skycoord, wcs=None): + skycoord = skycoord.transform_to(self.frame) + return ( + (self.lon_min <= skycoord.data.lon) & + (self.lon_max >= skycoord.data.lon) & + (self.lat_min <= skycoord.data.lat) & + (self.lat_max >= skycoord.data.lat) + ) diff --git a/regions/shapes/tests/test_circle.py b/regions/shapes/tests/test_circle.py index 8505cd20..dec69334 100644 --- a/regions/shapes/tests/test_circle.py +++ b/regions/shapes/tests/test_circle.py @@ -12,7 +12,7 @@ from ...core import PixCoord from ...tests.helpers import make_simple_wcs -from ..circle import CirclePixelRegion, CircleSkyRegion +from ..circle import CirclePixelRegion, CircleSkyRegion, CircleSphericalRegion from .utils import HAS_MATPLOTLIB # noqa from .test_common import BaseTestPixelRegion, BaseTestSkyRegion @@ -101,3 +101,14 @@ def test_contains(self, wcs): position = SkyCoord([1, 3] * u.deg, [2, 4] * u.deg) # 1,2 is outside, 3,4 is the center and is inside assert all(self.reg.contains(position, wcs) == np.array([False, True], dtype='bool')) + +class TestCircleSphericalRegion: + @staticmethod + def test_contains(wcs): + # TODO: pick a wcs and point where sky and spherical differs + reg_sky = CircleSkyRegion(SkyCoord(0 * u.deg, 80 * u.deg), radius=12 * u.deg) + reg_sph = CircleSphericalRegion(SkyCoord(0 * u.deg, 80 * u.deg), radius=12 * u.deg) + point = SkyCoord(180 * u.deg, 89 * u.deg) + assert reg_sky.contains(point, wcs) + assert reg_sph.contains(point, wcs) + assert reg_sph.contains(point) diff --git a/regions/shapes/tests/test_range.py b/regions/shapes/tests/test_range.py new file mode 100644 index 00000000..3dfc8290 --- /dev/null +++ b/regions/shapes/tests/test_range.py @@ -0,0 +1,23 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from numpy.testing import assert_equal +import astropy.units as u +from astropy.coordinates import SkyCoord +from ..range import RangeSphericalRegion + +# TODO: use test cases with finite & infinite values from here: +# http://www.ivoa.net/documents/SIA/20151223/REC-SIA-2.0-20151223.html#toc12 + +def test_range_default(): + region = RangeSphericalRegion() + assert region.contains(SkyCoord("0d 0d")) + +def test_range_finite(): + region = RangeSphericalRegion(12 * u.deg, 12.5 * u.deg, 34 * u.deg, 36.0 * u.deg) + assert region.contains(SkyCoord("12d 34d")) + assert not region.contains(SkyCoord("11d 34d")) + assert not region.contains(SkyCoord("13d 34d")) + assert not region.contains(SkyCoord("12d 33d")) + assert not region.contains(SkyCoord("12d 37d")) + + mask = region.contains(SkyCoord([12, 13] * u.deg, [34, 34] * u.deg)) + assert_equal(mask, [True, False])