diff --git a/CHANGES.rst b/CHANGES.rst index 174c6661..0eda5a52 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -7,6 +7,27 @@ General New Features ------------ +- Regions now includes ``SphericalSkyRegion`` ("region-on-celestial-sphere"), + complementing the implicitly planar ``SkyRegion`` ("region-on-images"). + ``SphericalSkyRegion`` does not require a WCS to determine whether points + are contained within the region or not (unlike ``SkyRegion``). + Additionally, ``SphericalSkyRegion`` classes include the method + ``transform_to`` to transform the regions between different + celestial coordinate reference frames. + It is also possible to transform between spherical and planar + (sky or pixel) regions (with ``to_sky``, ``to_pixel``, and ``to_spherical_sky``, + as appropriate), with the option to capture boundary distortions due to + projection effects between spherical and planar geometries + or to ignore them (e.g., assuming a circle stays a perfect circle). + Current spherical shapes supported include: ``CircleSphericalSkyRegion``, + ``CircleAnnulusSphericalSkyRegion``, ``RangeSphericalSkyRegion`` + (i.e., bounded by ranges of longitude and/or latitude), and + ``LuneSphericalSkyRegion`` (a slice between two great circles, + such as between two lines of longitude). + Support for additional spherical shapes, and for all cases of + planar <-> spherical transformation (where well defined) may be added + at a future date. [#618] + Bug Fixes --------- diff --git a/docs/compound.rst b/docs/compound.rst index 041bcfb8..e8e1d6dc 100644 --- a/docs/compound.rst +++ b/docs/compound.rst @@ -2,8 +2,8 @@ Combining Regions ================= There are a few ways to combine any two `~regions.Region` objects -into a compound region, i.e., a `~regions.CompoundPixelRegion` or -`~regions.CompoundSkyRegion` object. +into a compound region, i.e., a `~regions.CompoundPixelRegion`, +`~regions.CompoundSkyRegion`, or `~regions.CompoundSphericalSkyRegion` object. Let's start by defining two sky regions:: @@ -120,9 +120,14 @@ operator or the :meth:`~regions.Region.symmetric_difference` method:: Example Illustrating Compound Regions ------------------------------------- +The following examples demonstrate how to combine planar sky regions +and spherical sky regions, with the same circle centers and radii. + .. plot:: :include-source: false + # planar + import matplotlib.pyplot as plt import numpy as np from astropy.coordinates import Angle, SkyCoord @@ -141,6 +146,7 @@ Example Illustrating Compound Regions # remove sources dataset.image.data = np.zeros_like(dataset.image.data) + #---------------------------------------- # define 2 sky circles circle1 = CircleSkyRegion( center=SkyCoord(20, 0, unit='deg', frame='galactic'), @@ -156,6 +162,7 @@ Example Illustrating Compound Regions coords = np.array(np.meshgrid(lon, lat)).T.reshape(-1, 2) skycoords = SkyCoord(coords, unit='deg', frame='galactic') + #---------------------------------------- # get events in AND and XOR compound_and = circle1 & circle2 compound_xor = circle1 ^ circle2 @@ -167,6 +174,7 @@ Example Illustrating Compound Regions # plot fig = plt.figure() + fig.set_size_inches(7,3.5) ax = fig.add_axes([0.15, 0.1, 0.8, 0.8], projection=wcs, aspect='equal') ax.scatter(skycoords.l.value, skycoords.b.value, label='all', @@ -185,3 +193,80 @@ Example Illustrating Compound Regions ax.set_xlim(-0.5, dataset.config['shape'][1] - 0.5) ax.set_ylim(-0.5, dataset.config['shape'][0] - 0.5) + ax.set_title("Planar SkyRegions") + + +.. plot:: + :include-source: false + + # spherical + + import matplotlib.pyplot as plt + import numpy as np + from astropy.coordinates import Angle, SkyCoord + + from regions import CircleSphericalSkyRegion, make_example_dataset + + # load example dataset to get skymap + config = dict(crval=(0, 0), + crpix=(180, 90), + cdelt=(-1, 1), + shape=(180, 360)) + + dataset = make_example_dataset(data='simulated', config=config) + wcs = dataset.wcs + + # remove sources + dataset.image.data = np.zeros_like(dataset.image.data) + + #---------------------------------------- + # define 2 spherical sky circles + sph_circle1 = CircleSphericalSkyRegion( + center=SkyCoord(20, 0, unit='deg', frame='galactic'), + radius=Angle('30 deg')) + + sph_circle2 = CircleSphericalSkyRegion( + center=SkyCoord(50, 45, unit='deg', frame='galactic'), + radius=Angle('30 deg')) + + # define skycoords + lon = np.arange(-180, 181, 10) + lat = np.arange(-90, 91, 10) + coords = np.array(np.meshgrid(lon, lat)).T.reshape(-1, 2) + skycoords = SkyCoord(coords, unit='deg', frame='galactic') + + #---------------------------------------- + # get events in AND and XOR + sph_compound_and = sph_circle1 & sph_circle2 + sph_compound_xor = sph_circle1 ^ sph_circle2 + + sph_mask_and = sph_compound_and.contains(skycoords) + sph_skycoords_and = skycoords[sph_mask_and] + sph_mask_xor = sph_compound_xor.contains(skycoords) + sph_skycoords_xor = skycoords[sph_mask_xor] + + # plot + fig = plt.figure() + fig.set_size_inches(7,3.5) + ax = fig.add_axes([0.15, 0.1, 0.8, 0.8], projection=wcs, aspect='equal') + + ax.scatter(skycoords.l.value, skycoords.b.value, label='all', + transform=ax.get_transform('galactic')) + ax.scatter(sph_skycoords_xor.l.value, sph_skycoords_xor.b.value, color='orange', + label='xor', transform=ax.get_transform('galactic')) + ax.scatter(sph_skycoords_and.l.value, sph_skycoords_and.b.value, color='magenta', + label='and', transform=ax.get_transform('galactic')) + + boundary_kwargs = dict( + include_boundary_distortions=True, discretize_kwargs={"n_points":1000} + ) + sph_circle1.to_pixel(wcs=wcs,**boundary_kwargs).plot(ax=ax, edgecolor='green', facecolor='none', + alpha=0.8, lw=3) + sph_circle2.to_pixel(wcs=wcs,**boundary_kwargs).plot(ax=ax, edgecolor='red', facecolor='none', + alpha=0.8, lw=3) + + ax.legend(loc='lower right') + + ax.set_xlim(-0.5, dataset.config['shape'][1] - 0.5) + ax.set_ylim(-0.5, dataset.config['shape'][0] - 0.5) + ax.set_title("Spherical SkyRegions") diff --git a/docs/contains.rst b/docs/contains.rst index 9b50e2e5..784ac112 100644 --- a/docs/contains.rst +++ b/docs/contains.rst @@ -1,7 +1,10 @@ Checking for Points Inside Regions ================================== -Let's start by defining both a sky and pixel region:: +Points Inside Planar Regions +---------------------------- + +Let's start by defining both a planar sky and pixel region:: >>> from astropy.coordinates import Angle, SkyCoord >>> from regions import CircleSkyRegion, PixCoord, CirclePixelRegion @@ -60,3 +63,34 @@ Note that `regions.SkyRegion.contains` requires a WCS to be passed:: >>> skycoord = SkyCoord([50, 50], [10, 60], unit='deg') >>> sky_region.contains(skycoord, wcs) array([False, True]) + + +Points Inside Spherical Regions +------------------------------- + +For `~regions.SphericalSkyRegion` objects, checking whether point(s) are +contained inside that region requires no other input --- since these +regions are defined with a the spherical geometry, and not a projected geometry +(as captured through the projection encoded in a WCS) as in +`~regions.SkyRegion`. + +Let's define a spherical sky region:: + + >>> from regions import CircleSphericalSkyRegion + + >>> sph_sky_center = SkyCoord(42, 43, unit='deg') + >>> sph_sky_radius = Angle(25, 'deg') + >>> sph_sky_region = CircleSphericalSkyRegion(sph_sky_center, + ... sph_sky_radius) + >>> print(sph_sky_region) + Region: CircleSphericalSkyRegion + center: + radius: 25.0 deg + +Use the `~regions.SphericalSkyRegion.contains()` method to determine which +point(s) lie inside or outside the region:: + + >>> skycoord = SkyCoord([50, 50], [10, 60], unit='deg') + >>> sph_sky_region.contains(skycoord) + array([False, True]) diff --git a/docs/getting_started.rst b/docs/getting_started.rst index c67b08cd..6ea49481 100644 --- a/docs/getting_started.rst +++ b/docs/getting_started.rst @@ -10,14 +10,22 @@ Introduction The Regions package provides classes to represent: -* Regions defined using pixel coordinates (e.g., +* Regions defined using pixel coordinates (a "region-on-image"; e.g., `~regions.CirclePixelRegion`) * Regions defined using celestial coordinates, but still in an Euclidean - geometry (e.g., `~regions.CircleSkyRegion`) + geometry (i.e., a planar projection, as a "region-on-image"; + e.g., `~regions.CircleSkyRegion`) +* Regions defined using celestial coordinates, and with a spherical + geometry (a "region-on-celestial-sphere"; e.g., `~regions.CircleSphericalSkyRegion`) -To transform between sky and pixel regions, a `world coordinate system +To transform between (planar) sky and pixel regions, a `world coordinate system `_ object (e.g., -`astropy.wcs.WCS`) is needed. +`astropy.wcs.WCS`) is needed. To transform between spherical and planar (sky or pixel) +regions, in addition to a `wcs +`_, it is also +necessary to specify whether or not boundary distortions should be included +(capturing the projection effects inherent in projection-to-spherical +transformations, or the inverse). Regions also provides a unified interface for reading, writing, parsing, and serializing regions data in different formats, including @@ -187,10 +195,10 @@ Sky Regions Sky regions are regions that are defined using celestial coordinates. Please note they are **not** defined as regions on the celestial sphere, -but rather are meant to represent shapes on an image. They simply use -sky coordinates instead of pixel coordinates to define their position. -The remaining shape parameters are converted to pixels using the pixel -scale of the image. +but rather are meant to represent shapes on an image ("region-on-image"). +They simply use sky coordinates instead of pixel coordinates to define +their position. The remaining shape parameters are converted to pixels +using the pixel scale of the image. Let's create a sky region: @@ -236,3 +244,67 @@ You can access its properties via attributes: See the :ref:`shapes` documentation for the complete list of pixel-based regions and to learn more about :class:`~regions.Region` objects and their capabilities. + + +Spherical Sky Regions +--------------------- + +Spherical sky regions are defined using celestial coordinates, +and **are** defined as regions on the celestial sphere +("regions-on-celestial-sphere", in contrast to the planar Sky Regions). +In order to transform between spherical and planar ("region-on-image") regions, +the planar projection (encoded in a `world coordinate system +`_ object; e.g., +`astropy.wcs.WCS`) must be specified, along with a specification of +whether or not boundary distortions should be included. +These distortions (implemented through discrete boundary sampling) +capture the impact of the spherical-to-planar (or vice versa) projection. +However, it is possible to ignore these distortions (e.g., +transforming a spherical circle to a planar circle). + +Spherical sky regions are created using celestial coordinates (as +`~astropy.coordinates.SkyCoord`) and angular distances, +for instance specified as + +.. code-block:: python + + >>> from astropy.coordinates import Angle, SkyCoord + >>> from regions import CircleSphericalSkyRegion + >>> center = SkyCoord(42, 43, unit='deg') + >>> radius = Angle(3, 'deg') + >>> region = CircleSphericalSkyRegion(center, radius) + +Alternatively, one can define the radius using a +`~astropy.units.Quantity` object with angular units: + +.. code-block:: python + + >>> import astropy.units as u + >>> from regions import CircleSphericalSkyRegion + >>> center = SkyCoord(42, 43, unit='deg') + >>> radius = 3.0 * u.deg + >>> region = CircleSphericalSkyRegion(center, radius) + +You can print the region to get some information about its properties: + +.. code-block:: python + + >>> print(region) + Region: CircleSphericalSkyRegion + center: + radius: 3.0 deg + +You can access its properties via attributes: + +.. code-block:: python + + >>> region.center + + >>> region.radius + + +See the :ref:`shapes` documentation for the complete list of pixel-based +regions and to learn more about :class:`~regions.Region` objects and +their capabilities. diff --git a/docs/index.rst b/docs/index.rst index 9a488157..cf7cdaa3 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -32,22 +32,43 @@ Getting Started User Documentation ================== +Defining Regions +---------------- + .. toctree:: :maxdepth: 1 shapes metadata + +Using Regions +------------- + +.. toctree:: + :maxdepth: 1 + contains compound masks + spherical_frame_transform + spherical_bounding_regions plotting + +I/O & Reference +--------------------- + +.. toctree:: + :maxdepth: 1 + region_io shapely api -Example -======= +.. _index_examples: + +Examples +======== The following example shows how to read a DS9 region file and plot the regions on an image using Matplotlib. @@ -72,3 +93,63 @@ regions on an image using Matplotlib. regions = Regions.read(region_file, format='ds9') for i, region in enumerate(regions): region.plot(ax=ax) + +The next example demonstrates how to plot spherical regions (e.g., a circle +and a longitude/latitude "range") on a full sky image (seen in the +`visualization demo here`_) using Matplotlib. + +.. plot:: + + from astropy.wcs import WCS + from astropy.io import fits + from astropy.utils.data import get_pkg_data_filename + import astropy.units as u + from astropy.coordinates import SkyCoord + from astropy.visualization.wcsaxes.frame import EllipticalFrame + from matplotlib import patheffects + from matplotlib import pyplot as plt + + from regions import CircleSphericalSkyRegion, RangeSphericalSkyRegion + + image_file = get_pkg_data_filename('allsky/allsky_rosat.fits') + image_data = fits.getdata(image_file, ext=0, memmap=False) + image_hdr = fits.getheader(image_file) + wcs = WCS(image_hdr) + + + fig, ax = plt.subplots(figsize=(8,8), + subplot_kw=dict(projection=wcs, + frame_class=EllipticalFrame)) + + path_effects=[patheffects.withStroke(linewidth=3, foreground='black')] + ax.coords.grid(color='white') + ax.coords['glon'].set_ticklabel(color='white', path_effects=path_effects) + ax.set_axisbelow(True) + ax.set_xlabel(r"Galactic $\ell$", labelpad=10) + ax.set_ylabel(r"Galactic $b$") + im = ax.imshow(image_data, vmin=0., vmax=300., cmap='gray') + + # Clip the image to the frame + im.set_clip_path(ax.coords.frame.patch) + + # Image is in galactic coordinates + circ = CircleSphericalSkyRegion(SkyCoord(80,45,unit=u.deg,frame="galactic"), + 25*u.deg) + circ.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, edgecolor='red', + facecolor='none', lw=2) + + range = RangeSphericalSkyRegion( + frame="galactic", + longitude_range=[210,330]*u.deg, + latitude_range=[-55,-15]*u.deg, + ) + range.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":250} + ).plot(ax=ax, edgecolor='blue', + facecolor='none', lw=2) diff --git a/docs/plotting.rst b/docs/plotting.rst index 85988325..c9317d7e 100644 --- a/docs/plotting.rst +++ b/docs/plotting.rst @@ -1,6 +1,9 @@ Plotting Regions with Matplotlib ================================ +Plotting Pixel Regions +---------------------- + Some `~regions.PixelRegion` objects have an ``as_artist()`` method that returns an equivalent `matplotlib.patches` object. For example :meth:`regions.CirclePixelRegion.as_artist` returns a @@ -78,3 +81,55 @@ to convert it to a pixel region (using a WCS object): >>> sky_region = CircleSkyRegion(sky_center, sky_radius) >>> pixel_region = sky_region.to_pixel(wcs) >>> pixel_region.plot() + + +Plotting Spherical Sky Regions +------------------------------ + +Similarly, `~regions.SphericalSkyRegion` objects do not have an +``as_artist()`` or ``plot()`` method. To plot a `~regions.SphericalSkyRegion` +object, you will need to convert it to a pixel region (using a WCS object). +Boundary distortions effects can also be included in this conversion (by +setting the ``include_boundary_distortions`` keyword), to capture +the effects of projecting from spherical to a planar geometry. +(See the second example in :ref:`index_examples`.) + +It is also possible to use the coordinates of a discretized +`~regions.SphericalSkyRegion` to show the region's boundary in a figure. + +.. plot:: + :include-source: + + from astropy.coordinates import Angle, SkyCoord + import matplotlib.pyplot as plt + + from regions import CircleSphericalSkyRegion, make_example_dataset + + dataset = make_example_dataset(data='simulated') + wcs = dataset.wcs + + sph_sky_center = SkyCoord(42, 30, unit='deg', frame='galactic') + sph_sky_radius = Angle(12, 'deg') + sph_sky_region = CircleSphericalSkyRegion(sph_sky_center, sph_sky_radius) + + fig, ax = plt.subplots(subplot_kw=dict(projection=wcs)) + ax.grid(True) + ax.set_xlabel(r"Galactic $\ell$") + ax.set_ylabel(r"Galactic $b$") + + sph_sky_region.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, color='tab:red', lw=3) + + sph_sky_center2 = SkyCoord(42, 43, unit='deg', frame='galactic') + sph_sky_radius2 = Angle(25, 'deg') + sph_sky_region2 = CircleSphericalSkyRegion(sph_sky_center2, sph_sky_radius2) + poly_sph_sky2 = sph_sky_region2.discretize_boundary(n_points=1000) + ax.plot( + poly_sph_sky2.vertices.l, + poly_sph_sky2.vertices.b, + lw=2, color="tab:blue", + transform=ax.get_transform('galactic'), + ) diff --git a/docs/references.txt b/docs/references.txt index b60d3108..7a52d4c8 100644 --- a/docs/references.txt +++ b/docs/references.txt @@ -6,3 +6,4 @@ .. _GitHub repository: https://github.com/astropy/regions .. _Region documentation: https://astropy-regions.readthedocs.io/en/latest/ .. _coordinated package: https://www.astropy.org/affiliated/index.html#coordinated-packages +.. _visualization demo here: https://docs.astropy.org/en/stable/visualization/wcsaxes/custom_frames.html diff --git a/docs/region_io.rst b/docs/region_io.rst index 3aced3cb..54b36bbd 100644 --- a/docs/region_io.rst +++ b/docs/region_io.rst @@ -181,3 +181,10 @@ Region File Format Limitations ds9_io fits_io + + +Spherical Regions Limitations +============================= + +Reading/writing and serializing `~regions.SphericalSkyRegion` objects +is not currently implemented. diff --git a/docs/shapes.rst b/docs/shapes.rst index 7ef42f9f..b2bc36ec 100644 --- a/docs/shapes.rst +++ b/docs/shapes.rst @@ -11,10 +11,15 @@ Region Shapes ============= -Regions provides `~regions.Region` objects, defined in pixel or sky -coordinates, representing shapes such as circles, ellipses, rectangles, -polygons, lines, and points. There are also regions defining circular, -elliptical, and rectangular annuli. +Regions provides `~regions.Region` objects, representing +projected planar "region-in-image" shapes (defined in pixel +or sky coordinates), and also representing spherical +"region-on-celestial-sphere" shapes (defined using sky coordinates). + +The provided `~regions.Region` objects represent shapes such as circles, +ellipses, rectangles, polygons, lines, and points. There are also +regions defining circular, elliptical, and rectangular annuli, +and a longitude/latitude range (in spherical geometry). Defining Shapes @@ -26,34 +31,43 @@ that's currently supported. Circle ****** -`~regions.CircleSkyRegion` and `~regions.CirclePixelRegion` +`~regions.CircleSkyRegion`, `~regions.CirclePixelRegion`, and +`~regions.CircleSphericalSkyRegion` .. code-block:: python >>> from astropy.coordinates import SkyCoord >>> from astropy import units as u >>> from regions import PixCoord - >>> from regions import CircleSkyRegion, CirclePixelRegion + >>> from regions import (CircleSkyRegion, CirclePixelRegion, + ... CircleSphericalSkyRegion) >>> center_sky = SkyCoord(42, 43, unit='deg', frame='fk5') >>> region_sky = CircleSkyRegion(center=center_sky, radius=3 * u.deg) + >>> region_sph_sky = CircleSphericalSkyRegion(center=center_sky, radius=3 * u.deg) >>> region_pix = CirclePixelRegion(center=PixCoord(x=42, y=43), ... radius=4.2) -`~regions.CircleAnnulusSkyRegion` and `~regions.CircleAnnulusPixelRegion` +`~regions.CircleAnnulusSkyRegion`, `~regions.CircleAnnulusPixelRegion`, +and `~regions.CircleAnnulusSphericalSkyRegion` .. code-block:: python >>> from astropy.coordinates import SkyCoord >>> from astropy import units as u >>> from regions import PixCoord - >>> from regions import CircleAnnulusSkyRegion, CircleAnnulusPixelRegion + >>> from regions import (CircleAnnulusSkyRegion, + ... CircleAnnulusPixelRegion, + ... CircleAnnulusSphericalSkyRegion) >>> center_sky = SkyCoord(42, 43, unit='deg', frame='fk5') >>> region_sky = CircleAnnulusSkyRegion(center=center_sky, ... inner_radius=3 * u.deg, ... outer_radius=4 * u.deg) + >>> region_sph_sky = CircleAnnulusSphericalSkyRegion(center=center_sky, + ... inner_radius=3 * u.deg, + ... outer_radius=4 * u.deg) >>> region_pix = CircleAnnulusPixelRegion(center=PixCoord(x=42, y=43), ... inner_radius=4.2, ... outer_radius=5.2) @@ -162,23 +176,46 @@ Polygon .. This is not available yet, for now see `spherical_geometry`_ for .. spherical polygons and `Shapely`_ for pixel polygons. -`~regions.PolygonSkyRegion`, `~regions.PolygonPixelRegion`, and +`~regions.PolygonSkyRegion`, `~regions.PolygonPixelRegion`, +`~regions.PolygonSphericalSkyRegion`, and `~regions.RegularPolygonPixelRegion` .. code-block:: python >>> from astropy.coordinates import SkyCoord - >>> from regions import PixCoord, PolygonSkyRegion, PolygonPixelRegion + >>> from regions import (PixCoord, + ... PolygonSkyRegion, + ... PolygonPixelRegion, + ... PolygonSphericalSkyRegion) >>> from regions import RegularPolygonPixelRegion >>> vertices = SkyCoord([1, 2, 2], [1, 1, 2], unit='deg', frame='fk5') >>> region_sky = PolygonSkyRegion(vertices=vertices) + >>> region_sph_sky = PolygonSphericalSkyRegion(vertices=vertices) >>> vertices = PixCoord(x=[1, 2, 2], y=[1, 1, 2]) >>> region_pix = PolygonPixelRegion(vertices=vertices) >>> center = PixCoord(25, 25) >>> region2_pix = RegularPolygonPixelRegion(center, 6, 15) +Range +***** + +`~regions.RangeSphericalSkyRegion` (Range has no direct analog in +planar geometry. Instead, this shape can be transformed by discretizing +to a polygon and transforming the polygon.) + +.. code-block:: python + + >>> from astropy.coordinates import SkyCoord + >>> from astropy import units as u + >>> from regions import RangeSphericalSkyRegion + + >>> sph_range = RangeSphericalSkyRegion(frame="galactic", + ... longitude_range=[315,45]*u.deg, + ... latitude_range=[0,45]*u.deg) + + Point ***** @@ -232,14 +269,32 @@ The text regions can be used to annotate a text string on an image. Region Transformations ---------------------- -For every region shape there are two classes, one representing a "sky -region" and another representing a "pixel region" on a given image. A -key feature of the regions package is that it is possible to convert +For nearly all region shapes there are three classes. +First, for planar projections (e.g., "regions-on-image", +with Euclidean geometry), +there is one class representing a "sky region" and +another representing a "pixel region" on a given image. +There are also spherical region classes, +representing a "celestial sphere region" (with spherical geometry). + +(The spherical class for some shapes is not currently implemented. +Additionally, some spherical classes do not have projected/planar +analogs, including `~regions.RangeSphericalSkyRegion` +and `~regions.LuneSphericalSkyRegion`.) + +A key feature of the regions package is that it is possible to convert back and forth between sky and image regions given a WCS object (e.g., -`astropy.wcs.WCS`). +`astropy.wcs.WCS`). For conversions to and from spherical sky regions, +it is also necessary to specify how boundary distortions +(from projection effects) should be treated. + -As an example, let's use the :class:`~regions.CircleSkyRegion`, a sky -circle region: +Planar sky and pixel region transformations +******************************************* + + +As an example, let's use the :class:`~regions.CircleSkyRegion`, a +planar sky circle region: .. code-block:: python @@ -276,6 +331,106 @@ to a :class:`~regions.SkyRegion`, call the radius: 18.55481729935556 arcsec +Spherical to planar region transformations +****************************************** + +To demonstrate the transformation from spherical to planar regions, +let's now use a `~regions.CircleSphericalSkyRegion`, a sky circle +defined in spherical geometry: + + +.. code-block:: python + + >>> from regions import CircleSphericalSkyRegion + + >>> center = SkyCoord(50, 10, unit='deg') + >>> radius = Angle(30, 'arcsec') + >>> sph_sky_reg = CircleSphericalSkyRegion(center, radius) + + +To convert it to a planar sky circle region (i.e., +:class:`~regions.CircleSkyRegion`), call the +:meth:`~regions.SphericalSkyRegion.to_sky` method. If boundary distortions +are ignored (i.e., a spherical circle is transformed to a planar circle), +then a WCS object is not necessary: + +.. code-block:: python + + >>> sky_reg = sph_sky_reg.to_sky(include_boundary_distortions=False) + >>> print(sky_reg) # doctest: +FLOAT_CMP + Region: CircleSkyRegion + center: + radius: 30.0 arcsec + +To convert to a planar sky region that accounts for +boundary distortions arising from projection effects, +a WCS object must be passed (defining the projection). +In this case, the spherical region boundary is discretized +into a polygon before transforming, and a `~regions.PolygonSkyRegion` +instance is returned. +The number of points per side in the discretization (for a circle, +the number of points around the circumference) +can be specified through optional keywords that are passed +to the `~regions.SphericalSkyRegion.discretize_boundary()` method. + +.. code-block:: python + + >>> sky_reg = sph_sky_reg.to_sky(wcs=wcs, + ... include_boundary_distortions=True, + ... discretize_kwargs={"n_points": 10}) + >>> print(sky_reg) # doctest: +FLOAT_CMP + Region: PolygonSkyRegion + vertices: + + +Similarly, spherical sky regions can be converted to pixel regions +by specifying a WCS object and whether to ignore or include boundary distortions +(returning either a `~regions.CirclePixelRegion` or `~regions.PolygonPixelRegion`, +respectively). + +.. code-block:: python + + >>> pix_reg = sph_sky_reg.to_pixel(wcs=wcs, + ... include_boundary_distortions=False) + >>> print(pix_reg) # doctest: +FLOAT_CMP + Region: CirclePixelRegion + center: PixCoord(x=55.352057112146014, y=40.095831389269705) + radius: 0.010259141134880476 + + >>> pix_reg2 = sph_sky_reg.to_pixel(wcs=wcs, + ... include_boundary_distortions=True, + ... discretize_kwargs={"n_points": 10}) + >>> print(pix_reg2) # doctest: +FLOAT_CMP + Region: PolygonPixelRegion + vertices: PixCoord(x=[55.35441629 55.36250693 55.366607 55.36514936 + 55.35869012 55.34969718 55.34160657 55.33750862 55.33896755 55.34542545], + y=[40.09920163 40.0934575 40.08862017 40.08653731 40.08800466 40.09246189 + 40.09820641 40.10304382 40.10512634 40.1036587 ]) + + +Planar regions can also be transformed to spherical regions, with the +`~regions.SkyRegion.to_spherical_sky` and `~regions.PixelRegion.to_spherical_sky` +methods. + +.. code-block:: python + + >>> center = SkyCoord(50, 10, unit='deg') + >>> radius = Angle(30, 'arcsec') + >>> sky_reg = CircleSkyRegion(center, radius) + >>> sph_sky_reg = sky_reg.to_spherical_sky(include_boundary_distortions=False) + >>> print(sph_sky_reg) # doctest: +FLOAT_CMP + Region: CircleSphericalSkyRegion + center: + radius: 30.0 arcsec + + .. _regions-as_mpl_selector: Selectors for Regions diff --git a/docs/spherical_bounding_regions.rst b/docs/spherical_bounding_regions.rst new file mode 100644 index 00000000..46868db0 --- /dev/null +++ b/docs/spherical_bounding_regions.rst @@ -0,0 +1,256 @@ + +.. _spherical_bounding_regions: + +Bounding Regions for On-Sky Spatial Searches +============================================ + +The `~regions.SphericalSkyRegion` classes also provide +`~regions.SphericalSkyRegion.bounding_circle` and +`~regions.SphericalSkyRegion.bounding_lonlat` properties, +yielding the spherical circle and longitude and latitude ranges +enclosing the region, respectively. + +This bounding information can be helpful when querying databases +to select objects with tools such as +`Astroquery `_ +or `PyVO's TAP service `_. + +For many regions of interest, the spherical region properties can +directly provide the query search boundary +(e.g., the center and radius for a circle for a cone search, +or the vertices for a polygon region search). +However, this depends on what query shapes are supported by +the database. + +When a particular region shape is not supported by a +database, these bounding properties can be used to define a +"padded" search, and the results can then be downselected to only +those within the region using the `~regions.SphericalSkyRegion.contains()`. + + +(In particular, the bounding properties may be helpful when a +coordinate frame transformation is necessary to match the database +coordinate system; see :ref:`spherical_frame_transform`.) + + +As an example, let's create two spherical regions, a circle and range, +and determine their bounding longitude/latitude spans and circle: + +.. code-block:: python + + >>> from astropy.coordinates import SkyCoord + >>> from astropy import units as u + >>> from regions import CircleSphericalSkyRegion, RangeSphericalSkyRegion + >>> sph_circ = CircleSphericalSkyRegion(SkyCoord(100,-40,unit=u.deg, + ... frame="galactic"), + ... 30*u.deg) + >>> print(sph_circ.bounding_lonlat) # doctest: +FLOAT_CMP + (, ) + >>> sph_range = RangeSphericalSkyRegion(frame="galactic", + ... longitude_range=[315,45]*u.deg, + ... latitude_range=[0,45]*u.deg) + >>> print(sph_range.bounding_circle) # doctest: +FLOAT_CMP + Region: CircleSphericalSkyRegion + center: + radius: 47.2657903991002 deg + +(The bounding circle of a circle is always itself, and in the original +coordinate frame, the bounding longitude/latitude of a +spherical range is equivalent to the region.) + +Let's now transform these regions into the ICRS frame, +and obtain the new bounds, including the new frame longitude/latitude +bounds on the transformed Range region: + + +.. code-block:: python + + >>> sph_circ_transf = sph_circ.transform_to("icrs") + >>> print(sph_circ_transf.bounding_lonlat) # doctest: +FLOAT_CMP + (, ) + >>> sph_range_transf = sph_range.transform_to("icrs") + >>> print(sph_range_transf.bounding_lonlat) # doctest: +FLOAT_CMP + (, ) + >>> print(sph_range_transf.bounding_circle) # doctest: +FLOAT_CMP + Region: CircleSphericalSkyRegion + center: + radius: 47.26579039910021 deg + + +These bounding circles and longitude/latitude are visualized below for both +the original and transformed coordinate frames. + + +.. plot:: + :include-source: false + + from astropy.coordinates import SkyCoord, Angle + from astropy.visualization.wcsaxes.frame import EllipticalFrame + from astropy import units as u + + import matplotlib.pyplot as plt + from matplotlib.lines import Line2D + + from regions import (CircleSphericalSkyRegion, + RangeSphericalSkyRegion, + make_example_dataset) + + + dataset = make_example_dataset(data='simulated') + wcs = dataset.wcs + + dataset_icrs = make_example_dataset(data='simulated', config={'ctype': + ('RA---AIT', + 'DEC--AIT')}) + wcs_icrs = dataset_icrs.wcs + + sph_circ = CircleSphericalSkyRegion(SkyCoord(100,-40, + unit=u.deg, + frame="galactic"), + 30*u.deg) + sph_range = RangeSphericalSkyRegion(frame="galactic", + longitude_range=[315,45]*u.deg, + latitude_range=[0,45]*u.deg) + sph_circ_transf = sph_circ.transform_to("icrs") + sph_range_transf = sph_range.transform_to("icrs") + + + fig = plt.figure() + fig.set_size_inches(7,7) + + axes = [] + axes.append(fig.add_axes([0.15, 0.575, 0.8, 0.425], + projection=wcs, + frame_class=EllipticalFrame)) + axes.append(fig.add_axes([0.15, 0.05, 0.8, 0.425], + projection=wcs_icrs, + frame_class=EllipticalFrame)) + + ax = axes[0] + ax.coords.grid(color='gray') + ax.set_xlabel(r"Galactic $\ell$", labelpad=10) + ax.set_ylabel(r"Galactic $b$") + ax.set_title("Galactic coordinates", pad=5) + + patch = sph_circ.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, color='tab:blue', lw=3) + + sph_range.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":250} + ).plot(ax=ax, color='tab:red', lw=3) + + bound_color = 'tab:blue' + bound_lw = 0.75 + bound_zord = 2 + bound_lon, bound_lat = sph_circ.bounding_lonlat + for i in range(2): + # Need to manually "super sample" to get correct distortion. + # Unfortunately axv/hline works for just "aitoff" projection, + # not doing a WCS it seems. + npt = 250 + xarr = np.repeat(bound_lon[i].degree, npt) + yarr = np.linspace(-90,90,num=npt,endpoint=True) + l1 = Line2D(xarr, yarr, ls='--', color=bound_color, + lw=bound_lw, zorder=bound_zord, + transform=ax.get_transform('galactic')) + xarr = np.linspace(-180,180,num=npt,endpoint=True) + yarr = np.repeat(bound_lat[i].degree, npt) + l2 = Line2D(xarr, yarr, ls='--', color=bound_color, + lw=bound_lw, zorder=bound_zord, + transform=ax.get_transform('galactic')) + ax.add_artist(l1) + ax.add_artist(l2) + + + bound_color = 'tab:red' + sph_range.bounding_circle.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, color='tab:red', lw=bound_lw, ls='--', zorder=bound_zord) + + patch.set_clip_path(ax.coords.frame.patch) + + ax.set_xlim(20,340) + ax.set_ylim(10,170) + + ax = axes[1] + ax.coords.grid(color='gray') + ax.set_xlabel("RA", labelpad=10) + ax.set_ylabel("Dec") + ax.set_title("ICRS coordinates", pad=5) + + patch = sph_circ_transf.to_pixel( + wcs=wcs_icrs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, color='tab:blue', lw=3) + + sph_range_transf.to_pixel( + wcs=wcs_icrs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":250} + ).plot(ax=ax, color='tab:red', lw=3) + + + bound_color = 'tab:blue' + bound_lw = 0.75 + bound_zord = 2 + bound_lon, bound_lat = sph_circ_transf.bounding_lonlat + for i in range(2): + # Need to manually "super sample" to get correct distortion. + # Unfortunately axv/hline works for just "aitoff" projection, + # not doing a WCS it seems. + npt = 250 + xarr = np.repeat(bound_lon[i].degree, npt) + yarr = np.linspace(-90,90,num=npt,endpoint=True) + l1 = Line2D(xarr, yarr, ls='--', color=bound_color, + lw=bound_lw, zorder=bound_zord, + transform=ax.get_transform('icrs')) + xarr = np.linspace(-180,180,num=npt,endpoint=True) + yarr = np.repeat(bound_lat[i].degree, npt) + l2 = Line2D(xarr, yarr, ls='--', color=bound_color, + lw=bound_lw, zorder=bound_zord, + transform=ax.get_transform('icrs')) + ax.add_artist(l1) + ax.add_artist(l2) + + + bound_color = 'tab:red' + bound_lon, bound_lat = sph_range_transf.bounding_lonlat + for i in range(2): + # Need to manually "super sample" to get correct distortion. + # Unfortunately axv/hline works for just "aitoff" projection, + # not doing a WCS it seems. + npt = 250 + xarr = np.repeat(bound_lon[i].degree, npt) + yarr = np.linspace(-90,90,num=npt,endpoint=True) + l1 = Line2D(xarr, yarr, ls='--', color=bound_color, + lw=bound_lw, zorder=bound_zord, + transform=ax.get_transform('icrs')) + xarr = np.linspace(-180,180,num=npt,endpoint=True) + yarr = np.repeat(bound_lat[i].degree, npt) + l2 = Line2D(xarr, yarr, ls='--', color=bound_color, + lw=bound_lw, zorder=bound_zord, + transform=ax.get_transform('icrs')) + ax.add_artist(l1) + ax.add_artist(l2) + sph_range_transf.bounding_circle.to_pixel( + wcs=wcs_icrs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, color='tab:red', lw=bound_lw, ls='--', zorder=bound_zord) + + + patch.set_clip_path(ax.coords.frame.patch) + + ax.set_xlim(20,340) + ax.set_ylim(10,170) + ax.coords[0].set_format_unit(u.deg) diff --git a/docs/spherical_frame_transform.rst b/docs/spherical_frame_transform.rst new file mode 100644 index 00000000..86943c46 --- /dev/null +++ b/docs/spherical_frame_transform.rst @@ -0,0 +1,175 @@ +.. include:: references.txt + +.. _spherical_frame_transform: + +Coordinate Frame Transformation of Spherical Regions +==================================================== + +The `~astropy.coordinates.SkyCoord` class provides functionality for +representing celestial coordinates and for transforming between +different coordinate frames (see the :ref:`astropy-coordinates` documentation +for full details). + +Since the `~regions.SphericalSkyRegion` class represents geometric shapes on the +celestial sphere, these idealized shapes can also be represented in, and +transformed between, any spherical coordinate reference frame. +Transforming spherical regions can be useful in cases where a region of +interest is defined in one coordinate frame, but a query +searching for targets in this region needs to be specified +in second coordinate frame. + +As an example, let's start by defining two sky regions in the Galactic +coordinate frame: + +.. code-block:: python + + >>> from astropy.coordinates import SkyCoord + >>> from astropy import units as u + >>> from regions import CircleSphericalSkyRegion, RangeSphericalSkyRegion + >>> sph_circ = CircleSphericalSkyRegion(SkyCoord(100,-30,unit=u.deg, + ... frame="galactic"), + ... 30*u.deg) + >>> sph_range = RangeSphericalSkyRegion(frame="galactic", + ... longitude_range=[315,45]*u.deg, + ... latitude_range=[0,45]*u.deg) + >>> print(sph_circ) # doctest: +FLOAT_CMP + Region: CircleSphericalSkyRegion + center: + radius: 30.0 deg + >>> print(sph_range) # doctest: +FLOAT_CMP + Region: RangeSphericalSkyRegion + frame: galactic + longitude_range: [315. 45.] deg + latitude_range: [ 0. 45.] deg + + + +To convert these regions to the ICRS coordinate frame, the +`~regions.SphericalSkyRegion.transform_to()` method is used. + +.. code-block:: python + + >>> sph_circ_transf = sph_circ.transform_to("icrs") + >>> sph_range_transf = sph_range.transform_to("icrs") + >>> print(sph_circ_transf) # doctest: +FLOAT_CMP + Region: CircleSphericalSkyRegion + center: + radius: 30.0 deg + >>> print(sph_range_transf) # doctest: +FLOAT_CMP + Region: RangeSphericalSkyRegion + frame: icrs + longitude_bounds: , center_gc2=)> + latitude_bounds: , inner_radius=45.0 deg, outer_radius=90.0 deg)> + + +Note that the Range region boundaries cannot be simply described +by longitude/latitude boundaries in a transformed frame, so +the underlying circular annulus and lune boundaries (capturing +the original frame latitude and longitude bounds, respectively) +are used to describe this Range regions after a transformation. + +These original and transformed regions are shown below on +full-sky projections for the Galactic and ICRS coordinate +reference frames (respectively). + +.. plot:: + :include-source: false + + from astropy.coordinates import SkyCoord, Angle + from astropy.visualization.wcsaxes.frame import EllipticalFrame + from astropy import units as u + + import matplotlib.pyplot as plt + + from regions import (CircleSphericalSkyRegion, + RangeSphericalSkyRegion, + make_example_dataset) + + + dataset = make_example_dataset(data='simulated') + wcs = dataset.wcs + + dataset_icrs = make_example_dataset(data='simulated', config={'ctype': + ('RA---AIT', + 'DEC--AIT')}) + wcs_icrs = dataset_icrs.wcs + + sph_circ = CircleSphericalSkyRegion(SkyCoord(100,-30, + unit=u.deg, + frame="galactic"), + 30*u.deg) + sph_range = RangeSphericalSkyRegion(frame="galactic", + longitude_range=[315,45]*u.deg, + latitude_range=[0,45]*u.deg) + sph_circ_transf = sph_circ.transform_to("icrs") + sph_range_transf = sph_range.transform_to("icrs") + + + fig = plt.figure() + fig.set_size_inches(7,7) + + axes = [] + axes.append(fig.add_axes([0.15, 0.575, 0.8, 0.425], + projection=wcs, + frame_class=EllipticalFrame)) + axes.append(fig.add_axes([0.15, 0.05, 0.8, 0.425], + projection=wcs_icrs, + frame_class=EllipticalFrame)) + + ax = axes[0] + ax.coords.grid(color='black') + ax.set_xlabel(r"Galactic $\ell$", labelpad=10) + ax.set_ylabel(r"Galactic $b$") + ax.set_title("Galactic coordinates", pad=5) + + overlay = ax.get_coords_overlay('icrs') + overlay.grid(color='gray', ls='dotted') + + patch = sph_circ.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, color='tab:blue', lw=3) + + sph_range.to_pixel( + wcs=wcs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":250} + ).plot(ax=ax, color='tab:red', lw=3) + + patch.set_clip_path(ax.coords.frame.patch) + + ax.set_xlim(20,340) + ax.set_ylim(10,170) + + ax = axes[1] + ax.coords.grid(color='gray', ls='dotted') + ax.set_xlabel("RA", labelpad=10) + ax.set_ylabel("Dec") + ax.set_title("ICRS coordinates", pad=5) + + overlay = ax.get_coords_overlay('galactic') + overlay.grid(color='black', ls='solid') + + patch = sph_circ_transf.to_pixel( + wcs=wcs_icrs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":1000} + ).plot(ax=ax, color='tab:blue', lw=3) + + sph_range_transf.to_pixel( + wcs=wcs_icrs, + include_boundary_distortions=True, + discretize_kwargs={"n_points":250} + ).plot(ax=ax, color='tab:red', lw=3) + + patch.set_clip_path(ax.coords.frame.patch) + + ax.set_xlim(20,340) + ax.set_ylim(10,170) + ax.coords[0].set_format_unit(u.deg) diff --git a/regions/_utils/spherical_helpers.py b/regions/_utils/spherical_helpers.py new file mode 100644 index 00000000..a6b3c4c7 --- /dev/null +++ b/regions/_utils/spherical_helpers.py @@ -0,0 +1,569 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +This module provides spherical sky region calculation helper tools. +""" + +import astropy.units as u +import numpy as np +from astropy.coordinates import (Latitude, Longitude, SkyCoord, + SphericalRepresentation, + UnitSphericalRepresentation, + cartesian_to_spherical) + +__all__ = [] + + +def cross_product_skycoord2skycoord(c1, c2): + """ + Compute cross product of two sky coordinates (from a spherical + representation), returning a third sky coordinate (on a spherical + representation). + + Parameters + ---------- + c1 : `~astropy.coordinates.SkyCoord` + The first SkyCoord. + + c2 : `~astropy.coordinates.SkyCoord` + The second SkyCoord. + + Returns + ------- + `~astropy.coordinates.SkyCoord` + The cross product as a SkyCoord, with a spherical representation. + """ + cross = c1.frame.represent_as('cartesian').cross(c2.frame.represent_as('cartesian')) + c_cart = cross / cross.norm() + _, lat, lon = cartesian_to_spherical(c_cart.x, c_cart.y, c_cart.z) + return SkyCoord(lon, lat, frame=c1.frame) + + +def cross_product_sum_skycoord2skycoord(coos): + """ + Cross product sum of vertices, assuming vertices are in CW order. + + Use to determine minimum distance between all vertices, as a + centroid definition. + + Parameters + ---------- + coos : `~astropy.coordinates.SkyCoord` + The vertices as a SkyCoord. + + Returns + ------- + `~astropy.coordinates.SkyCoord` + The minimum distance centroid cross product as a SkyCoord, + with a spherical representation. + """ + # verts are in CW order: + verts_cart = coos.frame.represent_as('cartesian') + crosssum = None + + for i in range(len(verts_cart)): + if crosssum is None: + crosssum = verts_cart[i - 1].cross(verts_cart[i]) + else: + crosssum += verts_cart[i - 1].cross(verts_cart[i]) + + # Normalize sum of cross products to get cartesian representation of + # centroid === minimum distance to other points location + c_cart = crosssum / crosssum.norm() + _, lat, lon = cartesian_to_spherical(c_cart.x, c_cart.y, c_cart.z) + + # Ensure internal data representation format is consistent + # with input coordinate convention: + unit = coos[0].represent_as('spherical').lon.unit + return SkyCoord(lon.to(unit), lat.to(unit), frame=coos.frame) + + +def bounding_lonlat_poles_processing(region, lons_arr, lats_arr, inner_region=None): + """ + Check if region covers either pole & modify latitude bounds + accordingly. + + Parameters + ---------- + region : `~regions.SphericalSkyRegion` + The SphericalSkyRegion region. + + lons_arr : list-like [`~astropy.coordinates.Longitude`] + Initial estimate of longitude bounds. + + lats_arr : list-like [`~astropy.coordinates.Latitude`] + Initial estimate of latitude bounds. + + inner_region : `~regions.SphericalSkyRegion`, optional + The inner region, for an annulus (which modifies the pole logic). + Default is None (for a simply-connected region). + + Returns + ------- + lons_arr : list-like [`~astropy.coordinates.Longitude`] or None + Corrected longitude bounds. Is None if the region goes over + the pole (as the region contains all longitudes). + + lats_arr : list-like [`~astropy.coordinates.Latitude`] + Corrected latitude bounds. + """ + # Check if shape covers either pole & modify lats arr accordingly: + poles = SkyCoord([0, 0], [-90, 90], unit=u.deg, frame=region.frame) + + # ------------------------------------ + # Simply connected region + if inner_region is None: + pole_contains = region.contains(poles) + if np.any(pole_contains): + lons_arr = None + # S pole: + if pole_contains[0]: + lats_arr[0] = -90 * u.deg + # N pole + if pole_contains[1]: + lats_arr[1] = 90 * u.deg + + if lons_arr is not None: + return Longitude(lons_arr), Latitude(lats_arr) + + return None, Latitude(lats_arr) + + # ------------------------------------ + # Inner region set: annulus logic: + pole_contains = inner_region.contains(poles) + if np.any(pole_contains): + lats_raw_inner = get_circle_latitude_tangent_limits(inner_region.center, + inner_region.radius) + + # S pole: + if pole_contains[0]: + # Change lower lats bound to be the minimum of the + # inner boundary itself + lats_arr[0] = lats_raw_inner[0] + # N pole + if pole_contains[1]: + # Change upper lats bound to be the maximum of the + # inner boundary itself + lats_arr[1] = lats_raw_inner[1] + + if lons_arr is not None: + return Longitude(lons_arr), Latitude(lats_arr) + + return None, Latitude(lats_arr) + + +def _get_circle_latitude_tangent_points(center, radius): + # Get the points on an arbitrary circle that + # intersect the tangent latitude circle limits. + + # This includes lon values, with +180deg folding + # for over the pole cases + + crepr = center.represent_as('spherical') + + lons = [crepr.lon, crepr.lon] + lats = [crepr.lat - radius, crepr.lat + radius] + + # Fold if over poles: + if lats[0] < -90 * u.deg: + lats[0] = -180 * u.deg - lats[0] + lons[0] = lons[0] + 180 * u.deg + + if lats[1] > 90 * u.deg: + lats[1] = 180 * u.deg - lats[1] + lons[1] = lons[1] + 180 * u.deg + + return SkyCoord(lons, lats, frame=center.frame) + + +def _get_circle_longitude_tangent_points(center, radius): + # Get the points on an arbitrary circle that + # intersect the tangent longitude circle limits. + + # This includes lat values = 0deg + # (all longitude circles have centers on equator) + + crepr = center.represent_as('spherical') + + lats = [0 * u.deg, 0 * u.deg] + + lats_ref = np.array( + [(crepr.lat - radius).to(u.deg).value, (crepr.lat + radius).to(u.deg).value] + ) + lons = [] + + if np.any(np.abs(lats_ref) > 90): + return None + + if np.any(np.abs(lats_ref) == 90): + return SkyCoord( + [crepr.lon - 90 * u.deg, crepr.lon + 90 * u.deg], lats, frame=center.frame + ) + + for sgn in [-1, 1]: + # Do 1 then -1, because of lon increasing to east + lon_gc = crepr.lon - sgn * ( + np.arccos( + np.sin(radius.to(u.radian).value) / np.cos(crepr.lat.to(u.radian).value) + ) + * u.radian + ).to(u.deg) + + lons.append(lon_gc + sgn * 90 * u.deg) + + return SkyCoord(lons, lats, frame=center.frame) + + +def get_circle_latitude_tangent_limits(center, radius): + """ + For an arbitrary spherical circle with center (lon0, lat0) and + radius R0, get the tangent latitude circle limits (encoded at + latitude values, not radii of the latitude circles). + + Use these to determine the latitude bounding limits from those tangents. + (the points where the latitude circles intersect this circle) + (Tangent circles have Rlat = 90 - (lat0 +- R0), + equivalent to latitude limits lat0 +- R0) + + Ignores any issues with "over the pole" bounds -- + this only computes the values of the tangent circles. + Other processing will handle over the pole bounds logic. + + Parameters + ---------- + center : `~astropy.coordinates.SkyCoord` + The circle center as a SkyCoord. + + radius : `~astropy.coordinates.Angle` or `~astropy.Quantity` + The circle radius (as an |Angle| or |Quantity| with angular units). + + Returns + ------- + latitude_limits : `~astropy.coordinates.Latitude` + Length two |Latitude| with the latitudes of the tangent circles. + """ + tan_lat_pts = _get_circle_latitude_tangent_points(center, radius) + tan_lat_pts = tan_lat_pts.represent_as('spherical') + + return Latitude(tan_lat_pts.lat).to(u.deg) + + +def get_circle_longitude_tangent_limits(center, radius): + """ + For an arbitrary spherical circle with center (lon0, lat0) and + radius R0, get the longitudes of the centers of the tangent + "longitude great circles" (as all longitude circles have + latitude=0deg). + + Use these to determine the longitude bounds for this circle (the + points where the longitude GCs intersect the circle). + + If | lat0 +- R0 | > 90deg: lon_bounds = None (Crosses pole, so no + tangent longitude great circle exists -- all longitude lines [great + circles] cross this circle.) + + If | lat0 +- R0 | = 90deg: lon_bounds = [lon0-90,lon0+90] (Touches + either pole, so spans full 180 of longitudes centered on lon0) + + If | lat0 +- R0 | < 90deg: lon0 +- arccos[ sin(R0) / cos(lat0) ] + + Parameters + ---------- + center : `~astropy.coordinates.SkyCoord` + The circle center as a SkyCoord. + + radius : `~astropy.coordinates.Angle` or `~astropy.Quantity` + The circle radius (as an |Angle| or |Quantity| with angular units). + + Returns + ------- + longitude_limits : `~astropy.coordinates.Longitude` + Length two |Longitude| with the longitudes of the centers + of the tangent longitude great circles. + """ + tan_lon_pts = _get_circle_longitude_tangent_points(center, radius) + + if tan_lon_pts is None: + return None + + tan_lon_pts = tan_lon_pts.represent_as('spherical') + return Longitude(tan_lon_pts.lon).to(u.deg) + + +def _add_tan_pts_if_in_pa_range( + coord_list, + tan_pts, + gc, + wrap_ang, + pas_verts_wrap, + coord=None, +): + pas_tan_pts = gc.center.position_angle(tan_pts).to(u.deg) + + # CHECK RANGES: + # To handle possible cases of lon values "wrapping around" across + # the standard 360->0 wrap, wrap lon values + coord values + # around the first entry angle + pas_tan_pts_wrap = pas_tan_pts.wrap_at(wrap_ang) + + in_range = (pas_tan_pts_wrap >= pas_verts_wrap[0]) & ( + pas_tan_pts_wrap <= pas_verts_wrap[1] + ) + + if np.any(in_range): + tptrepr = tan_pts.represent_as('spherical') + coord_list = np.append(coord_list, getattr(tptrepr, coord)[in_range]) + + return coord_list + + +def _check_edge_lt_pi(pas_verts, wrap_ang): + pas_verts_wrap = pas_verts.wrap_at(wrap_ang) + + is_valid_arc_length = ((pas_verts_wrap[1] - pas_verts_wrap[0]).to(u.deg) <= 180 * u.deg) + + return pas_verts_wrap, is_valid_arc_length + + +def _validate_vertices_ordering(verts, gc, gc_center=None): + if gc is not None: + gc_center = gc.center + pas_verts = gc_center.position_angle(verts).to(u.deg) + + wrap_ang = pas_verts[0] + + # Check ordering of vertices pas: + + # Principle: ALL EDGES must be <= 180 deg of length + # So difference of pa[1] - pa[0] <= 180 deg + # If not, swap the order. + + pas_verts_wrap, is_valid_arc_length = _check_edge_lt_pi(pas_verts, wrap_ang) + + if not is_valid_arc_length: + wrap_ang_opp = pas_verts[-1] + pas_verts_wrap_opp, is_valid_arc_length_opp = _check_edge_lt_pi( + pas_verts[::-1], wrap_ang_opp + ) + if is_valid_arc_length_opp: + return pas_verts_wrap_opp, wrap_ang_opp + + raise ValueError('Invalid arc') # should never occur + + return pas_verts_wrap, wrap_ang + + +def _validate_lon_bounds_ordering(lons_arr, centroid): + # Invert longitude order if centroid is outside of range: + wrap_ang = lons_arr[0] + centroid_lon_wrap = centroid.represent_as('spherical').lon.wrap_at(wrap_ang) + in_range_lon = (centroid_lon_wrap >= lons_arr[0].wrap_at(wrap_ang)) & ( + centroid_lon_wrap <= lons_arr[1].wrap_at(wrap_ang) + ) + if not in_range_lon: + lons_arr = lons_arr[::-1] + + return lons_arr + + +def get_edge_raw_lonlat_bounds_circ_edges(vertices, centroid, gcs): + """ + Get the raw longitude / latitude bounds from the circle edges of + spherical sky region. + + Parameters + ---------- + vertices : `~astropy.coordinates.SkyCoord` + The vertices as a SkyCoord. + + centroid : `~astropy.coordinates.SkyCoord` + The polygon centroid as a SkyCoord. + + gcs : `~regions.CircleSphericalSkyRegion` + The circle boundaries as CircleSphericalSkyRegion instances. + + Returns + ------- + longitude_limits, latitude_limits: `~astropy.coordinates.Longitude` + Length two |Longitude| and |Latitude| with the computed + longitude/latitude bounds from the polygon edges. + """ + # Consider lon/lat of vertices: may produce min/max bounds: + vrepr = vertices.represent_as('spherical') + + # lons_list = vrepr.lon + # lats_list = vrepr.lat + + # Special handling: + # Exclude vertices from longitude bounds if any is on a pole + lons_list = [] + lats_list = [] + for v in vrepr: + if np.abs(v.lat.to(u.deg).deg) < 90: + lons_list.append(v.lon) + lats_list.append(v.lat) + lons_list = Longitude(lons_list, unit=u.radian) + lats_list = Latitude(lats_list, unit=u.radian) + + # Need to also check for "out bulging" from edges, + # as far as latitude/lon bounds: + # eg, 2 vertices at ~60deg: the gc arc goes ~closer to the pole; + # a circle centered close to the pole but not extending over it: + # WIDE lon bounds + + # Requires checking if the "bounding" points are *ON* the polygons. + + for i, gc in enumerate(gcs): + # PAs from gc center to vertices: + verts = SkyCoord(np.concatenate([[vertices[i - 1]], [vertices[i]]])) + + pas_verts_wrap, wrap_ang = _validate_vertices_ordering(verts, gc) + + # -------------------------------------------------------- + # Latitude tangent points from bound circle as len 2 SkyCoord: + # Only add to the list if the tangent point is located along this edge + + tan_lat_pts = _get_circle_latitude_tangent_points(gc.center, gc.radius) + + lats_list = _add_tan_pts_if_in_pa_range( + lats_list, tan_lat_pts, gc, wrap_ang, pas_verts_wrap, + coord='lat' + ) + + # -------------------------------------------------------- + # Longitude tangent points from bound circle as len 2 SkyCoord: + # Only add to the list if the tangent point is located along this edge + + tan_lon_pts = _get_circle_longitude_tangent_points(gc.center, gc.radius) + if tan_lon_pts is not None: + lons_list = _add_tan_pts_if_in_pa_range( + lons_list, tan_lon_pts, gc, wrap_ang, pas_verts_wrap, + coord='lon' + ) + + lons_arr = [lons_list.min(), lons_list.max()] + lats_arr = [lats_list.min(), lats_list.max()] + + # -------------------------------------------------------- + # Invert longitude order if centroid is outside of range: + lons_arr = _validate_lon_bounds_ordering(lons_arr, centroid) + + return Longitude(lons_arr).to(u.deg), Latitude(lats_arr).to(u.deg) + + +def _discretize_edge_boundary(vertices, circ, n_points, + circ_center=None, circ_radius=None): + + if circ is not None: + circ_center = circ.center + circ_radius = circ.radius + + # Discretize an edge boundary defined by a circle, geodetic or not: + # either great circle arc, or a span of a non-great circle + # (e.g., constant lat edges of RangeSphericalSkyRegion) + + # For every edge boundary: determine range of PAs spanned by lines + # connecting circle center to the two vertices bounding that edge: + pas_verts = circ.center.position_angle(vertices).to(u.deg) + + pas_verts_wrap, wrap_ang = _validate_vertices_ordering( + vertices, circ, gc_center=circ_center, + ) + + # Need to wrap angles to calculate span: wrap at lower value: + pas_verts_wrap = pas_verts.wrap_at(wrap_ang) + pa_span = pas_verts_wrap[1] - pas_verts_wrap[0] + + # Sample angle range over pa_span with Npoints, and offset + # to start at pas_verts[0] + theta = np.linspace(0, 1, num=n_points, endpoint=False) * pa_span + pas_verts[0] + + # Calculate directional offsets to get boundary discretization, + # with vertices as SkyCoords + bound_verts = circ.center.directional_offset_by(theta, circ_radius) + + return bound_verts + + +def discretize_line_boundary(coords, n_points): + """ + Discretize all edge boundaries for spherical sky regions. + + Parameters + ---------- + coords : `~astropy.coordinates.SkyCoord` + The start, end of the line as a SkyCoord. + + n_points : int + The number of points for discretization along each edge. + + Returns + ------- + edge_bound_verts : `~astropy.coordinates.SkyCoord` + The discretized boundary edge vertices. + """ + # Iterate over full set of vertices & boundary circles for + # a region (polygon or range) + + gc_center = cross_product_skycoord2skycoord(coords[0], coords[1]) + gc_radius = 90 * u.deg + + bound_verts = _discretize_edge_boundary(coords, None, n_points, + circ_center=gc_center, circ_radius=gc_radius) + + return SkyCoord( + SkyCoord( + bound_verts, + representation_type=UnitSphericalRepresentation, + ), + representation_type=SphericalRepresentation, + ) + + +def discretize_all_edge_boundaries(vertices, circs, n_points): + """ + Discretize all edge boundaries for spherical sky regions. + + Parameters + ---------- + vertices : `~astropy.coordinates.SkyCoord` + The vertices of the spherical sky region. + + circs : `~regions.CircleSphericalSkyRegion` + The circle boundaries as CircleSphericalSkyRegion instances. + + n_points : int + The number of points for discretization along each edge. + + Returns + ------- + all_edge_bound_verts : `~astropy.coordinates.SkyCoord` + The discretized boundary edge vertices. + """ + # Iterate over full set of vertices & boundary circles for + # a region (polygon or range) + + all_edge_bound_verts = None + for i, circ in enumerate(circs): + # PAs from gc center to vertices: + verts = SkyCoord(np.concatenate([[vertices[i - 1]], [vertices[i]]])) + + bound_verts = _discretize_edge_boundary(verts, circ, n_points) + + if all_edge_bound_verts is None: + all_edge_bound_verts = bound_verts + else: + all_edge_bound_verts = SkyCoord(np.concatenate( + [all_edge_bound_verts.copy(), bound_verts] + )) + # For some reason concatenate is adding distances, + # so use remove those by running from + # UnitSpherical->SphericalRepresentation... + all_edge_bound_verts = SkyCoord( + SkyCoord( + all_edge_bound_verts, + representation_type=UnitSphericalRepresentation, + ), + representation_type=SphericalRepresentation, + ) + + return all_edge_bound_verts diff --git a/regions/core/attributes.py b/regions/core/attributes.py index a05f21cd..1accd286 100644 --- a/regions/core/attributes.py +++ b/regions/core/attributes.py @@ -154,6 +154,38 @@ def _validate(self, value): 'scalar angle') +class TwoValAngleorNone(RegionAttribute): + """ + Descriptor class to check that value is a list-like array of length + 2 or None, with angles having either `~astropy.coordinates.Angle` or + `~astropy.units.Quantity` with angular units. + """ + + def _validate(self, value): + if value is None: + pass + elif isinstance(value, Quantity): + if len(value) != 2: + raise ValueError(f'{self.name!r} must be length 2') + + if not value.unit.physical_type == 'angle': + raise ValueError(f'{self.name!r} must have angular units') + else: + try: + if len(value) == 2: + # Check if list-like with individual entries having angular units: + for val in value: + if isinstance(val, Quantity): + if not val.unit.physical_type == 'angle': + raise ValueError(f'{self.name!r} must have angular units') + else: + raise ValueError(f'{self.name!r} must be an angle of length 2') + else: + raise ValueError(f'{self.name!r} must be an angle of length 2') + except TypeError: + raise ValueError(f'{self.name!r} must be an angle of length 2') from None + + class RegionType(RegionAttribute): """ Descriptor class to check the region type of value. diff --git a/regions/core/compound.py b/regions/core/compound.py index 58cabeb7..71204b1b 100644 --- a/regions/core/compound.py +++ b/regions/core/compound.py @@ -1,13 +1,16 @@ # Licensed under a 3-clause BSD style license - see LICENSE.rst import operator as op +import astropy.units as u import numpy as np +from astropy.coordinates import Latitude, Longitude +from regions._utils.spherical_helpers import bounding_lonlat_poles_processing from regions.core.attributes import RegionType -from regions.core.core import PixelRegion, SkyRegion +from regions.core.core import PixelRegion, SkyRegion, SphericalSkyRegion from regions.core.mask import RegionMask -__all__ = ['CompoundPixelRegion', 'CompoundSkyRegion'] +__all__ = ['CompoundPixelRegion', 'CompoundSkyRegion', 'CompoundSphericalSkyRegion'] class CompoundPixelRegion(PixelRegion): @@ -94,6 +97,26 @@ def to_sky(self, wcs): region2=skyreg2, meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + sphreg1 = self.region1.to_spherical_sky( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + sphreg2 = self.region2.to_spherical_sky( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + return CompoundSphericalSkyRegion( + region1=sphreg1, + region2=sphreg2, + operator=self.operator, + meta=self.meta.copy(), + visual=self.visual.copy(), + ) + @staticmethod def _make_annulus_path(patch_inner, patch_outer): """ @@ -252,5 +275,234 @@ def to_pixel(self, wcs): region2=pixreg2, meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + sphreg1 = self.region1.to_spherical_sky( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + sphreg2 = self.region2.to_spherical_sky( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + return CompoundSphericalSkyRegion( + region1=sphreg1, + region2=sphreg2, + operator=self.operator, + meta=self.meta.copy(), + visual=self.visual.copy(), + ) + def as_artist(self, ax, **kwargs): raise NotImplementedError + + +class CompoundSphericalSkyRegion(SphericalSkyRegion): + """ + A class that represents the logical combination of two regions in + spherical sky coordinates. + + Parameters + ---------- + region1 : `~regions.SphericalSkyRegion` + First spherical sky region. + region2 : `~regions.SphericalSkyRegion` + Second spherical sky region. + operator : callable + A callable binary operator. + meta : `~regions.RegionMeta`, optional + A dictionary that stores the meta attributes of this region. + visual : `~regions.RegionVisual`, optional + A dictionary that stores the visual meta attributes of this + region. + """ + + _params = ('region1', 'region2', 'operator') + region1 = RegionType('region1', SphericalSkyRegion) + region2 = RegionType('region2', SphericalSkyRegion) + + def __init__(self, region1, region2, operator, meta=None, visual=None): + if not callable(operator): + raise TypeError('operator must be callable') + + self.region1 = region1 + self.region2 = region2 + if meta is None: + self.meta = region1.meta + else: + self.meta = meta + if visual is None: + self.visual = region1.visual + else: + self.visual = visual + self._operator = operator + + @property + def operator(self): + return self._operator + + @property + def frame(self): + return self.region1.frame + + @property + def bounding_circle(self): + from regions.shapes.circle import CircleSphericalSkyRegion + + if self.operator in [op.or_]: + # Union: + # Compute connection between constituent two bounding circles, + # then total "span" on either side, then half and get center. + + circ1 = self.region1.bounding_circle + circ2 = self.region2.bounding_circle + + pa = circ1.center.position_angle(circ2.center) + + e1 = circ1.center.directional_offset_by(180 * u.deg + pa, circ1.radius) + e2 = circ2.center.directional_offset_by(pa, circ2.radius) + sep = e1.separation(e2) + + cnew = e1.directional_offset_by(pa, sep / 2) + + return CircleSphericalSkyRegion(cnew, sep / 2) + + # Disjoint set / XOR: + # Not clean + # Intersection: + # NOT anywhere near as clean. Much better to define for each subclass + # that has a compound region as their internal region logic. + # TODO: + # General solution requires finding intersections of boundaries, + # and then making a bounding circle. Doable with composed + # circle boundaries & polygon-specific logic from spherical-geometry + # Maybe better to not make bounding_circle and bounding_lonlat + # abstract classes of the base SphericalSkyRegion, and only + # add as properties to individual explicitly defined classes? + raise NotImplementedError + + def _get_edge_raw_lonlat_bounds(self): + if self.operator in [op.or_]: + lons1, lats1 = self.region1.bounding_lonlat + lons2, lats2 = self.region2.bounding_lonlat + if (lons1 is None) | (lons2 is None): + return (None, + Latitude([np.min([lats1[0].to(u.deg).value, + lats2[0].to(u.deg).value]) * u.deg, + np.max([lats1[1].to(u.deg).value, + lats2[1].to(u.deg).value]) * u.deg]).to(u.deg)) + # Both lon ranges defined: + return (Longitude([np.min([lons1[0].to(u.deg).value, + lons2[0].to(u.deg).value]) * u.deg, + np.max([lons1[1].to(u.deg).value, + lons2[1].to(u.deg).value]) * u.deg]).to(u.deg), + Latitude([np.min([lats1[0].to(u.deg).value, + lats2[0].to(u.deg).value]) * u.deg, + np.max([lats1[1].to(u.deg).value, + lats2[1].to(u.deg).value]) * u.deg]).to(u.deg)) + + # XOR, Intersection: not implemented + raise NotImplementedError + + @property + def bounding_lonlat(self): + if self.operator in [op.or_]: + # Union: + # Find min max from constituent bounds + + lons_arr, lats_arr = self._get_edge_raw_lonlat_bounds() + + # Check if shape covers either pole & modify lats arr accordingly: + lons_arr, lats_arr = bounding_lonlat_poles_processing( + self, lons_arr, lats_arr + ) + + return lons_arr, lats_arr + + # Disjoint set / XOR: + # Not clean + # Intersection: + # NOT anywhere near as clean. Much better to define for each subclass + # that has a compound region as their internal region logic. + # TODO: + # General solution requires finding intersections of boundaries, + # and then making a bounding lonlat. Doable with composed + # circle boundaries & polygon-specific logic from spherical-geometry + # Maybe better to not make bounding_circle and bounding_lonlat + # abstract classes of the base SphericalSkyRegion, and only + # add as properties to individual explicitly defined classes? + raise NotImplementedError + + def contains(self, coord): + in_reg = self.operator( + self.region1.contains(coord), self.region2.contains(coord) + ) + if self.meta.get('include', True): + return in_reg + else: + return np.logical_not(in_reg) + + def transform_to(self, frame, merge_attributes=True): + # Transform the consitituent regions then recompose: + return self.__class__( + self.region1.transform_to(frame, merge_attributes=merge_attributes), + self.region2.transform_to(frame, merge_attributes=merge_attributes), + self.operator, + meta=self.meta, + visual=self.visual, + ) + + def discretize_boundary(self, n_points=100): + return CompoundSphericalSkyRegion( + self.region1.discretize_boundary(n_points=n_points), + self.region2.discretize_boundary(n_points=n_points), + operator=self.operator, + meta=self.meta.copy(), + visual=self.visual.copy() + ) + + def to_sky( + self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None, + ): + planarreg1 = self.region1.to_sky( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + planarreg2 = self.region2.to_sky( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + return CompoundSkyRegion( + region1=planarreg1, + region2=planarreg2, + operator=self.operator, + meta=self.meta.copy(), + visual=self.visual.copy(), + ) + + def to_pixel( + self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None, + ): + pixreg1 = self.region1.to_pixel( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + pixreg2 = self.region2.to_pixel( + wcs=wcs, + include_boundary_distortions=include_boundary_distortions, + discretize_kwargs=discretize_kwargs, + ) + return CompoundPixelRegion( + region1=pixreg1, + region2=pixreg2, + operator=self.operator, + meta=self.meta.copy(), + visual=self.visual.copy(), + ) diff --git a/regions/core/core.py b/regions/core/core.py index 5b9ea305..1f6e3c0c 100644 --- a/regions/core/core.py +++ b/regions/core/core.py @@ -4,12 +4,15 @@ import operator import numpy as np +from astropy.coordinates.sky_coordinate_parsers import _get_frame_class +from regions._utils.spherical_helpers import bounding_lonlat_poles_processing from regions.core.metadata import RegionMeta, RegionVisual from regions.core.pixcoord import PixCoord from regions.core.registry import RegionsRegistry -__all__ = ['Region', 'PixelRegion', 'SkyRegion'] +__all__ = ['Region', 'PixelRegion', 'SkyRegion', + 'SphericalSkyRegion', 'ComplexSphericalSkyRegion'] __doctest_skip__ = ['Region.serialize', 'Region.write'] @@ -316,6 +319,39 @@ def to_sky(self, wcs): """ raise NotImplementedError + @abc.abstractmethod + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + """ + Convert to an equivalent spherical `~regions.SphericalSkyRegion` + instance. + + Parameters + ---------- + wcs : `~astropy.wcs.WCS` instance, optional + The world coordinate system transformation to use to convert + between sky and pixel coordinates. Required if transforming + with boundary distortions (if ``include_boundary_distortions`` is True). + Ignored if boundary distortions not included. + + include_boundary_distortions : bool, optional + If True, accounts for boundary boundary distortions in spherical to planar + conversions, by discretizing the boundary and converting the boundary polygon. + Default is False, which converts to an equivalent idealized shape. + + discretize_kwargs : dict, optional + Optional keyword arguments to pass to discretize_boundary() method + if including boundary distortions. + + Returns + ------- + spherical_sky_region : `~regions.SphericalSkyRegion` + A spherical sky region, with an equivalent shape (if + ``include_boundary_distortions`` is False), or a discretized polygon of + the boundary (if ``include_boundary_distortions`` is True). + """ + raise NotImplementedError + @property @abc.abstractmethod def area(self): @@ -529,3 +565,431 @@ def to_pixel(self, wcs): A pixel region. """ raise NotImplementedError + + @abc.abstractmethod + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + """ + Convert to an equivalent spherical `~regions.SphericalSkyRegion` + instance. + + Parameters + ---------- + wcs : `~astropy.wcs.WCS` instance, optional + The world coordinate system transformation to use to convert + between sky and pixel coordinates. Required if transforming + with boundary distortions (if ``include_boundary_distortions`` is True). + Ignored if boundary distortions not included. + + include_boundary_distortions : bool, optional + If True, accounts for boundary boundary distortions in spherical to planar + conversions, by discretizing the boundary and converting the boundary polygon. + Default is False, which converts to an equivalent idealized shape. + + discretize_kwargs : dict, optional + Optional keyword arguments to pass to discretize_boundary() method + if including boundary distortions. + + Returns + ------- + spherical_sky_region : `~regions.SphericalSkyRegion` + A spherical sky region, with an equivalent shape (if + ``include_boundary_distortions`` is False), or a discretized polygon of + the boundary (if ``include_boundary_distortions`` is True). + """ + raise NotImplementedError + + +class SphericalSkyRegion(Region): + """ + Base class for all spherical sky regions (compared to the implicitly + planar SkyRegions). + """ + + def intersection(self, other): + """ + Return a region representing the intersection of this region + with ``other``. + + Parameters + ---------- + other : `~regions.SphericalSkyRegion` + The other region to use for the intersection. + """ + from .compound import CompoundSphericalSkyRegion + + return CompoundSphericalSkyRegion( + region1=self, region2=other, operator=operator.and_ + ) + + def symmetric_difference(self, other): + """ + Return the union of the two regions minus any areas contained in + the intersection of the two regions. + + Parameters + ---------- + other : `~regions.SphericalSkyRegion` + The other region to use for the symmetric difference. + """ + from .compound import CompoundSphericalSkyRegion + + return CompoundSphericalSkyRegion( + region1=self, region2=other, operator=operator.xor + ) + + def union(self, other): + """ + Return a region representing the union of this region with + ``other``. + + Parameters + ---------- + other : `~regions.SphericalSkyRegion` + The other region to use for the union. + """ + from .compound import CompoundSphericalSkyRegion + + return CompoundSphericalSkyRegion( + region1=self, region2=other, operator=operator.or_ + ) + + @staticmethod + def _validate_frame(frame): + # TODO: handle offset origin transformations + + frame = ( + _get_frame_class(frame) if isinstance(frame, str) else frame + ) + return frame + + @property + def frame(self): + """ + Coordinate frame of the region. + """ + if 'center' in self._params: + return self.center.frame + elif 'vertices' in self._params: + return self.vertices[0].frame + else: + raise AttributeError( + "Either 'center' or 'vertices' must be an attribute/property " + 'of the SphericalSkyRegion.' + ) + + @property + @abc.abstractmethod + def bounding_circle(self): + """ + Bounding circle for the spherical sky region. + + Returns + ------- + circle_sph_sky_region: `~regions.CircleSphericalSkyRegion` + A circle spherical sky region object. + """ + raise NotImplementedError + + def _validate_lonlat_bounds(self, lons_arr, lats_arr, inner_region=None): + # Check if shape covers either pole & modify lats arr accordingly: + lons_arr, lats_arr = bounding_lonlat_poles_processing( + self, lons_arr, lats_arr, inner_region=inner_region + ) + return lons_arr, lats_arr + + @property + @abc.abstractmethod + def bounding_lonlat(self): + """ + Bounding longitude and latitude of the spherical sky region, in + the region's frame. + + Returns + ------- + lons_arr : list of `~astropy.coordinates.Longitude` + List of lower, upper boundary longitude values. + + lons_arr : list of `~astropy.coordinates.Latitude` + List of lower, upper boundary latitude values. + """ + raise NotImplementedError + + @abc.abstractmethod + def contains(self, coord): + """ + Check whether a sky coordinate falls inside the spherical sky + region. + + Parameters + ---------- + coord : `~astropy.coordinates.SkyCoord` + The position or positions to check. + """ + raise NotImplementedError + + @abc.abstractmethod + def transform_to(self, frame, merge_attributes=True): + """ + Transform the `SphericalSkyRegion` instance into another + instance with a different coordinate reference frame. + + The precise frame transformed to depends on ``merge_attributes``. + If `False`, the destination frame is used exactly as passed in. + But this is often not quite what one wants. E.g., suppose one wants to + transform an ICRS coordinate that has an obstime attribute to FK4; in + this case, one likely would want to use this information. Thus, the + default for ``merge_attributes`` is `True`, in which the precedence is + as follows: (1) explicitly set (i.e., non-default) values in the + destination frame; (2) explicitly set values in the source; (3) default + value in the destination frame. + + Note that in either case, any explicitly set attributes on the source + |SkyCoord| that are not part of the destination frame's definition are + kept (stored on the resulting |SkyCoord|), and thus one can round-trip + (e.g., from FK4 to ICRS to FK4 without losing obstime). + + Parameters + ---------- + frame : str, or `~astropy.coordinates.BaseCoordinateFrame` class or instance + The frame to transform this coordinate into. + merge_attributes : bool, optional + Whether the default attributes in the destination frame are allowed + to be overridden by explicitly set attributes in the source + (see note above; default: `True`). + + Returns + ------- + sph_sky_region : `~regions.SphericalSkyRegion` + A new spherical sky region represented in the ``frame`` frame. + """ + raise NotImplementedError + + @abc.abstractmethod + def discretize_boundary(self, n_points=10): + """ + Discretize the boundary into a PolygonSphericalSkyRegion, as an + approximation where all sides follow great circles. + + Parameters + ---------- + n_points : int, optional + Number of points along the region's boundary. + + Returns + ------- + poly_sky_region: `~regions.PolygonSphericalSkyRegion` + Spherical sky polygon object. + """ + raise NotImplementedError + + @abc.abstractmethod + def to_sky( + self, wcs=None, include_boundary_distortions=False, discretize_kwargs=None + ): + """ + Convert to a planar `~regions.SkyRegion` instance. + + Parameters + ---------- + wcs : `~astropy.wcs.WCS` instance, optional + The world coordinate system transformation to use to convert + between sky and pixel coordinates. Required if transforming + with boundary distortions (if ``include_boundary_distortions`` is True). + Ignored if boundary distortions not included. + + include_boundary_distortions : bool, optional + If True, accounts for boundary boundary distortions in spherical to planar + conversions, by discretizing the boundary and converting the boundary polygon. + Default is False, which converts to an equivalent idealized shape. + + discretize_kwargs : dict, optional + Optional keyword arguments to pass to discretize_boundary() method + if including boundary distortions. + + Returns + ------- + sky_region : `~regions.SkyRegion` + A planar sky region, with an equivalent shape (if + ``include_boundary_distortions`` is False), or a discretized polygon of + the boundary (if ``include_boundary_distortions`` is True). + """ + raise NotImplementedError + + @abc.abstractmethod + def to_pixel( + self, wcs=None, include_boundary_distortions=False, discretize_kwargs=None + ): + """ + Convert to a planar `~regions.PixelRegion` instance. + + Parameters + ---------- + wcs : `~astropy.wcs.WCS` instance, optional + The world coordinate system transformation to use to convert + between sky and pixel coordinates. Required if transforming + with boundary distortions (if ``include_boundary_distortions`` is True). + Ignored if boundary distortions not included. + + include_boundary_distortions : bool, optional + If True, accounts for boundary boundary distortions in spherical to planar + conversions, by discretizing the boundary and converting the boundary polygon. + Default is False, which converts to an equivalent idealized shape. + + discretize_kwargs : dict, optional + Optional keyword arguments to pass to discretize_boundary() method + if including boundary distortions. + + Returns + ------- + pixel_region : `~regions.PixelRegion` + A pixel region, with an equivalent shape (if + ``include_boundary_distortions`` is False), or a discretized polygon of + the boundary (if ``include_boundary_distortions`` is True). + """ + raise NotImplementedError + + def write(self, filename, format=None, overwrite=False, **kwargs): + # Note explicitly that this is not yet implemented. + raise NotImplementedError + + def serialize(self, format=None, **kwargs): + # Note explicitly that this is not yet implemented. + raise NotImplementedError + + +class ComplexSphericalSkyRegion(SphericalSkyRegion): + """ + Base class for complex cases, where the definitional parameters do + not / cannot transform between coordinate frames (including + RangeSphericalSkyRegion). + """ + + # Because the parameters don't transform, + # these object cannot be defined simply as "compound regions" + + # So some repr and str methods are changed. + + # (For now, this is just RangeSphericalSkyRegion, + # but breaking out these methods/properties to streamline). + + # Instead, a variable set _boundaries is used, containing the + # subclass boundary properties. These boundary properties are, in order of preference, + # (a) accessed from `_[boundname]`, for transformed instances of subclass shapes that + # don't simply map in parameterization from frame to frame (eg, Range) + # (b) derived on-the-fly + + _boundaries = () + + def __repr__(self): + prefix = f'{self.__class__.__name__}' + cls_info = [] + + _do_params_info = False + if (self._params is not None) and ((len(self._params) > 1) | (self._params[0] != 'frame')): + # Only do param-based repr if params other than "frame" are set + # If only param is "frame", make repr with boundary info + _do_params_info = True + + if _do_params_info: + for param in self._params: + if param == 'text': + # place quotes around text value + keyval = f'{param}={getattr(self, param)!r}' + elif param == 'frame': + attr = getattr(self, param) + keyval = f'{param}={getattr(attr, "name", repr(attr))}' + else: + keyval = f'{param}={getattr(self, param)}' + cls_info.append(keyval) + + cls_info = ',\n'.join(cls_info) + return f'<{prefix}(\n{cls_info}\n)>' + else: + # First check if "frame" in self._params: + if (self._params is not None) and self._params[0] == 'frame': + param = 'frame' + attr = getattr(self, param) + keyval = f'{param}={getattr(attr, "name", repr(attr))}' + cls_info.append(keyval) + + # If "params" is None, eg for a "transformed" + # complex shape which is no longer described by the "high-level" + # region parameters, instead give information on the regions + # contained within _boundaries info: + for bound in self._boundaries: + keyval = f'{bound}={repr(getattr(self, bound))}' + cls_info.append(keyval) + + cls_info = ',\n'.join(cls_info) + return f'<{prefix}(\n{cls_info}\n)>' + + def __str__(self): + cls_info = [('Region', self.__class__.__name__)] + + _do_params_info = False + if (self._params is not None) and ((len(self._params) > 1) | (self._params[0] != 'frame')): + # Only do param-based str if params other than "frame" are set. + # If only param is "frame", make str with boundary info + _do_params_info = True + + if _do_params_info: + for param in self._params: + if param in ['text']: + keyval = (param, repr(getattr(self, param))) + + elif param == 'frame': + attr = getattr(self, param) + keyval = (param, getattr(attr, 'name', repr(attr))) + else: + keyval = (param, getattr(self, param)) + cls_info.append(keyval) + return '\n'.join([f'{key}: {val}' for key, val in cls_info]) + else: + # First check if "frame" in self._params: + if (self._params is not None) and (self._params[0] == 'frame'): + param = 'frame' + attr = getattr(self, param) + keyval = (param, getattr(attr, 'name', repr(attr))) + cls_info.append(keyval) + + # If "params" is None, eg for a transformed complex shape, + # instead give information on the regions contained + # within _boundaries info: + for bound in self._boundaries: + keyval = (bound, repr(getattr(self, bound))) + cls_info.append(keyval) + + return '\n'.join([f'{key}: {val}' for key, val in cls_info]) + + @property + @abc.abstractmethod + def _compound_region(self): + """ + Compound region containing composite boundaries for this complex + region shape. + """ + + def copy(self, **changes): + # Boundaries: only copy stashed internal attributes "_{bound name}", + # as otherwise these are derived on-the-fly. + boundaries_interalattr = [f"_{bn}" for bn in list(self._boundaries)] + fields = boundaries_interalattr + [ + '_frame', '_vertices', + '_is_original_frame', '_params', + ] + ['meta', 'visual'] + if self._params is not None: + fields += list(self._params) + + for field in fields: + if (field not in changes) & hasattr(self, field): + changes[field] = copy.deepcopy(getattr(self, field)) + + return self.__class__(**changes) + + @property + def frame(self): + frame = getattr(self, '_frame', None) + if frame is None: + return self._compound_region.frame + return frame diff --git a/regions/core/tests/test_compound.py b/regions/core/tests/test_compound.py index bf720d0e..8522ec3a 100644 --- a/regions/core/tests/test_compound.py +++ b/regions/core/tests/test_compound.py @@ -8,15 +8,25 @@ import astropy.units as u import numpy as np import pytest -from astropy.coordinates import SkyCoord +from astropy.coordinates import Latitude, Longitude, SkyCoord from numpy.testing import assert_allclose -from regions.core import (CompoundPixelRegion, CompoundSkyRegion, PixCoord, +from regions._utils.examples import make_example_dataset +from regions.core import (CompoundPixelRegion, CompoundSkyRegion, + CompoundSphericalSkyRegion, PixCoord, RegionBoundingBox) -from regions.shapes import CirclePixelRegion, CircleSkyRegion +from regions.shapes import (CirclePixelRegion, CircleSkyRegion, + CircleSphericalSkyRegion) from regions.tests.helpers import make_simple_wcs +@pytest.fixture(scope='session', name='wcs') +def fixture_wcs(): + config = dict(crpix=(18, 9), cdelt=(-10, 10), shape=(18, 36)) + dataset = make_example_dataset(config=config) + return dataset.wcs + + class TestCompoundPixel: """ Test compound pixel regions. @@ -60,13 +70,19 @@ def test_rotate(self): assert reg.region2.meta == {} assert reg.region2.visual == {} - def test_union(self): + def test_union(self, wcs): union = self.c1 | self.c2 assert isinstance(union, CompoundPixelRegion) mask = union.to_mask() assert_allclose(mask.data[:, 7], [0, 0, 1, 1, 1, 1, 1, 0, 0]) assert_allclose(mask.data[:, 6], [0, 1, 1, 1, 1, 1, 1, 1, 0]) + union_to_sphsky = union.to_spherical_sky(wcs) + assert isinstance(union_to_sphsky, CompoundSphericalSkyRegion) + + union_to_sky = union.to_sky(wcs) + assert isinstance(union_to_sky, CompoundSkyRegion) + def test_intersection(self): intersection = self.c1 & self.c2 mask = intersection.to_mask() @@ -169,3 +185,78 @@ def test_compound_sky(): meta = {'test_meta': True} compound = CompoundSkyRegion(c1, c2, operator.and_, meta=meta) assert compound.meta['test_meta'] + + union_to_sphsky = union.to_spherical_sky(wcs) + assert isinstance(union_to_sphsky, CompoundSphericalSkyRegion) + + union_to_pixel = union.to_pixel(wcs) + assert isinstance(union_to_pixel, CompoundPixelRegion) + + +def test_compound_spherical_sky(wcs): + skycoord1 = SkyCoord(0 * u.deg, 0 * u.deg, frame='galactic') + c1 = CircleSphericalSkyRegion(skycoord1, 1 * u.deg) + + skycoord2 = SkyCoord(1 * u.deg, 1 * u.deg, frame='galactic') + c2 = CircleSphericalSkyRegion(skycoord2, 0.5 * u.deg) + + test_coord1 = SkyCoord(1.2 * u.deg, 1.2 * u.deg, frame='galactic') + test_coord2 = SkyCoord(0.5 * u.deg, 0.5 * u.deg, frame='galactic') + test_coord3 = SkyCoord(0.7 * u.deg, 0.7 * u.deg, frame='galactic') + test_coord4 = SkyCoord(2 * u.deg, 5 * u.deg, frame='galactic') + + assert c2.contains(test_coord1) and not c1.contains(test_coord1) + assert not c2.contains(test_coord2) and c1.contains(test_coord2) + assert c1.contains(test_coord3) and c2.contains(test_coord3) + assert (not c2.contains(test_coord4) + and not c1.contains(test_coord4)) + + coords = SkyCoord([test_coord1, test_coord2, test_coord3, test_coord4], + frame='galactic') + + union = c1 | c2 + assert (union.contains(coords) == [True, True, True, False]).all() + + intersection = c1 & c2 + assert ((intersection.contains(coords) + == [False, False, True, False]).all()) + + diff = c1 ^ c2 + assert (diff.contains(coords) == [True, True, False, False]).all() + + c3 = CircleSphericalSkyRegion(test_coord4, 0.1 * u.deg) + union = c1 | c2 | c3 + assert (union.contains(coords) == [True, True, True, True]).all() + + intersection = c1 & c2 & c3 + assert ((intersection.contains(coords) + == [False, False, False, False]).all()) + + diff = c1 ^ c2 ^ c3 + assert (diff.contains(coords) == [True, True, False, True]).all() + assert 'Compound' in str(union) + + union_to_sky = union.to_sky(wcs) + assert isinstance(union_to_sky, CompoundSkyRegion) + + union_to_pixel = union.to_pixel(wcs) + assert isinstance(union_to_pixel, CompoundPixelRegion) + + assert isinstance(union.bounding_circle, CircleSphericalSkyRegion) + + bound_lonlat = union.bounding_lonlat + assert isinstance(bound_lonlat[0], Longitude) + assert isinstance(bound_lonlat[1], Latitude) + + assert isinstance(union.transform_to('icrs'), + CompoundSphericalSkyRegion) + + assert isinstance(union.discretize_boundary(n_points=2), + CompoundSphericalSkyRegion) + + skycoord4 = SkyCoord(1 * u.deg, 89 * u.deg, frame='galactic') + c4 = CircleSphericalSkyRegion(skycoord4, 2.5 * u.deg) + union = c1 | c4 + bound_lonlat = union.bounding_lonlat + assert bound_lonlat[0] is None + assert isinstance(bound_lonlat[1], Latitude) diff --git a/regions/shapes/__init__.py b/regions/shapes/__init__.py index 2efdfcad..56b7f37d 100644 --- a/regions/shapes/__init__.py +++ b/regions/shapes/__init__.py @@ -7,7 +7,9 @@ from .circle import * # noqa: F401, F403 from .ellipse import * # noqa: F401, F403 from .line import * # noqa: F401, F403 +from .lune import * # noqa: F401, F403 from .point import * # noqa: F401, F403 from .polygon import * # noqa: F401, F403 +from .range import * # noqa: F401, F403 from .rectangle import * # noqa: F401, F403 from .text import * # noqa: F401, F403 diff --git a/regions/shapes/annulus.py b/regions/shapes/annulus.py index 22e72956..fe22f0e1 100644 --- a/regions/shapes/annulus.py +++ b/regions/shapes/annulus.py @@ -13,17 +13,19 @@ RegionMetaDescr, RegionVisualDescr, ScalarAngle, ScalarPixCoord, ScalarSkyCoord) -from regions.core.compound import CompoundPixelRegion -from regions.core.core import PixelRegion, SkyRegion +from regions.core.compound import (CompoundPixelRegion, + CompoundSphericalSkyRegion) +from regions.core.core import PixelRegion, SkyRegion, SphericalSkyRegion from regions.core.metadata import RegionMeta, RegionVisual from regions.core.pixcoord import PixCoord -from regions.shapes.circle import CirclePixelRegion +from regions.shapes.circle import CirclePixelRegion, CircleSphericalSkyRegion from regions.shapes.ellipse import EllipsePixelRegion, EllipseSkyRegion from regions.shapes.rectangle import RectanglePixelRegion, RectangleSkyRegion -__all__ = ['AnnulusPixelRegion', 'AsymmetricAnnulusPixelRegion', +__all__ = ['AnnulusPixelRegion', 'AnnulusSphericalSkyRegion', + 'AsymmetricAnnulusPixelRegion', 'AsymmetricAnnulusSkyRegion', - 'CircleAnnulusPixelRegion', 'CircleAnnulusSkyRegion', + 'CircleAnnulusPixelRegion', 'CircleAnnulusSkyRegion', 'CircleAnnulusSphericalSkyRegion', 'EllipseAnnulusPixelRegion', 'EllipseAnnulusSkyRegion', 'RectangleAnnulusPixelRegion', 'RectangleAnnulusSkyRegion'] @@ -80,6 +82,43 @@ def rotate(self, center, angle): return self.copy(**changes) +class AnnulusSphericalSkyRegion(SphericalSkyRegion, abc.ABC): + """ + A base class for spherical sky annulus regions. + """ + + @property + def _compound_region(self): + return CompoundSphericalSkyRegion(self._inner_region, self._outer_region, + operator.xor, self.meta, self.visual) + + @property + def frame(self): + return self._inner_region.frame + + @property + def bounding_circle(self): + return self._outer_region.bounding_circle + + @property + def bounding_lonlat(self): + # Bounding lonlat of inner, outer regions + # Check of inner region goes over the pole so + # a more limited lat range is needed + lons_arr, lats_arr = self._outer_region.bounding_lonlat + + # Check if shape covers either pole & modify lats arr accordingly, + # accounting for annular geometry: + lons_arr, lats_arr = self._validate_lonlat_bounds( + lons_arr, lats_arr, inner_region=self._inner_region + ) + + return lons_arr, lats_arr + + def contains(self, coord): + return self._compound_region.contains(coord) + + class CircleAnnulusPixelRegion(AnnulusPixelRegion): """ A circular annulus in pixel coordinates. @@ -157,6 +196,31 @@ def to_sky(self, wcs): return CircleAnnulusSkyRegion(center, inner_radius, outer_radius, self.meta.copy(), self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires planar to spherical projection (using WCS) and discretization + # Will require implementing discretization in pixel space + # to get correct handling of distortions. + raise NotImplementedError + + # ### Potential solution: + # # Leverage polygon class to_spherical_sky() functionality without + # # distortions, as the distortions were already computed in creating + # # that polygon approximation + # return self.discretize_boundary(**discretize_kwargs).to_spherical_sky( + # wcs=wcs, include_boundary_distortions=False + # ) + + return self.to_sky(wcs).to_spherical_sky() + class CircleAnnulusSkyRegion(SkyRegion): """ @@ -205,6 +269,171 @@ def to_pixel(self, wcs): meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires planar to spherical projection (using WCS) and discretization + # Will require implementing discretization in pixel space + # to get correct handling of distortions. + raise NotImplementedError + + # ### Potential solution: + # # Leverage polygon class to_spherical_sky() functionality without + # # distortions, as the distortions were already computed in creating + # # that polygon approximation + # return self.to_pixel(wcs).discretize_boundary(**discretize_kwargs).to_spherical_sky( + # wcs=wcs, include_boundary_distortions=False + # ) + + return CircleAnnulusSphericalSkyRegion( + self.center, self.inner_radius, self.outer_radius, + self.meta.copy(), self.visual.copy() + ) + + +class CircleAnnulusSphericalSkyRegion(AnnulusSphericalSkyRegion): + """ + Class for a circular annulus sky region, where the circular annulus + is interpreted within a spherical geometry reference frame. + + Parameters + ---------- + center : `~astropy.coordinates.SkyCoord` + The center position. + inner_radius : `~astropy.units.Quantity` + The inner radius in angular units. + outer_radius : `~astropy.units.Quantity` + The outer radius in angular units. + meta : `~regions.RegionMeta` or `dict`, optional + A dictionary that stores the meta attributes of the region. + visual : `~regions.RegionVisual` or `dict`, optional + A dictionary that stores the visual meta attributes of the + region. + """ + + _component_class = CircleSphericalSkyRegion + _params = ('center', 'inner_radius', 'outer_radius') + center = ScalarSkyCoord('The center position as a |SkyCoord|.') + inner_radius = PositiveScalarAngle('The inner radius as a |Quantity| ' + 'angle.') + outer_radius = PositiveScalarAngle('The outer radius as a |Quantity| ' + 'angle.') + meta = RegionMetaDescr('The meta attributes as a |RegionMeta|') + visual = RegionVisualDescr('The visual attributes as a |RegionVisual|.') + + def __init__(self, center, inner_radius, outer_radius, meta=None, + visual=None): + self.center = center + self.inner_radius = inner_radius + self.outer_radius = outer_radius + self.meta = meta or RegionMeta() + self.visual = visual or RegionVisual() + + if inner_radius >= outer_radius: + raise ValueError('outer_radius must be greater than inner_radius') + + @property + def _inner_region(self): + return self._component_class(self.center, self.inner_radius, + self.meta, self.visual) + + @property + def _outer_region(self): + return self._component_class(self.center, self.outer_radius, + self.meta, self.visual) + + def transform_to(self, frame, merge_attributes=True): + frame = self._validate_frame(frame) + + # Only center transforms, radii preserved + center_transf = self.center.transform_to(frame, merge_attributes=merge_attributes) + + return CircleAnnulusSphericalSkyRegion( + center_transf, + self.inner_radius.copy(), + self.outer_radius.copy(), + self.meta.copy(), + self.visual.copy() + ) + + def discretize_boundary(self, n_points=100): + return CompoundSphericalSkyRegion( + self._inner_region.discretize_boundary(n_points=n_points), + self._outer_region.discretize_boundary(n_points=n_points), + operator=operator.xor, + meta=self.meta.copy(), + visual=self.visual.copy() + ) + + def to_sky( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None + ): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + # Use to_pixel(), then apply "small angle approx" to get planar sky. + return self.to_pixel( + include_boundary_distortions=include_boundary_distortions, + wcs=wcs, + discretize_kwargs=discretize_kwargs, + ).to_sky(wcs) + + return CircleAnnulusSkyRegion( + self.center.copy(), + self.inner_radius.copy(), + self.outer_radius.copy(), + meta=self.meta.copy(), + visual=self.visual.copy() + ) + + def to_pixel( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None, + ): + if include_boundary_distortions: + from .polygon import PolygonPixelRegion + + if discretize_kwargs is None: + discretize_kwargs = {} + + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + polygonized = self.discretize_boundary(**discretize_kwargs) + + inner_vertices = wcs.world_to_pixel(polygonized.region1.vertices) + outer_vertices = wcs.world_to_pixel(polygonized.region2.vertices) + + return CompoundPixelRegion( + PolygonPixelRegion(PixCoord(*inner_vertices)), + PolygonPixelRegion(PixCoord(*outer_vertices)), + operator=operator.xor, + meta=self.meta.copy(), + visual=self.visual.copy() + ) + + return self.to_sky().to_pixel(wcs) + class AsymmetricAnnulusPixelRegion(AnnulusPixelRegion): """ @@ -458,6 +687,10 @@ def to_sky(self, wcs): meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + class EllipseAnnulusSkyRegion(AsymmetricAnnulusSkyRegion): """ @@ -516,6 +749,10 @@ def to_pixel(self, wcs): meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + class RectangleAnnulusPixelRegion(AsymmetricAnnulusPixelRegion): """ @@ -599,6 +836,10 @@ def to_sky(self, wcs): meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + class RectangleAnnulusSkyRegion(AsymmetricAnnulusSkyRegion): """ @@ -657,3 +898,7 @@ def to_pixel(self, wcs): return RectangleAnnulusPixelRegion(*self.to_pixel_args(wcs), meta=self.meta.copy(), visual=self.visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError diff --git a/regions/shapes/circle.py b/regions/shapes/circle.py index c06d33dd..0b131c3b 100644 --- a/regions/shapes/circle.py +++ b/regions/shapes/circle.py @@ -10,17 +10,19 @@ from astropy.coordinates import Angle from regions._geometry import circular_overlap_grid +from regions._utils.spherical_helpers import ( + get_circle_latitude_tangent_limits, get_circle_longitude_tangent_limits) from regions._utils.wcs_helpers import pixel_scale_angle_at_skycoord from regions.core.attributes import (PositiveScalar, PositiveScalarAngle, RegionMetaDescr, RegionVisualDescr, ScalarPixCoord, ScalarSkyCoord) from regions.core.bounding_box import RegionBoundingBox -from regions.core.core import PixelRegion, SkyRegion +from regions.core.core import PixelRegion, SkyRegion, SphericalSkyRegion from regions.core.mask import RegionMask from regions.core.metadata import RegionMeta, RegionVisual from regions.core.pixcoord import PixCoord -__all__ = ['CirclePixelRegion', 'CircleSkyRegion'] +__all__ = ['CirclePixelRegion', 'CircleSkyRegion', 'CircleSphericalSkyRegion'] class CirclePixelRegion(PixelRegion): @@ -92,6 +94,31 @@ def to_sky(self, wcs): return CircleSkyRegion(center, radius, meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires planar to spherical projection (using WCS) and discretization + # Will require implementing discretization in pixel space + # to get correct handling of distortions. + raise NotImplementedError + + # ### Potential solution: + # # Leverage polygon class to_spherical_sky() functionality without + # # distortions, as the distortions were already computed in creating + # # that polygon approximation + # return self.discretize_boundary(**discretize_kwargs).to_spherical_sky( + # wcs=wcs, include_boundary_distortions=False + # ) + + return self.to_sky(wcs).to_spherical_sky() + @property def bounding_box(self): """ @@ -215,3 +242,155 @@ def to_pixel(self, wcs): radius = (self.radius / pixscale).to(u.pix).value return CirclePixelRegion(center, radius, meta=self.meta.copy(), visual=self.visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires planar to spherical projection (using WCS) and discretization + # Will require implementing discretization in pixel space + # to get correct handling of distortions. + raise NotImplementedError + + # ### Potential solution: + # # Leverage polygon class to_spherical_sky() functionality without + # # distortions, as the distortions were already computed in creating + # # that polygon approximation + # return self.to_pixel(wcs).discretize_boundary(**discretize_kwargs).to_spherical_sky( + # wcs=wcs, include_boundary_distortions=False + # ) + + return CircleSphericalSkyRegion( + self.center, self.radius, meta=self.meta, visual=self.visual + ) + + +class CircleSphericalSkyRegion(SphericalSkyRegion): + """ + Class for a circular sky region, where the circle is interpreted + within a spherical geometry reference frame. + + This region is very much akin to the `~regions.CircleSkyRegion` + class (and borrows internal attributes following the same structure). + + Parameters + ---------- + center : `~astropy.coordinates.SkyCoord` + The center position. + radius : `~astropy.units.Quantity` + The radius in angular units. + meta : `~regions.RegionMeta` or `dict`, optional + A dictionary that stores the meta attributes of the region. + visual : `~regions.RegionVisual` or `dict`, optional + A dictionary that stores the visual meta attributes of the + region. + """ + + _params = ('center', 'radius') + center = ScalarSkyCoord('The center position as a |SkyCoord|.') + radius = PositiveScalarAngle('The radius as a |Quantity| angle.') + meta = RegionMetaDescr('The meta attributes as a |RegionMeta|') + visual = RegionVisualDescr('The visual attributes as a |RegionVisual|.') + + def __init__(self, center, radius, meta=None, visual=None): + self.center = center + self.radius = radius + self.meta = meta or RegionMeta() + self.visual = visual or RegionVisual() + + def contains(self, coord): + in_circle = self.center.separation(coord) < self.radius + if self.meta.get('include', True): + return in_circle + else: + return np.logical_not(in_circle) + + @property + def bounding_circle(self): + return self.copy() + + @property + def bounding_lonlat(self): + lons_arr = get_circle_longitude_tangent_limits(self.center, self.radius) + lats_arr = get_circle_latitude_tangent_limits(self.center, self.radius) + + lons_arr, lats_arr = self._validate_lonlat_bounds(lons_arr, lats_arr) + + return lons_arr, lats_arr + + def transform_to(self, frame, merge_attributes=True): + frame = self._validate_frame(frame) + + center_transf = self.center.transform_to( + frame, merge_attributes=merge_attributes + ) + + return CircleSphericalSkyRegion( + center_transf, + self.radius.copy(), + self.meta.copy(), + self.visual.copy() + ) + + def discretize_boundary(self, n_points=100): + # Avoid circular imports: + from .polygon import PolygonSphericalSkyRegion + theta = np.linspace(0, 1, num=n_points, endpoint=False) * 360 * u.deg + # Need to invert order because of CW convention: + bound_verts = self.center.directional_offset_by(theta[::-1], self.radius) + return PolygonSphericalSkyRegion(bound_verts) + + def to_sky( + self, wcs=None, include_boundary_distortions=False, discretize_kwargs=None + ): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + # Use to_pixel(), then apply "small angle approx" to get planar sky. + return self.to_pixel( + include_boundary_distortions=include_boundary_distortions, + wcs=wcs, + discretize_kwargs=discretize_kwargs, + ).to_sky(wcs) + + return CircleSkyRegion( + self.center, self.radius, meta=self.meta, visual=self.visual + ) + + def to_pixel( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None, + ): + if include_boundary_distortions: + from .polygon import PolygonPixelRegion + + if discretize_kwargs is None: + discretize_kwargs = {} + + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + verts = wcs.world_to_pixel( + self.discretize_boundary(**discretize_kwargs).vertices + ) + return PolygonPixelRegion( + PixCoord(*verts), meta=self.meta.copy(), visual=self.visual.copy() + ) + + return self.to_sky().to_pixel(wcs) diff --git a/regions/shapes/ellipse.py b/regions/shapes/ellipse.py index c4bfca37..6e6669ca 100644 --- a/regions/shapes/ellipse.py +++ b/regions/shapes/ellipse.py @@ -121,6 +121,10 @@ def to_sky(self, wcs): meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + @property def bounding_box(self): """ @@ -376,3 +380,7 @@ def to_pixel(self, wcs): return EllipsePixelRegion(center, width, height, angle=angle, meta=self.meta.copy(), visual=self.visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError diff --git a/regions/shapes/line.py b/regions/shapes/line.py index 3e045117..303b2ef3 100644 --- a/regions/shapes/line.py +++ b/regions/shapes/line.py @@ -84,6 +84,10 @@ def to_sky(self, wcs): return LineSkyRegion(start, end, meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + @property def bounding_box(self): xmin = min(self.start.x, self.end.x) @@ -197,3 +201,7 @@ def to_pixel(self, wcs): end = PixCoord(end_x, end_y) return LinePixelRegion(start, end, meta=self.meta.copy(), visual=self.visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError diff --git a/regions/shapes/lune.py b/regions/shapes/lune.py new file mode 100644 index 00000000..bbfcdfca --- /dev/null +++ b/regions/shapes/lune.py @@ -0,0 +1,240 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +This module defines a lune (a "wedge" bounded by two great circles) +spherical sky region. +""" + +import operator + +import astropy.units as u +from astropy.coordinates import SkyCoord + +from regions._utils.spherical_helpers import ( + cross_product_skycoord2skycoord, cross_product_sum_skycoord2skycoord, + discretize_all_edge_boundaries, get_edge_raw_lonlat_bounds_circ_edges) +from regions.core.attributes import (RegionMetaDescr, RegionVisualDescr, + ScalarSkyCoord) +from regions.core.compound import CompoundSphericalSkyRegion +from regions.core.core import SphericalSkyRegion +from regions.core.metadata import RegionMeta, RegionVisual +from regions.shapes.circle import CircleSphericalSkyRegion +from regions.shapes.polygon import PolygonSphericalSkyRegion + +__all__ = ['LuneSphericalSkyRegion'] + + +class LuneSphericalSkyRegion(SphericalSkyRegion): + """ + Class for a spherical "lune", the intersection between to great + circles in spherical geometry. + + Composed of two CircleSphericalSkyRegions. + + Parameters + ---------- + center_gc1 : `~astropy.coordinates.SkyCoord` + The center position of the first great circle. + center_gc2 : `~astropy.coordinates.SkyCoord` + The center position of the second great circle. + meta : `~regions.RegionMeta` or `dict`, optional + A dictionary that stores the meta attributes of the region. + visual : `~regions.RegionVisual` or `dict`, optional + A dictionary that stores the visual meta attributes of the + region. + """ + + _params = ('center_gc1', 'center_gc2') + center_gc1 = ScalarSkyCoord('The center position of the first great circle as a |SkyCoord|.') + center_gc2 = ScalarSkyCoord('The center position of the second great circle as a |SkyCoord|.') + + meta = RegionMetaDescr('The meta attributes as a |RegionMeta|') + visual = RegionVisualDescr('The visual attributes as a |RegionVisual|.') + + def __init__(self, center_gc1, center_gc2, meta=None, visual=None): + self.center_gc1 = center_gc1 + self.center_gc2 = center_gc2 + + self.meta = meta or RegionMeta() + self.visual = visual or RegionVisual() + + @property + def _circle_1(self): + return CircleSphericalSkyRegion( + self.center_gc1, 90 * u.deg, + self.meta, self.visual + ) + + @property + def _circle_2(self): + return CircleSphericalSkyRegion( + self.center_gc2, 90 * u.deg, + self.meta, self.visual + ) + + @property + def _edge_circs(self): + """ + Get list of the great circles defining the lune boundaries. + """ + return [self._circle_1, self._circle_2] + + @property + def _compound_region(self): + return CompoundSphericalSkyRegion(self._circle_1, self._circle_2, + operator.and_, self.meta, self.visual) + + @property + def frame(self): + return self.center_gc1.frame + + @property + def centroid(self): + """ + Region centroid. + """ + # Cross product to get poles: + c_pole = cross_product_skycoord2skycoord(self.center_gc2, self.center_gc1) + c_pole_a = cross_product_skycoord2skycoord(self.center_gc1, self.center_gc2) + + # Use centers + poles to get centroid from cross products of vertices + # with cartesian representation + # Assume CW order is gc1, c_pole, gc2, c_pole_a + + c_gc1_sph = self.center_gc1.represent_as('spherical') + c_gc2_sph = self.center_gc2.represent_as('spherical') + c_p_sph = c_pole.represent_as('spherical') + c_p_a_sph = c_pole_a.represent_as('spherical') + + verts = SkyCoord( + [ + c_gc2_sph.lon, + c_p_a_sph.lon, + c_gc1_sph.lon, + c_p_sph.lon, + ], + [ + c_gc2_sph.lat, + c_p_a_sph.lat, + c_gc1_sph.lat, + c_p_sph.lat, + ], + frame=self.frame + ) + + return cross_product_sum_skycoord2skycoord(verts) + + @property + def vertices(self): + """ + Region vertices. + + For a lune, these are the intersections of the two bounding + great circles, and are antipodes. + """ + # Cross product to get poles: + c_pole = cross_product_skycoord2skycoord(self.center_gc2, self.center_gc1) + c_pole_a = cross_product_skycoord2skycoord(self.center_gc1, self.center_gc2) + + c_p_sph = c_pole.represent_as('spherical') + c_p_a_sph = c_pole_a.represent_as('spherical') + + verts = SkyCoord( + [ + c_p_sph.lon, + c_p_a_sph.lon, + ], + [ + c_p_sph.lat, + c_p_a_sph.lat, + ], + frame=self.frame + ) + + return verts + + @property + def bounding_circle(self): + return CircleSphericalSkyRegion(self.centroid, 90 * u.deg) + + @property + def bounding_lonlat(self): + lons_arr, lats_arr = get_edge_raw_lonlat_bounds_circ_edges( + self.vertices, self.centroid, self._edge_circs + ) + + lons_arr, lats_arr = self._validate_lonlat_bounds(lons_arr, lats_arr) + + return lons_arr, lats_arr + + def contains(self, coord): + return self._compound_region.contains(coord) + + def transform_to(self, frame, merge_attributes=True): + frame = self._validate_frame(frame) + + center_gc1_transf = self.center_gc1.transform_to(frame, merge_attributes=merge_attributes) + center_gc2_transf = self.center_gc2.transform_to(frame, merge_attributes=merge_attributes) + + return LuneSphericalSkyRegion( + center_gc1_transf, + center_gc2_transf, + self.meta.copy(), + self.visual.copy() + ) + + def discretize_boundary(self, n_points=100): + bound_verts = discretize_all_edge_boundaries( + self.vertices, self._edge_circs, n_points + ) + return PolygonSphericalSkyRegion(bound_verts) + + def to_sky( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None + ): + if not include_boundary_distortions: + raise ValueError( + 'Invalid parameter: `include_boundary_distortions=False`!\n' + 'Transforming lune to planar sky region is only possible by ' + 'including boundary distortions, as there is no analogous sky region.' + ) + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + # Use to_pixel(), then apply "small angle approx" to get planar sky. + return self.to_pixel( + include_boundary_distortions=include_boundary_distortions, + wcs=wcs, + discretize_kwargs=discretize_kwargs, + ).to_sky(wcs) + + def to_pixel( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None, + ): + if not include_boundary_distortions: + raise ValueError( + 'Invalid parameter: `include_boundary_distortions=False`!\n' + 'Transforming range to planar pixel region is only possible by ' + 'including boundary distortions, as there is no analogous pixel region.' + ) + + if discretize_kwargs is None: + discretize_kwargs = {} + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + raise NotImplementedError diff --git a/regions/shapes/point.py b/regions/shapes/point.py index 9aecece1..ef083c19 100644 --- a/regions/shapes/point.py +++ b/regions/shapes/point.py @@ -89,6 +89,10 @@ def to_sky(self, wcs): return PointSkyRegion(center, meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + @property def bounding_box(self): return RegionBoundingBox.from_float(self.center.x, self.center.x, @@ -183,3 +187,7 @@ def to_pixel(self, wcs): center = PixCoord(center_x, center_y) return PointPixelRegion(center, meta=self.meta.copy(), visual=self.visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError diff --git a/regions/shapes/polygon.py b/regions/shapes/polygon.py index 8889564b..12e70c75 100644 --- a/regions/shapes/polygon.py +++ b/regions/shapes/polygon.py @@ -2,24 +2,31 @@ """ This module defines polygon regions in both pixel and sky coordinates. """ +import operator import astropy.units as u import numpy as np +from astropy.coordinates import SkyCoord +from astropy.stats import circmean from regions._geometry import polygonal_overlap_grid from regions._geometry.pnpoly import points_in_polygon +from regions._utils.spherical_helpers import ( + cross_product_skycoord2skycoord, cross_product_sum_skycoord2skycoord, + discretize_all_edge_boundaries, get_edge_raw_lonlat_bounds_circ_edges) from regions.core.attributes import (OneDPixCoord, OneDSkyCoord, PositiveScalar, RegionMetaDescr, RegionVisualDescr, ScalarAngle, ScalarPixCoord) from regions.core.bounding_box import RegionBoundingBox -from regions.core.core import PixelRegion, SkyRegion +from regions.core.compound import CompoundSphericalSkyRegion +from regions.core.core import PixelRegion, SkyRegion, SphericalSkyRegion from regions.core.mask import RegionMask from regions.core.metadata import RegionMeta, RegionVisual from regions.core.pixcoord import PixCoord __all__ = ['PolygonPixelRegion', 'RegularPolygonPixelRegion', - 'PolygonSkyRegion'] + 'PolygonSkyRegion', 'PolygonSphericalSkyRegion'] class PolygonPixelRegion(PixelRegion): @@ -112,6 +119,35 @@ def to_sky(self, wcs): return PolygonSkyRegion(vertices=vertices_sky, meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires planar to spherical projection (using WCS) and discretization + # Will require implementing discretization in pixel space + # to get correct handling of distortions. + raise NotImplementedError + + # ### Potential solution: + # # Leverage polygon class to_spherical_sky() functionality without + # # distortions, as the distortions were already computed in creating + # # that polygon approximation + # return self.to_pixel(wcs).discretize_boundary(**discretize_kwargs).to_spherical_sky( + # wcs=wcs, include_boundary_distortions=False + # ) + + # TODO: ensure vertices are in CW order, + # as implicitly a planar -> spherical polygon will not be the + # "large" complement polygon on the sphere? + + return self.to_sky(wcs).to_spherical_sky() + @property def bounding_box(self): xmin = self.vertices.x.min() @@ -383,3 +419,239 @@ def to_pixel(self, wcs): vertices_pix = PixCoord(x, y) return PolygonPixelRegion(vertices_pix, meta=self.meta.copy(), visual=self.visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires planar to spherical projection (using WCS) and discretization + # Will require implementing discretization in pixel space + # to get correct handling of distortions. + raise NotImplementedError + + # ### Potential solution: + # # Leverage polygon class to_spherical_sky() functionality without + # # distortions, as the distortions were already computed in creating + # # that polygon approximation + # return self.to_pixel(wcs).discretize_boundary(**discretize_kwargs).to_spherical_sky( + # wcs=wcs, include_boundary_distortions=False + # ) + + # TODO: ensure vertices are in CW order, + # as implicitly a planar -> spherical polygon will not be the + # "large" complement polygon on the sphere? + + return PolygonSphericalSkyRegion( + self.vertices, + self.meta.copy(), + self.visual.copy() + ) + + +class PolygonSphericalSkyRegion(SphericalSkyRegion): + """ + A spherical polygon defined using vertices in sky coordinates. + + Internal "contains" logic could be changed to + leverage spherical-geometry package functionality. + + This region is very much akin to the `~regions.PolygonPixelRegion` + class (and borrows internal attributes following the same structure). + + Parameters + ---------- + vertices : `~astropy.coordinates.SkyCoord` + The vertices of the polygon. Assumed to be in CW order. + meta : `~regions.RegionMeta` or `dict`, optional + A dictionary that stores the meta attributes of the region. + visual : `~regions.RegionVisual` or `dict`, optional + A dictionary that stores the visual meta attributes of the + region. + """ + + _params = ('vertices',) + vertices = OneDSkyCoord('The vertices of the polygon as a |SkyCoord| ' + 'array.') + meta = RegionMetaDescr('The meta attributes as a |RegionMeta|') + visual = RegionVisualDescr('The visual attributes as a |RegionVisual|.') + + def __init__(self, vertices, meta=None, visual=None): + self.vertices = vertices + self.meta = meta or RegionMeta() + self.visual = visual or RegionVisual() + + @property + def _edge_circs(self): + """ + Get list of the great circles defining the polygon boundaries. + """ + from .circle import CircleSphericalSkyRegion + gcs = [] + for i in range(len(self.vertices)): + # verts are in CW order: Cross product to get bounding great circle centers + # Compute GCs and stack into a compound set: + c_gc = cross_product_skycoord2skycoord(self.vertices[i - 1], self.vertices[i]) + gcs.append(CircleSphericalSkyRegion(c_gc, 90 * u.deg)) + return gcs + + @property + def _compound_region(self): + # Need N great circles to define boundaries for an N-sided polygon. + # verts are in CW order: Cross product to get bounding great circle centers + # Compute GCs and stack into a compound set: + compreg = None + gcs = self._edge_circs + for gc in gcs: + if compreg is None: + compreg = gc + else: + compreg = CompoundSphericalSkyRegion( + compreg, gc, operator.and_, self.meta, self.visual, + ) + + return compreg + + @property + def centroid(self): + """ + Region centroid. + + Defined as the point equidistant from all vertices. + + However, if this is not contained within the polygon, instead + use the average of the vertices' positions. + """ + # Calculate from cross products of vertices with cartesian representation + # verts are in CW order: + # Minimum distance + centroid_mindist = self.centroid_mindist + + if not self.contains(centroid_mindist): + return self.centroid_avg + + return centroid_mindist + + @property + def centroid_mindist(self): + """ + Region centroid, defined as the point equidistant from all + vertices (minimum distance). + """ + # Calculate from cross products of vertices with cartesian representation + # verts are in CW order: + # Minimum distance + return cross_product_sum_skycoord2skycoord(self.vertices) + + @property + def centroid_avg(self): + """ + Region centroid, taking the average of the vertices' + coordinates. + """ + verts_sph = self.vertices.represent_as('spherical') + + lons = verts_sph.lon + lats = verts_sph.lat + + lon = circmean(lons) + lat = np.average(lats) + return SkyCoord(lon, lat, frame=self.frame) + + @property + def bounding_circle(self): + from .circle import CircleSphericalSkyRegion + cent = self.centroid + seps = cent.separation(self.vertices) + return CircleSphericalSkyRegion(center=cent, radius=np.max(seps)) + + @property + def bounding_lonlat(self): + lons_arr, lats_arr = get_edge_raw_lonlat_bounds_circ_edges( + self.vertices, self.centroid, self._edge_circs + ) + + lons_arr, lats_arr = self._validate_lonlat_bounds(lons_arr, lats_arr) + + return lons_arr, lats_arr + + def contains(self, coord): + return self._compound_region.contains(coord) + + def transform_to(self, frame, merge_attributes=True): + frame = self._validate_frame(frame) + + # Only center transforms, radii preserved + verts_transf = self.vertices.transform_to(frame, merge_attributes=merge_attributes) + + return PolygonSphericalSkyRegion( + verts_transf, + self.meta.copy(), + self.visual.copy() + ) + + def discretize_boundary(self, n_points=10): + bound_verts = discretize_all_edge_boundaries( + self.vertices, self._edge_circs, n_points + ) + return PolygonSphericalSkyRegion(bound_verts) + + def to_sky( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None + ): + + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + # Use to_pixel(), then apply "small angle approx" to get planar sky. + return self.to_pixel( + include_boundary_distortions=include_boundary_distortions, + wcs=wcs, + discretize_kwargs=discretize_kwargs, + ).to_sky(wcs) + + return PolygonSkyRegion( + self.vertices, + meta=self.meta, + visual=self.visual + ) + + def to_pixel( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None, + ): + + if discretize_kwargs is None: + discretize_kwargs = {} + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + verts = wcs.world_to_pixel( + self.discretize_boundary(**discretize_kwargs).vertices + ) + + return PolygonPixelRegion( + PixCoord(*verts), meta=self.meta.copy(), + visual=self.visual.copy() + ) + + return self.to_sky().to_pixel(wcs) diff --git a/regions/shapes/range.py b/regions/shapes/range.py new file mode 100644 index 00000000..405248b9 --- /dev/null +++ b/regions/shapes/range.py @@ -0,0 +1,893 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +This module defines a range spherical sky region (bounded by longitude +and/or latitude ranges). +""" + +import copy +import operator + +import astropy.units as u +import numpy as np +from astropy.coordinates import Angle, SkyCoord +from astropy.stats import circmean + +from regions._utils.spherical_helpers import ( + cross_product_skycoord2skycoord, cross_product_sum_skycoord2skycoord, + discretize_all_edge_boundaries, get_edge_raw_lonlat_bounds_circ_edges) +from regions.core.attributes import TwoValAngleorNone +from regions.core.compound import CompoundSphericalSkyRegion +from regions.core.core import ComplexSphericalSkyRegion +from regions.core.metadata import RegionMeta, RegionVisual +from regions.core.pixcoord import PixCoord +from regions.shapes.annulus import CircleAnnulusSphericalSkyRegion +from regions.shapes.circle import CircleSphericalSkyRegion +from regions.shapes.lune import LuneSphericalSkyRegion +from regions.shapes.polygon import (PolygonPixelRegion, + PolygonSphericalSkyRegion) +from regions.shapes.whole_sky import WholeSphericalSkyRegion + +__all__ = ['RangeSphericalSkyRegion'] + + +class RangeSphericalSkyRegion(ComplexSphericalSkyRegion): + """ + Sky Region defined within a range of longitude and/or latitudes. At + least one set of longitude or latitude bounds must be set. + + If latitude_range[0] > latitude_range[1], will wrap over the poles + (eg, will *exclude* the central latitude range.) + + Parameters + ---------- + frame : `~astropy.coordinates.BaseCoordinateFrame` + The coordinate frame for the specified range region. + + longitude_range : length 2 list-like of `~astropy.coordinates.Angle` or + `~astropy.units.Quantity` or None + Longitude range in region. Spans from first to second entries + (in increasing longitude value), so inverting the order will select + the complement longitude range. + Longitude values must be within [0, 360] degrees (or equivalent). + + latitude_range : length 2 list-like of `~astropy.coordinates.Angle` or + `~astropy.units.Quantity` or None + Latitude range in region. Spans from first to second entries + (in increasing latitude value), so inverting the order will select + the complement latitude range. + I.e., if latitude_range[0] > latitude_range[1], will wrap over the poles, + and will instead exclude the central latitude range. + Longitude values must be within [-90, 90] degrees (or equivalent). + + meta : `~regions.RegionMeta` or `dict`, optional + A dictionary that stores the meta attributes of the region. + + visual : `~regions.RegionVisual` or `dict`, optional + A dictionary that stores the visual meta attributes of the + region. + + **kwargs : Keyword arguments + Additional keyword arguments passed when created new instances + after coordinate transformations. + """ + + _params = ('frame', 'longitude_range', 'latitude_range') + _boundaries = ('longitude_bounds', 'latitude_bounds') + longitude_range = TwoValAngleorNone( + 'The range of longitude values, as a length-2 |Quantity| angle or list or as None' + ) + latitude_range = TwoValAngleorNone( + 'The range of longitude values, as a length-2 |Quantity| angle or list or as None' + ) + + def __init__( + self, + frame=None, + longitude_range=None, + latitude_range=None, + meta=None, + visual=None, + **kwargs + ): + # Validate inputs: + if ((longitude_range is None) & (latitude_range is None) + & (kwargs.get('_longitude_bounds', None) is None) + & (kwargs.get('_latitude_bounds', None) is None)): + # Non-frame transformation instantiation is only intended to be used + # by setting longitude_range and latitude_range, so note those in the error message + raise ValueError('A range for at least one of longitude and latitude must be set') + + self.longitude_range = longitude_range + self.latitude_range = latitude_range + + # Validate lon/lat range inputs, to ensure ranges have expected definitions. + # Run this after setting the attribute, which performs initial validation using + # TwoValAngleorNone attribute class + self._validate_longitude_range(self.longitude_range) + self._validate_latitude_range(self.latitude_range) + + self._longitude_bounds = None + self._latitude_bounds = None + self._vertices = None + self._is_original_frame = kwargs.pop('_is_original_frame', True) + + # Check if "_frame" is passed as a kwarg, from a copy: + _frame = None + if kwargs.get('_frame', 'missing') != 'missing': + _frame = kwargs.pop('_frame', None) + + # Passing frame overrides any value passed in "_frame" in the kwargs, + # so only if "frame"=None is "_frame" used: + if frame is None: + frame = _frame + + # Validate frame and get into standard format: + frame = self._validate_frame(frame) + + # -------------------- + if ( + ((longitude_range is None) + & (latitude_range is None)) + & ((kwargs.get('_longitude_bounds', None) is not None) + | (kwargs.get('_latitude_bounds', None) is not None)) + ): + # Create class by directly passing long/lat bounds + # -- for cases when creating transformed RangeSphericalSkyRegion instances + + self._longitude_bounds = kwargs.pop('_longitude_bounds', None) + self._latitude_bounds = kwargs.pop('_latitude_bounds', None) + + # Also directly input _vertices, as these are only computed in the + # original frame + self._vertices = kwargs.pop('_vertices', None) + + # Also set _params to to just ("frame") for this transformed instance. + self._params = ('frame',) + + # Validate direct bounds, vertices inputs: + self._validate_inputs_lonlat_bounds_vertices() + + # Stash frame in private attribute to access when building on-the-fly + # derived attributes & bounds: + self._frame = frame + + self.meta = meta or RegionMeta() + self.visual = visual or RegionVisual() + + def __eq__(self, other): + # Equality operator for RangeSphericalSkyRegion. + + # Modified from default, because the lon/lat ranges can be none + # -- requiring comparison of boundaries. + + # All Region properties are compared for strict equality except + # for Quantity parameters, which allow for different units if they + # are directly convertible. + + if not isinstance(other, self.__class__): + return False + + meta_params = ['meta', 'visual'] + intern_params = ['_frame', '_vertices', '_is_original_frame', '_params'] + params_list = ( + list(self._params) + [f"_{bn}" for bn in list(self._boundaries)] + + intern_params + meta_params + ) + + # now check the parameter values + # Note that Quantity comparisons allow for different units + # if they directly convertible (e.g., 1. * u.deg == 60. * u.arcmin) + for param in params_list: + # np.any is used for SkyCoord array comparisons + if np.any(getattr(self, param) != getattr(other, param)): + return False + # Note frame comparisons are done directly, so no + # catching of TypeError is necessary + + return True + + def _validate_longitude_range(self, longitude_range): + # Specifying longitude_range as an TwoValAngleorNone attribute already validates + # for input as either length-2 list-like with angular units, or None + + # Here, only check that if not None, + # the values are within [0,360] deg or equivalent. + if longitude_range is not None: + for i in range(2): + if not ((longitude_range[i].to(u.deg).value >= 0) + & (longitude_range[i].to(u.deg).value <= 360)): + raise ValueError('Longitude values must be within [0, 360] degrees or ' + 'equivalent!') + + def _validate_latitude_range(self, latitude_range): + # Specifying latitude_range as an TwoValAngleorNone attribute already validates + # for input as either length-2 list-like with angular units, or None + + # Here, only check that if not None, + # the values are within [-90,90] deg or equivalent. + if latitude_range is not None: + for i in range(2): + if not ((latitude_range[i].to(u.deg).value >= -90) + & (latitude_range[i].to(u.deg).value <= 90)): + raise ValueError('Latitude values must be within [-90, 90] degrees or ' + 'equivalent!') + + def _validate_inputs_lonlat_bounds_vertices(self): + # Validate direct bounds, vertices inputs: + + # Lon bounds can be lune or None + if not (isinstance(self._longitude_bounds, LuneSphericalSkyRegion) + | (self._longitude_bounds is None)): + raise ValueError('Invalid direct longitude bounds input! Must be ' + 'LuneSphericalSkyRegion or None.') + + # Lat bounds can be circle, annulus, compound xor of whole sky and annulus, or None + if not (isinstance(self._latitude_bounds, (CircleAnnulusSphericalSkyRegion, + CircleSphericalSkyRegion)) + | (isinstance(self._latitude_bounds, CompoundSphericalSkyRegion) + & isinstance(getattr(self._latitude_bounds, 'region1', None), + WholeSphericalSkyRegion) + & isinstance(getattr(self._latitude_bounds, 'region2', None), + CircleAnnulusSphericalSkyRegion) + & (getattr(self._latitude_bounds, 'operator', None) == operator.xor)) + | (self._latitude_bounds is None)): + raise ValueError('Invalid direct latitude bounds input! Must be one of: ' + 'CircleSphericalSkyRegion, CircleAnnulusSphericalSkyRegion, ' + 'CompoundSphericalSkyRegion (as WholeSphericalSkyRegion ' + '^ CircleAnnulusSphericalSkyRegion), or None.') + + # Vertices can be Skycoord of specified length or None, depending on bounds: + if self._bound_nverts > 0: + if (isinstance(self._vertices, SkyCoord) + & (len(self._vertices) != self._bound_nverts)): + raise ValueError('Invalid vertices direct input! Inconsistent with bounds.') + elif (self._vertices is not None) & (self._bound_nverts == 0): + raise ValueError('Invalid vertices direct input! Inconsistent with bounds.') + + # ------------------------------------------------------------------------------- + # ALWAYS derive boundaries on the fly, IF all _params are not None + # --- only a concern for RangeSphericalSkyRegion.... + + def _derive_longitude_bounds(self): + # Internal, on-the-fly boundaries determination. + # A concern if range values change after original _longitude_bounds + # are computed, if these were "static" for all cases + if self.longitude_range is None: + return None + + # Define GC centers based on longitude range + c1 = SkyCoord(self.longitude_range[0], 0 * u.deg, frame=self.frame) + c2 = SkyCoord(self.longitude_range[1], 0 * u.deg, frame=self.frame) + + # Cross product to get pole. + # Invert order if lon_range[1] > lon_range[0] + if self.longitude_range[0] > self.longitude_range[1]: + c_pole = cross_product_skycoord2skycoord(c2, c1) + else: + c_pole = cross_product_skycoord2skycoord(c1, c2) + + # GC centers from cross products: + center_gc1 = cross_product_skycoord2skycoord(c_pole, c1) + center_gc2 = cross_product_skycoord2skycoord(c2, c_pole) + + c_gc1_sph = center_gc1.represent_as('spherical') + c_gc2_sph = center_gc2.represent_as('spherical') + + # Squash small machine-rounding problems in cross product: + # latitude should be 0 for these longitude-describing GCs, + # as they pass through the poles. + # Only possible in the original frame, when defining from longitude_range + + center_gc1 = SkyCoord(c_gc1_sph.lon, 0 * u.deg, frame=self.frame) + center_gc2 = SkyCoord(c_gc2_sph.lon, 0 * u.deg, frame=self.frame) + + return LuneSphericalSkyRegion( + center_gc1, center_gc2, self.meta, self.visual + ) + + def _derive_latitute_bounds(self): + # Internal, on-the-fly boundaries determination. + # A concern if range values change after original _latitude_bounds + # are computed, if these were "static" for all cases + if self.latitude_range is None: + return None + + # Define north pole: 0 lon, 90 lat + c_pole = SkyCoord(0 * u.deg, 90 * u.deg, frame=self.frame) + # Define south pole: + c_pole_s = SkyCoord(0 * u.deg, -90 * u.deg, frame=self.frame) + + lat_arr = self.latitude_range.copy() + + # Logic if this is a "across-polar" region is handled by swapping + # the logic of the latitude boundary itself: + _wraps_poles = False + if lat_arr[0] > lat_arr[1]: + _wraps_poles = True + lat_arr = lat_arr[::-1] + + inner_rad = 90.0 * u.deg - lat_arr[1] + outer_rad = 90.0 * u.deg - lat_arr[0] + + # Special handling: check if either edge of the range is -90 or 90: + if ( + (lat_arr[1] == Angle(90 * u.deg)) + & (lat_arr[0] == Angle(-90 * u.deg)) + ): + # Whole sky: no constraint: + return None + + elif lat_arr[1] == Angle(90 * u.deg): + if _wraps_poles: + # Check if it wrapped pole, with input [90*u.deg, X*u.deg]: + # equivalent to [-90*u.deg, X*u.deg], and is a + # simple circle centered at S pole: + lat_bounds = CircleSphericalSkyRegion( + c_pole_s, lat_arr[0] + 90.0 * u.deg, + self.meta, self.visual + ) + else: + # Simple circle case: + lat_bounds = CircleSphericalSkyRegion( + c_pole, outer_rad, self.meta, self.visual + ) + return lat_bounds + elif lat_arr[0] == Angle(-90 * u.deg): + if _wraps_poles: + # Check if it wrapped pole, with input [X*u.deg, -90*u.deg], + # equivalent to [X*u.deg, 90*u.deg], and is a simple circle + # Radius from input X*u.deg constraint is in inner_rad after swap + lat_bounds = CircleSphericalSkyRegion( + c_pole, inner_rad, self.meta, self.visual + ) + else: + # Simple circle centered at S pole: + lat_bounds = CircleSphericalSkyRegion( + c_pole_s, lat_arr[1] + 90.0 * u.deg, + self.meta, self.visual + ) + return lat_bounds + + # True annular range: + if _wraps_poles: + # Negate by passing as xor with whole sky + return CompoundSphericalSkyRegion( + WholeSphericalSkyRegion(), + CircleAnnulusSphericalSkyRegion( + c_pole, inner_rad, outer_rad, + self.meta, self.visual + ), + operator.xor, self.meta, self.visual + ) + + # Not wrapping poles: + return CircleAnnulusSphericalSkyRegion( + c_pole, inner_rad, outer_rad, self.meta, self.visual + ) + + @property + def _bound_nverts(self): + # Number of vertices for boundary: + # lat-only: 0 (annulus) + # lon-only: 2 (lune) + # both: 4 (range), unless one touches the pole: 3 + + # Only compute params once: + lon_bounds = self.longitude_bounds + lat_bounds = self.latitude_bounds + + if (lon_bounds is not None) & (lat_bounds is not None): + if isinstance(lat_bounds, CompoundSphericalSkyRegion): + nverts = 6 + # pole-wrapping gives 2 disjoint triangles + elif isinstance(lat_bounds, CircleAnnulusSphericalSkyRegion): + # 4 vertices: + nverts = 4 + elif isinstance(lat_bounds, CircleSphericalSkyRegion): + nverts = 3 + elif (lon_bounds is not None) & (lat_bounds is None): + # 2 vertices: + nverts = 2 + else: + # Only lat_bounds: 0 vertices + nverts = 0 + return nverts + + @property + def _edge_circs(self): + """ + Get list of the circles defining the range boundaries. + + Only used if vertices=4, otherwise lune/annulus properties used. + """ + bound_nverts = self._bound_nverts + + if bound_nverts == 6: + raise NotImplementedError + + if bound_nverts == 4: + return [ + self.longitude_bounds._circle_1, + self.latitude_bounds._outer_region, + self.longitude_bounds._circle_2, + self.latitude_bounds._inner_region, + ] + if bound_nverts == 3: + return [ + self.longitude_bounds._circle_1, + self.latitude_bounds, + self.longitude_bounds._circle_2, + ] + + @property + def _compound_region(self): + if self.longitude_bounds is None: + return self.latitude_bounds + + if self.latitude_bounds is None: + return self.longitude_bounds + + return CompoundSphericalSkyRegion( + self.longitude_bounds, self.latitude_bounds, + operator.and_, self.meta, self.visual + ) + + @property + def longitude_bounds(self): + """ + Longitude bounds of the range spherical sky region. + + Defined as a spherical lune, or None if no longitude bound set. + """ + # If set, get bounds from internal attribute: + if getattr(self, '_longitude_bounds', None) is not None: + return self._longitude_bounds + + # Otherwise, check whether range is set. + # If set, derive bounds on-the-fly: + elif self.longitude_range is not None: + return self._derive_longitude_bounds() + + # If both are set to None, then the bounds are intended to be not set: + else: + return None + + @property + def latitude_bounds(self): + """ + Latitude bounds of the range spherical sky region. + + Defined as a spherical circular annulus, a circle (if ending at + a pole), or None if no latitude bound set. + """ + # If set, get bounds from internal attribute: + if getattr(self, '_latitude_bounds', None) is not None: + return self._latitude_bounds + + # Otherwise, check whether range is set. + # If set, derive bounds on-the-fly: + elif self.latitude_range is not None: + return self._derive_latitute_bounds() + + # If both are set to None, then the bounds are intended to be not set: + else: + return None + + @property + def centroid(self): + """ + Region centroid. + + If both longitude and latitude bounds are set: + + Defined as the point equidistant from all vertices. + + However, if this point falls outside of the region + (as in cases very "long and narrow" ranges), + the centroid is instead approximated as the average + of the longitude and latitudes of all vertices. + + If only longitude or only latitude bounds are set, + the centroid is taken as the centroid of that + boundary shape (i.e., the lune centroid or annulus center). + """ + return self.centroid_mindist + + @property + def centroid_avg(self): + """ + An approximate "centroid", taking the average of the vertices + longitude / latitudes, and invert as needed to define an + "inside" point for a wide angle longitude range. + + Will behave poorly for regions crossing the poles. + """ + bound_nverts = self._bound_nverts + + if bound_nverts == 6: + raise NotImplementedError + + if (bound_nverts == 4) | (bound_nverts == 3): + # Both lon and lat bounds + verts_sph = self.vertices.represent_as('spherical') + + lons = verts_sph.lon + lats = verts_sph.lat + + lon = circmean(lons) + lat = np.average(lats) + centroid = SkyCoord(lon, lat, frame=self.frame) + if ( + (not self.contains(centroid)) and ( + not self.longitude_bounds.contains(centroid) + ) + ): + # If centroid not contained in range, or in subset longitude bounds: + # modify longitude by adding 180 deg: + centroid = SkyCoord(lon + 180 * u.deg, lat, frame=self.frame) + + return centroid + elif bound_nverts == 2: + # Only lon bounds: just get centroid of lune: + return self.longitude_bounds.centroid + elif bound_nverts == 0: + # Only lat_bounds: just get center of annulus: + return self.latitude_bounds.center + + @property + def centroid_mindist(self): + """ + Region centroid. + + Defined as the point equidistant from all vertices. + """ + bound_nverts = self._bound_nverts + + if bound_nverts == 6: + raise NotImplementedError + + if (bound_nverts == 4) | (bound_nverts == 3): + # If both lon/lat bounds: + # Calculate from sum of cross products of vertices with cartesian representation + # verts are in CW order: + centroid_mindist = cross_product_sum_skycoord2skycoord(self.vertices) + + if ((not self.contains(centroid_mindist)) + and (not self.latitude_bounds.contains(centroid_mindist))): + # If it's not contained, check if it's because it's + # outside the latitude range: if so, instead just use avg centroid: + + # Problem of very elongated bounds at high latitude -- + # end up outside range. If so, fall back to average + return self.centroid_avg + + return centroid_mindist + elif bound_nverts == 2: + # Only lon bounds: just get centroid of lune: + return self.longitude_bounds.centroid + elif bound_nverts == 0: + # Only lat_bounds: just get center of annulus: + return self.latitude_bounds.center + + @property + def vertices(self): + """ + Region vertices. + + The result depends on which ranges are set: + + If both longitude and latitude ranges are set, returns 4 vertices + (the "corners" of the composite range) as a SkyCoord. + + If only a longitude range is set, returns 2 vertices (the + vertices of the longitude spherical lune boundary) as a SkyCoord. + + If only latitude ranges are set, returns None (no vertices + for the spherical circular annulus latitude bounds). + """ + if getattr(self, '_vertices', None) is not None: + return self._vertices + else: + bound_nverts = self._bound_nverts + + if bound_nverts == 6: + raise NotImplementedError + + if bound_nverts == 4: + # Both lon/lat bounds:on the fly construct from the lon/lat range attrib + # Longitude min is on the right, so derive CW from lower right corner. + + # Note this is never called for transformed instances, + # which have vertices info stashed in _vertices. + verts_lon = [ + self.longitude_range[0], + self.longitude_range[1], + self.longitude_range[1], + self.longitude_range[0] + ] + verts_lat = [ + self.latitude_range[0], + self.latitude_range[0], + self.latitude_range[1], + self.latitude_range[1], + ] + return SkyCoord(verts_lon, verts_lat, frame=self.frame) + + elif bound_nverts == 3: + # Both lon/lat bounds:on the fly construct from the lon/lat range attrib + # Longitude min is on the right, so derive CW from lower right corner. + # However, ONE of the vertices is on a pole. + + # Note this is never called for transformed instances, + # which have vertices info stashed in _vertices. + + if np.abs(self.latitude_range[0].to(u.deg).value) == 90: + verts_lon = [ + self.longitude_range[0], + self.longitude_range[1], + self.longitude_range[0] + ] + verts_lat = [ + self.latitude_range[0], + self.latitude_range[1], + self.latitude_range[1], + ] + + elif np.abs(self.latitude_range[1].to(u.deg).value) == 90: + verts_lon = [ + self.longitude_range[0], + self.longitude_range[1], + self.longitude_range[0] + ] + verts_lat = [ + self.latitude_range[0], + self.latitude_range[0], + self.latitude_range[1], + ] + return SkyCoord(verts_lon, verts_lat, frame=self.frame) + + elif bound_nverts == 2: + # 2 vertices: directly return from lune itself. + return self.longitude_bounds.vertices + elif bound_nverts == 0: + # Only lat_bounds: no vertices + return None + + @property + def bounding_circle(self): + bound_nverts = self._bound_nverts + if bound_nverts == 0: + # No vertices: only an annulus: + # Use the annulus bounding circle: + return self.latitude_bounds.bounding_circle + elif bound_nverts == 2: + # Lune: get lune bounding circle + return self.longitude_bounds.bounding_circle + elif bound_nverts == 4: + # Quadrilateral: use max of seps to vertices: + cent = self.centroid_mindist + seps = cent.separation(self.vertices) + return CircleSphericalSkyRegion(center=cent, radius=np.max(seps)) + elif bound_nverts == 6: + raise NotImplementedError + + @property + def bounding_lonlat(self): + bound_nverts = self._bound_nverts + if bound_nverts == 0: + # Only an annulus: Use the annulus bounding lonlat: + return self.latitude_bounds.bounding_lonlat + + if bound_nverts == 2: + # Lune: get lune bounding lonlat + return self.longitude_bounds.bounding_lonlat + + if (bound_nverts == 4) | (bound_nverts == 3): + # Only used for both lon+lat bounds, + # otherwise directly calls lon or lat bounding_lonlat() + lons_arr, lats_arr = get_edge_raw_lonlat_bounds_circ_edges( + self.vertices, self.centroid, self._edge_circs, + ) + + # If original frame: keep lon bounds if touching a pole + # Only validate and change if not the original frame + # (eg, when the region could cap the pole) + if not self._is_original_frame: + lons_arr, lats_arr = self._validate_lonlat_bounds(lons_arr, lats_arr) + + return lons_arr, lats_arr + + if bound_nverts == 6: + raise NotImplementedError + + def contains(self, coord): + if self._is_original_frame: + # Just compute directly from lon/lat ranges: + + # Ensure coordinate is in same frame: + c_repr = coord.represent_as('spherical') + + # Longitude range constraints: + if self.longitude_range is not None: + # Ensure lon range is wrapped at 360, matching internal SkyCoord convention + # To handle possible cases of lon ranges "wrapping around" across + # the standard 360->0 wrap, wrap lon ranges + coord values + # around the first entry angle + wrap_ang = self.longitude_range[0] + coord_lons = c_repr.lon.wrap_at(wrap_ang) + in_range_lon = ( + (coord_lons >= Angle(self.longitude_range[0]).wrap_at(wrap_ang)) + & (coord_lons <= Angle(self.longitude_range[1]).wrap_at(wrap_ang)) + ) + # No constraints + elif coord.isscalar: + in_range_lon = True + else: + in_range_lon = np.ones(coord.shape, dtype=bool) + + # Latitude range constraints: + if self.latitude_range is not None: + coord_lat = c_repr.lat + + # Check if latitude ranges are "simpler": + # Equivalent to no constraints, or simpler circles: + if ( + (self.latitude_range[1] == Angle(90 * u.deg)) + & (self.latitude_range[0] == Angle(-90 * u.deg)) + ): + # Whole sky: + if coord.isscalar: + in_range_lat = True + else: + in_range_lat = np.ones(coord.shape, dtype=bool) + elif ( + self.latitude_range[1] == Angle(90 * u.deg) + ): + # Range includes latitude up to N pole + in_range_lat = coord_lat >= self.latitude_range[0] + elif ( + self.latitude_range[0] == Angle(-90 * u.deg) + ): + # Range includes latitudes down to S pole + in_range_lat = coord_lat <= self.latitude_range[1] + else: + # Actual annulus + in_range_lat = ( + (coord_lat >= self.latitude_range[0]) + & (coord_lat <= self.latitude_range[1]) + ) + # No constraints + elif coord.isscalar: + in_range_lat = True + else: + in_range_lat = np.ones(coord.shape, dtype=bool) + + in_range = in_range_lon & in_range_lat + + if self.meta.get('include', True): + return in_range + else: + return np.logical_not(in_range) + + # If transformed: must use boundary compound region logic: + return self._compound_region.contains(coord) + + def transform_to(self, frame, merge_attributes=True): + frame = self._validate_frame(frame) + + transformed = {} + + # Internal attribute noting this is a transformed frame. + transformed['_is_original_frame'] = False + + fields = list(self._params) + for field in fields: + # High-level lon/lat attributes will not transform: set to None + transformed[field] = None + + boundaries = list(self._boundaries) + for bound in boundaries: + # Will be passing directly into internal attributes + # _longitude_bounds, _latitude_bounds + bound_orig = getattr(self, bound) + if bound_orig is not None: + transformed[f"_{bound}"] = bound_orig.transform_to( + frame, merge_attributes=merge_attributes + ) + else: + transformed[f"_{bound}"] = getattr(self, bound) + + # Also transform vertices: + for field in ['_vertices']: + verts = copy.deepcopy(getattr(self, field)) + if (verts is None) and self._is_original_frame: + # If internal attribute not set, and original frame, + # get the on-the-fly-derived verts. + verts = self.vertices + if verts is not None: + transformed[field] = verts.transform_to( + frame, merge_attributes=merge_attributes + ) + else: + transformed[field] = None + + # Add transformed frame to transformed dict: + transformed['frame'] = frame + + return self.__class__(**transformed) + + def discretize_boundary(self, n_points=10): + bound_nverts = self._bound_nverts + + if bound_nverts == 0: + return self.latitude_bounds.discretize_boundary(n_points=n_points) + elif bound_nverts == 2: + return self.longitude_bounds.discretize_boundary(n_points=n_points) + elif (bound_nverts == 4) | (bound_nverts == 3): + bound_verts = discretize_all_edge_boundaries( + self.vertices, self._edge_circs, n_points + ) + return PolygonSphericalSkyRegion(bound_verts) + elif bound_nverts == 6: + raise NotImplementedError + + def to_sky( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None + ): + if not include_boundary_distortions: + raise ValueError( + 'Invalid parameter: `include_boundary_distortions=False`!\n' + 'Transforming range to planar sky region is only possible by ' + 'including boundary distortions.' + ) + if discretize_kwargs is None: + discretize_kwargs = {} + + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + # Use to_pixel(), then apply "small angle approx" to get planar sky. + return self.to_pixel( + include_boundary_distortions=include_boundary_distortions, + wcs=wcs, + discretize_kwargs=discretize_kwargs, + ).to_sky(wcs) + + def to_pixel( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None, + ): + if not include_boundary_distortions: + raise ValueError( + 'Invalid parameter: `include_boundary_distortions=False`!\n' + 'Transforming range to planar pixel region is only possible by ' + 'including boundary distortions.' + ) + + if discretize_kwargs is None: + discretize_kwargs = {} + if include_boundary_distortions: + if wcs is None: + raise ValueError( + "'wcs' must be set if 'include_boundary_distortions'=True" + ) + # Requires spherical to planar projection (from WCS) and discretization + disc_bound = self.discretize_boundary(**discretize_kwargs) + # Anticipating complex, wrapped over the poles case: + if isinstance(disc_bound, CompoundSphericalSkyRegion): + return disc_bound.to_pixel(wcs) + + verts = wcs.world_to_pixel(disc_bound.vertices) + + return PolygonPixelRegion( + PixCoord(*verts), meta=self.meta.copy(), + visual=self.visual.copy() + ) diff --git a/regions/shapes/rectangle.py b/regions/shapes/rectangle.py index 773d9711..fab0e696 100644 --- a/regions/shapes/rectangle.py +++ b/regions/shapes/rectangle.py @@ -120,6 +120,10 @@ def to_sky(self, wcs): meta=self.meta.copy(), visual=self.visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + @property def bounding_box(self): w2 = self.width / 2. @@ -418,3 +422,7 @@ def to_pixel(self, wcs): return RectanglePixelRegion(center, width, height, angle=angle, meta=self.meta.copy(), visual=self.visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError diff --git a/regions/shapes/tests/test_annulus.py b/regions/shapes/tests/test_annulus.py index 3800d1ca..630869d1 100644 --- a/regions/shapes/tests/test_annulus.py +++ b/regions/shapes/tests/test_annulus.py @@ -2,22 +2,37 @@ import astropy.units as u import numpy as np -from astropy.coordinates import SkyCoord +import pytest +from astropy.coordinates import Latitude, Longitude, SkyCoord +from astropy.io import fits from astropy.tests.helper import assert_quantity_allclose +from astropy.utils.data import get_pkg_data_filename +from astropy.wcs import WCS from numpy.testing import assert_allclose from regions.core import PixCoord, RegionMeta, RegionVisual +from regions.core.compound import CompoundPixelRegion, CompoundSkyRegion from regions.shapes.annulus import (CircleAnnulusPixelRegion, CircleAnnulusSkyRegion, + CircleAnnulusSphericalSkyRegion, EllipseAnnulusPixelRegion, EllipseAnnulusSkyRegion, RectangleAnnulusPixelRegion, RectangleAnnulusSkyRegion) +from regions.shapes.circle import CircleSphericalSkyRegion from regions.shapes.tests.test_common import (BaseTestPixelRegion, - BaseTestSkyRegion) + BaseTestSkyRegion, + BaseTestSphericalSkyRegion) from regions.tests.helpers import make_simple_wcs +@pytest.fixture(scope='session', name='wcs') +def wcs_fixture(): + filename = get_pkg_data_filename('data/example_header.fits') + header = fits.getheader(filename) + return WCS(header) + + class TestCircleAnnulusPixelRegion(BaseTestPixelRegion): meta = RegionMeta({'text': 'test'}) visual = RegionVisual({'color': 'blue'}) @@ -69,6 +84,21 @@ def test_transformation(self): skyannulus = self.reg.to_sky(wcs=self.wcs) assert isinstance(skyannulus, CircleAnnulusSkyRegion) + def test_to_spherical_sky(self, wcs): + sphskyann = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=False) + assert isinstance(sphskyann, CircleAnnulusSphericalSkyRegion) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=True) + + def test_to_spherical_sky_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_spherical_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + def test_rotate(self): reg = self.reg.rotate(PixCoord(2, 3), 90 * u.deg) assert_allclose(reg.center.xy, (1, 4)) @@ -118,6 +148,21 @@ def test_transformation(self): pixannulus = self.reg.to_pixel(wcs=self.wcs) assert isinstance(pixannulus, CircleAnnulusPixelRegion) + def test_to_spherical_sky(self, wcs): + sphskyann = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=False) + assert isinstance(sphskyann, CircleAnnulusSphericalSkyRegion) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=True) + + def test_to_spherical_sky_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_spherical_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + def test_eq(self): reg = self.reg.copy() assert reg == self.reg @@ -125,6 +170,135 @@ def test_eq(self): assert reg != self.reg +class TestCircleAnnulusSphericalSkyRegion(BaseTestSphericalSkyRegion): + inside = [(3 * u.deg, 4.006 * u.deg)] + outside = [(3 * u.deg, 7 * u.deg)] + meta = RegionMeta({'text': 'test'}) + visual = RegionVisual({'color': 'blue'}) + reg = CircleAnnulusSphericalSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), + 20 * u.arcsec, 30 * u.arcsec, meta=meta, + visual=visual) + skycoord = SkyCoord(3 * u.deg, 4 * u.deg, frame='icrs') + wcs = make_simple_wcs(skycoord, 5 * u.arcsec, 20) + + expected_repr = (', inner_radius=20.0 ' + 'arcsec, outer_radius=30.0 arcsec)>') + expected_str = ('Region: CircleAnnulusSphericalSkyRegion\ncenter: \ninner_radius: ' + '20.0 arcsec\nouter_radius: 30.0 arcsec') + + def test_init(self): + assert_quantity_allclose(self.reg.center.ra, self.skycoord.ra) + assert_quantity_allclose(self.reg.inner_radius, 20 * u.arcsec) + assert_quantity_allclose(self.reg.outer_radius, 30 * u.arcsec) + + def test_valid_radii(self): + with pytest.raises(ValueError) as excinfo: + _ = CircleAnnulusSphericalSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), + 30 * u.arcsec, 20 * u.arcsec) + estr = 'outer_radius must be greater than inner_radius' + assert estr in str(excinfo.value) + + def test_copy(self): + reg = self.reg.copy() + assert_allclose(reg.center.ra.deg, 3) + assert_allclose(reg.inner_radius.to_value('arcsec'), 20) + assert_allclose(reg.outer_radius.to_value('arcsec'), 30) + assert reg.meta == self.meta + assert reg.visual == self.visual + + def test_contains(self): + assert not self.reg.contains(self.skycoord) + test_coord = SkyCoord(3 * u.deg, 10 * u.deg, frame='icrs') + assert not self.reg.contains(test_coord) + + def test_transformation(self): + pixannulus = self.reg.to_pixel(wcs=self.wcs) + assert isinstance(pixannulus, CircleAnnulusPixelRegion) + + skyannulus = self.reg.to_sky(wcs=self.wcs) + assert isinstance(skyannulus, CircleAnnulusSkyRegion) + + polysky = self.reg.to_sky(self.wcs, include_boundary_distortions=True) + assert isinstance(polysky, CompoundSkyRegion) + + polypix = self.reg.to_pixel(self.wcs, include_boundary_distortions=True) + assert isinstance(polypix, CompoundPixelRegion) + + def test_transformation_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_pixel(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + def test_frame_transformation(self): + reg2 = self.reg.transform_to('galactic') + assert reg2.center == self.reg.center.transform_to('galactic') + assert_allclose(reg2.inner_radius.to_value('arcsec'), 20) + assert isinstance(reg2, CircleAnnulusSphericalSkyRegion) + assert reg2.frame.name == 'galactic' + + def test_eq(self): + reg = self.reg.copy() + assert reg == self.reg + reg.inner_radius = 10 * u.arcsec + assert reg != self.reg + + def test_bounding_circle(self): + skycoord = SkyCoord(3 * u.deg, 4 * u.deg) + circ = CircleSphericalSkyRegion(skycoord, 30 * u.arcsec, + meta=self.meta, + visual=self.visual) + + bounding_circle = self.reg.bounding_circle + assert bounding_circle == circ + + def test_bounding_lonlat(self): + skycoord = SkyCoord(3 * u.deg, 0 * u.deg) + reg = CircleAnnulusSphericalSkyRegion(skycoord, + 20 * u.arcsec, + 30 * u.arcsec) + bounding_lonlat = reg.bounding_lonlat + + assert_quantity_allclose(bounding_lonlat[0], + Longitude([3. * u.deg - 30 * u.arcsec, + 3. * u.deg + 30 * u.arcsec])) + + assert_quantity_allclose(bounding_lonlat[1], + Latitude([-30 * u.arcsec, + 30 * u.arcsec])) + + skycoord2 = SkyCoord(3 * u.deg, 90 * u.deg) + reg2 = CircleAnnulusSphericalSkyRegion(skycoord2, + 20 * u.arcsec, + 30 * u.arcsec) + bounding_lonlat2 = reg2.bounding_lonlat + + assert bounding_lonlat2[0] is None + + assert_quantity_allclose(bounding_lonlat2[1], + Latitude([90 * u.deg - 30 * u.arcsec, + 90 * u.deg - 20 * u.arcsec])) + + skycoord3 = SkyCoord(3 * u.deg, -90 * u.deg) + reg3 = CircleAnnulusSphericalSkyRegion(skycoord3, + 20 * u.arcsec, + 30 * u.arcsec) + bounding_lonlat3 = reg3.bounding_lonlat + + assert bounding_lonlat3[0] is None + + assert_quantity_allclose(bounding_lonlat3[1], + Latitude([-90 * u.deg + 20 * u.arcsec, + -90 * u.deg + 30 * u.arcsec])) + + class TestEllipseAnnulusPixelRegion(BaseTestPixelRegion): meta = RegionMeta({'text': 'test'}) visual = RegionVisual({'color': 'blue'}) diff --git a/regions/shapes/tests/test_api.py b/regions/shapes/tests/test_api.py index f7c8c4c2..a8d52b18 100644 --- a/regions/shapes/tests/test_api.py +++ b/regions/shapes/tests/test_api.py @@ -7,24 +7,31 @@ import astropy.units as u import pytest -from astropy.coordinates import SkyCoord +from astropy.coordinates import Angle, SkyCoord from astropy.wcs import WCS -from regions.core.core import PixelRegion, Region, SkyRegion +from regions.core.core import (PixelRegion, Region, SkyRegion, + SphericalSkyRegion) from regions.core.mask import RegionMask from regions.core.metadata import RegionMeta, RegionVisual from regions.core.pixcoord import PixCoord from regions.shapes.annulus import (CircleAnnulusPixelRegion, CircleAnnulusSkyRegion, + CircleAnnulusSphericalSkyRegion, EllipseAnnulusPixelRegion, EllipseAnnulusSkyRegion, RectangleAnnulusPixelRegion, RectangleAnnulusSkyRegion) -from regions.shapes.circle import CirclePixelRegion, CircleSkyRegion +from regions.shapes.circle import (CirclePixelRegion, CircleSkyRegion, + CircleSphericalSkyRegion) from regions.shapes.ellipse import EllipsePixelRegion, EllipseSkyRegion +from regions.shapes.lune import LuneSphericalSkyRegion from regions.shapes.point import PointPixelRegion, PointSkyRegion -from regions.shapes.polygon import PolygonPixelRegion, PolygonSkyRegion +from regions.shapes.polygon import (PolygonPixelRegion, PolygonSkyRegion, + PolygonSphericalSkyRegion) +from regions.shapes.range import RangeSphericalSkyRegion from regions.shapes.rectangle import RectanglePixelRegion, RectangleSkyRegion +from regions.shapes.whole_sky import WholeSphericalSkyRegion PIXEL_REGIONS = [ CirclePixelRegion(PixCoord(3, 4), radius=5), @@ -51,6 +58,21 @@ 5 * u.deg, 5 * u.deg, 7 * u.deg), PointSkyRegion(SkyCoord(6 * u.deg, 5 * u.deg))] +SPHERICAL_SKY_REGIONS = [ + CircleSphericalSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), radius=5 * u.deg), + CircleAnnulusSphericalSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), 5 * u.deg, + 7 * u.deg), + LuneSphericalSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), + SkyCoord(153 * u.deg, 4 * u.deg)), + PolygonSphericalSkyRegion(SkyCoord([1, 4, 3] * u.deg, [2, 4, 4] * u.deg)), + RangeSphericalSkyRegion( + longitude_range=[0, 10] * u.deg, + latitude_range=[40, 60] * u.deg, + frame='icrs'), + WholeSphericalSkyRegion()] + +INCLUDE_BOUNDARY_DISTORTIONS = [True, False] + MASK_MODES = ['center', 'exact', 'subpixels'] COMMON_WCS = WCS(naxis=2) COMMON_WCS.wcs.ctype = 'RA---TAN', 'DEC--TAN' @@ -88,6 +110,41 @@ def test_pix_to_sky(region): pytest.xfail() +@pytest.mark.parametrize('region', PIXEL_REGIONS, ids=ids_func) +@pytest.mark.parametrize('include_dist', INCLUDE_BOUNDARY_DISTORTIONS, ids=ids_func) +def test_pix_to_spherical_sky(region, include_dist): + # TODO: remove expected failures when implemented + # Expected failure: + # No spherical Ellipse, EllipseAnnulus, Point, Rectangle, RectangleAnnulus + # Also expected failure: + # Boundary distortions not yet implemented for + # CirclePixelRegion, CircleAnnulusPixelRegion, PolygonPixelRegion + if ( + isinstance(region, + (EllipsePixelRegion, + EllipseAnnulusPixelRegion, + PointPixelRegion, + RectanglePixelRegion, + RectangleAnnulusPixelRegion)) + | ((isinstance(region, + (CirclePixelRegion, + CircleAnnulusPixelRegion, + PolygonPixelRegion))) & include_dist) + ): + with pytest.raises(NotImplementedError): + sph_sky_region = region.to_spherical_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + assert isinstance(sph_sky_region, SphericalSkyRegion) + else: + sph_sky_region = region.to_spherical_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + assert isinstance(sph_sky_region, SphericalSkyRegion) + + @pytest.mark.parametrize('region', PIXEL_REGIONS, ids=ids_func) @pytest.mark.parametrize('mode', MASK_MODES, ids=ids_func) def test_pix_to_mask(region, mode): @@ -115,6 +172,118 @@ def test_sky_to_pix(region): assert isinstance(pix_region, PixelRegion) +@pytest.mark.parametrize('region', SKY_REGIONS, ids=ids_func) +@pytest.mark.parametrize('include_dist', INCLUDE_BOUNDARY_DISTORTIONS, ids=ids_func) +def test_sky_to_spherical_sky(region, include_dist): + # TODO: remove expected failures when implemented + # Expected failure: + # No spherical Ellipse, EllipseAnnulus, Point, Rectangle, RectangleAnnulus + # Also expected failure: + # Boundary distortions not yet implemented for + # CircleSkyRegion, CircleAnnulusSkyRegion, PolygonSkyRegion + if ( + isinstance(region, + (EllipseSkyRegion, + EllipseAnnulusSkyRegion, + PointSkyRegion, + RectangleSkyRegion, + RectangleAnnulusSkyRegion)) + | ((isinstance(region, + (CircleSkyRegion, + CircleAnnulusSkyRegion, + PolygonSkyRegion))) & include_dist) + ): + with pytest.raises(NotImplementedError): + sph_sky_region = region.to_spherical_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + assert isinstance(sph_sky_region, SphericalSkyRegion) + else: + sph_sky_region = region.to_spherical_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + assert isinstance(sph_sky_region, SphericalSkyRegion) + + +@pytest.mark.parametrize('region', SPHERICAL_SKY_REGIONS, ids=ids_func) +@pytest.mark.parametrize('include_dist', INCLUDE_BOUNDARY_DISTORTIONS, ids=ids_func) +def test_spherical_sky_to_sky(region, include_dist): + # TODO: remove excepted failure for Lune+distortion once Lune handling is implemented + if isinstance(region, WholeSphericalSkyRegion): + with pytest.raises(NotImplementedError): + _ = region.to_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + + elif isinstance(region, LuneSphericalSkyRegion) & include_dist: + # Include distortions + with pytest.raises(NotImplementedError): + _ = region.to_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + + elif (isinstance(region, (LuneSphericalSkyRegion, RangeSphericalSkyRegion)) + & (not include_dist)): + # No distortions: not defined: + with pytest.raises(ValueError) as excinfo: + _ = region.to_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + + estr = 'Invalid parameter: `include_boundary_distortions=False`!' + assert estr in str(excinfo.value) + + else: + sky_region = region.to_sky( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + assert isinstance(sky_region, SkyRegion) + + +@pytest.mark.parametrize('region', SPHERICAL_SKY_REGIONS, ids=ids_func) +@pytest.mark.parametrize('include_dist', INCLUDE_BOUNDARY_DISTORTIONS, ids=ids_func) +def test_spherical_sky_to_pix(region, include_dist): + # TODO: remove excepted failure for Lune+distortion once Lune handling is implemented + if isinstance(region, WholeSphericalSkyRegion): + with pytest.raises(NotImplementedError): + _ = region.to_pixel( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + + elif isinstance(region, LuneSphericalSkyRegion) & include_dist: + # Include distortions + with pytest.raises(NotImplementedError): + _ = region.to_pixel( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + + elif (isinstance(region, (LuneSphericalSkyRegion, RangeSphericalSkyRegion)) + & (not include_dist)): + # No distortions: not defined: + with pytest.raises(ValueError) as excinfo: + _ = region.to_pixel( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + estr = 'Invalid parameter: `include_boundary_distortions=False`!' + assert estr in str(excinfo.value) + + else: + pixel_region = region.to_pixel( + COMMON_WCS, + include_boundary_distortions=include_dist + ) + assert isinstance(pixel_region, PixelRegion) + + @pytest.mark.parametrize('region', PIXEL_REGIONS, ids=ids_func) def test_attribute_validation_pixel_regions(region): invalid_values = dict(center=[PixCoord([1, 2], [2, 3]), 1, @@ -182,6 +351,49 @@ def test_attribute_validation_sky_regions(region): assert f'{attr!r} must' in str(excinfo.value) +@pytest.mark.parametrize('region', SPHERICAL_SKY_REGIONS, ids=ids_func) +def test_attribute_validation_spherical_sky_regions(region): + invalid_values = dict(center=[PixCoord([1, 2], [2, 3]), 1, + SkyCoord([1 * u.deg], [2 * u.deg]), + (10, 10), (10 * u.deg, 10 * u.deg)], + radius=[u.Quantity([1 * u.deg, 5 * u.deg]), + [1], SkyCoord(1 * u.deg, 2 * u.deg), + 1, 3 * u.km, 0.0 * u.deg, -10. * u.deg], + angle=[u.Quantity([1 * u.deg, 2 * u.deg]), 2, + SkyCoord(1 * u.deg, 2 * u.deg), 3. * u.km], + vertices=[u.Quantity('1deg'), 2, + SkyCoord(1 * u.deg, 2 * u.deg), + SkyCoord([[1 * u.deg, 2 * u.deg]], + [[2 * u.deg, 3 * u.deg]]), + 3 * u.km, (10, 10), + (10 * u.deg, 10 * u.deg)], + longitude_range=[1, [1], [1, 2], [1 * u.deg], + [3 * u.km, 5 * u.km], u.Quantity([1 * u.deg]), + u.Quantity([3 * u.km, 5 * u.km]), + Angle([1 * u.deg])]) + + invalid_values['width'] = invalid_values['radius'] + invalid_values['height'] = invalid_values['radius'] + invalid_values['inner_height'] = invalid_values['radius'] + invalid_values['inner_width'] = invalid_values['radius'] + invalid_values['outer_height'] = invalid_values['radius'] + invalid_values['outer_width'] = invalid_values['radius'] + invalid_values['inner_radius'] = invalid_values['radius'] + invalid_values['outer_radius'] = invalid_values['radius'] + invalid_values['start'] = invalid_values['center'] + invalid_values['end'] = invalid_values['radius'] + invalid_values['latitude_range'] = invalid_values['longitude_range'] + + region_attrs = [attr for attr, _ in vars(region).items()] + + for attr in invalid_values: + if hasattr(region, attr) & (attr in region_attrs): + for val in invalid_values.get(attr, None): + with pytest.raises(ValueError) as excinfo: + setattr(region, attr, val) + assert f'{attr!r} must' in str(excinfo.value) + + @pytest.mark.parametrize('region', PIXEL_REGIONS + SKY_REGIONS, ids=ids_func) def test_metadata(region): region.meta = {'text': 'hello'} diff --git a/regions/shapes/tests/test_circle.py b/regions/shapes/tests/test_circle.py index 3010f188..5a5d921b 100644 --- a/regions/shapes/tests/test_circle.py +++ b/regions/shapes/tests/test_circle.py @@ -3,7 +3,7 @@ import astropy.units as u import numpy as np import pytest -from astropy.coordinates import SkyCoord +from astropy.coordinates import Latitude, Longitude, SkyCoord from astropy.io import fits from astropy.tests.helper import assert_quantity_allclose from astropy.utils.data import get_pkg_data_filename @@ -12,9 +12,12 @@ from regions._utils.optional_deps import HAS_MATPLOTLIB from regions.core import PixCoord, RegionMeta, RegionVisual -from regions.shapes.circle import CirclePixelRegion, CircleSkyRegion +from regions.shapes.circle import (CirclePixelRegion, CircleSkyRegion, + CircleSphericalSkyRegion) +from regions.shapes.polygon import PolygonPixelRegion, PolygonSkyRegion from regions.shapes.tests.test_common import (BaseTestPixelRegion, - BaseTestSkyRegion) + BaseTestSkyRegion, + BaseTestSphericalSkyRegion) from regions.tests.helpers import make_simple_wcs @@ -59,6 +62,21 @@ def test_pix_sky_roundtrip(self): assert reg_new.meta['text'] != self.reg.meta['text'] assert reg_new.visual['color'] != self.reg.visual['color'] + def test_to_spherical_sky(self, wcs): + sphskycircle = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=False) + assert isinstance(sphskycircle, CircleSphericalSkyRegion) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=True) + + def test_to_spherical_sky_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_spherical_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + @pytest.mark.skipif(not HAS_MATPLOTLIB, reason='matplotlib is required') def test_as_artist(self): patch = self.reg.as_artist() @@ -117,6 +135,20 @@ def test_transformation(self, wcs): skycircle2.center.data.lat) assert_quantity_allclose(skycircle2.radius, skycircle.radius) + sphskycircle = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=False) + assert isinstance(sphskycircle, CircleSphericalSkyRegion) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=True) + + def test_to_spherical_sky_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_spherical_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + def test_dimension_center(self): center = SkyCoord([1, 2] * u.deg, [3, 4] * u.deg) radius = 2 * u.arcsec @@ -140,3 +172,111 @@ def test_eq(self): def test_zero_size(self): with pytest.raises(ValueError): CircleSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), 0. * u.arcsec) + + +class TestCircleSphericalSkyRegion(BaseTestSphericalSkyRegion): + inside = [(3 * u.deg, 4 * u.deg)] + outside = [(3 * u.deg, 0 * u.deg)] + meta = RegionMeta({'text': 'test'}) + visual = RegionVisual({'color': 'blue'}) + reg = CircleSphericalSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), 2 * u.arcsec, + meta=meta, visual=visual) + + expected_repr = (', radius=2.0 arcsec)>') + expected_str = ('Region: CircleSphericalSkyRegion\ncenter: \nradius: 2.0 ' + 'arcsec') + + def test_copy(self): + reg = self.reg.copy() + assert_allclose(reg.center.ra.deg, 3) + assert_allclose(reg.radius.to_value('arcsec'), 2) + assert reg.meta == self.meta + assert reg.visual == self.visual + + def test_transformation(self, wcs): + skycoord = SkyCoord(3 * u.deg, 4 * u.deg, frame='galactic') + sphskycircle = CircleSphericalSkyRegion(skycoord, 2 * u.arcsec) + + pixcircle = sphskycircle.to_pixel(wcs) + + assert_allclose(pixcircle.center.x, -50.5) + assert_allclose(pixcircle.center.y, 299.5) + assert_allclose(pixcircle.radius, 0.027777777777828305) + + skycircle = sphskycircle.to_sky(wcs) + assert isinstance(skycircle, CircleSkyRegion) + + sphskycircle2 = pixcircle.to_spherical_sky(wcs) + + assert_quantity_allclose(sphskycircle.center.data.lon, + sphskycircle2.center.data.lon) + assert_quantity_allclose(sphskycircle.center.data.lat, + sphskycircle2.center.data.lat) + assert_quantity_allclose(sphskycircle2.radius, sphskycircle.radius) + + polysky = sphskycircle.to_sky(wcs, include_boundary_distortions=True) + assert isinstance(polysky, PolygonSkyRegion) + + polypix = sphskycircle.to_pixel(wcs, include_boundary_distortions=True) + assert isinstance(polypix, PolygonPixelRegion) + + def test_transformation_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_pixel(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + def test_frame_transformation(self): + skycoord = SkyCoord(3 * u.deg, 4 * u.deg, frame='galactic') + reg = CircleSphericalSkyRegion(skycoord, 2 * u.arcsec) + + reg2 = reg.transform_to('icrs') + assert reg2.center == skycoord.transform_to('icrs') + assert_allclose(reg2.radius.to_value('arcsec'), 2) + assert isinstance(reg2, CircleSphericalSkyRegion) + assert reg2.frame.name == 'icrs' + + def test_dimension_center(self): + center = SkyCoord([1, 2] * u.deg, [3, 4] * u.deg) + radius = 2 * u.arcsec + with pytest.raises(ValueError) as excinfo: + CircleSphericalSkyRegion(center, radius) + estr = "'center' must be a scalar SkyCoord" + assert estr in str(excinfo.value) + + def test_eq(self): + reg = self.reg.copy() + assert reg == self.reg + reg.radius = 3 * u.arcsec + assert reg != self.reg + + def test_zero_size(self): + with pytest.raises(ValueError): + CircleSkyRegion(SkyCoord(3 * u.deg, 4 * u.deg), 0. * u.arcsec) + + def test_bounding_circle(self): + skycoord = SkyCoord(3 * u.deg, 4 * u.deg, frame='galactic') + reg = CircleSphericalSkyRegion(skycoord, 2 * u.arcsec) + + bounding_circle = reg.bounding_circle + assert bounding_circle == reg + + def test_bounding_lonlat(self): + skycoord = SkyCoord(3 * u.deg, 0 * u.deg, frame='galactic') + reg = CircleSphericalSkyRegion(skycoord, 2 * u.arcsec) + bounding_lonlat = reg.bounding_lonlat + + assert_quantity_allclose(bounding_lonlat[0], + Longitude([3. * u.deg - 2 * u.arcsec, + 3. * u.deg + 2 * u.arcsec])) + + assert_quantity_allclose(bounding_lonlat[1], + Latitude([-2 * u.arcsec, + 2 * u.arcsec])) diff --git a/regions/shapes/tests/test_common.py b/regions/shapes/tests/test_common.py index 74e6435f..5f4d7713 100644 --- a/regions/shapes/tests/test_common.py +++ b/regions/shapes/tests/test_common.py @@ -5,6 +5,7 @@ import numpy as np import pytest +from astropy.coordinates import SkyCoord from numpy.testing import assert_allclose, assert_equal from regions._utils.optional_deps import HAS_MATPLOTLIB @@ -83,3 +84,22 @@ def test_plot(self): class BaseTestSkyRegion(BaseTestRegion): # TODO: here we should add inside/outside tests as above pass + + +class BaseTestSphericalSkyRegion(BaseTestRegion): + def test_contains_scalar(self): + if len(self.inside) > 0: + skycoord = SkyCoord(*self.inside[0]) + assert self.reg.contains(skycoord) + + if len(self.outside) > 0: + skycoord = SkyCoord(*self.outside[0]) + assert not self.reg.contains(skycoord) + + def test_contains_array_1d(self): + coos = self.inside + self.outside + lon, lat = zip(*coos) + skycoord = SkyCoord(list(lon), list(lat)) + actual = self.reg.contains(skycoord) + assert_equal(actual[:len(self.inside)], True) + assert_equal(actual[len(self.inside):], False) diff --git a/regions/shapes/tests/test_lune.py b/regions/shapes/tests/test_lune.py new file mode 100644 index 00000000..9be99b51 --- /dev/null +++ b/regions/shapes/tests/test_lune.py @@ -0,0 +1,123 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import astropy.units as u +import pytest +from astropy.coordinates import Latitude, Longitude, SkyCoord +from astropy.io import fits +from astropy.tests.helper import assert_quantity_allclose +from astropy.utils.data import get_pkg_data_filename +from astropy.wcs import WCS +from numpy.testing import assert_allclose + +from regions.core import RegionMeta, RegionVisual +from regions.shapes.circle import CircleSphericalSkyRegion +from regions.shapes.lune import LuneSphericalSkyRegion +from regions.shapes.polygon import PolygonSphericalSkyRegion +from regions.shapes.tests.test_common import BaseTestSphericalSkyRegion + + +@pytest.fixture(scope='session', name='wcs') +def wcs_fixture(): + filename = get_pkg_data_filename('data/example_header.fits') + header = fits.getheader(filename) + return WCS(header) + + +class TestLuneSphericalSkyRegion(BaseTestSphericalSkyRegion): + inside = [(90 * u.deg, 0 * u.deg)] + outside = [(75 * u.deg, 0 * u.deg)] + meta = RegionMeta({'text': 'test'}) + visual = RegionVisual({'color': 'blue'}) + reg = LuneSphericalSkyRegion(SkyCoord(3 * u.deg, 0 * u.deg), + SkyCoord(178 * u.deg, 0 * u.deg), + meta=meta, visual=visual) + + expected_repr = (', center_gc2=)>') + expected_str = ('Region: LuneSphericalSkyRegion\ncenter_gc1: \ncenter_gc2: ') + + def test_copy(self): + reg = self.reg.copy() + assert_allclose(reg.center_gc1.ra.deg, 3) + assert_allclose(reg.center_gc2.ra.deg, 178) + assert reg.meta == self.meta + assert reg.visual == self.visual + + def test_transformation(self, wcs): + # Test error for no boundary distortions: + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_pixel(wcs) + estr = 'Invalid parameter: `include_boundary_distortions=False`!' + assert estr in str(excinfo.value) + + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_sky(wcs) + estr = 'Invalid parameter: `include_boundary_distortions=False`!' + assert estr in str(excinfo.value) + + # Test for not implemented with include_boundary_distortions=True + with pytest.raises(NotImplementedError): + _ = self.reg.to_sky(wcs, + include_boundary_distortions=True, + discretize_kwargs={'n_points': 2}) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_pixel(wcs, + include_boundary_distortions=True, + discretize_kwargs={'n_points': 2}) + + def test_frame_transformation(self): + reg2 = self.reg.transform_to('galactic') + reg2_cent = reg2.centroid + transf_reg_cent = self.reg.centroid.transform_to('galactic') + assert isinstance(reg2, LuneSphericalSkyRegion) + assert transf_reg_cent.frame.name == reg2_cent.frame.name + assert_quantity_allclose(reg2_cent.l, transf_reg_cent.l) + assert_quantity_allclose(reg2_cent.b, transf_reg_cent.b) + assert_allclose(reg2.center_gc1.l.deg, + self.reg.center_gc1.transform_to('galactic').l.deg) + assert reg2.frame.name == 'galactic' + + def test_eq(self): + reg = self.reg.copy() + assert reg == self.reg + reg.center_gc1 = SkyCoord(4 * u.deg, 0 * u.deg) + assert reg != self.reg + + def test_bounding_circle(self): + skycoord = SkyCoord(90.5 * u.deg, 0 * u.deg) + reg = CircleSphericalSkyRegion(skycoord, 90 * u.deg) + + bc = self.reg.bounding_circle + assert bc == reg + + def test_bounding_lonlat(self): + bounding_lonlat = self.reg.bounding_lonlat + + assert_quantity_allclose(bounding_lonlat[0], + Longitude([88. * u.deg, + 93. * u.deg])) + + assert_quantity_allclose(bounding_lonlat[1], + Latitude([-90 * u.deg, + 90 * u.deg])) + + reg2 = LuneSphericalSkyRegion(SkyCoord(3 * u.deg, 20 * u.deg), + SkyCoord(178 * u.deg, 20 * u.deg), + meta=self.meta, visual=self.visual) + bounding_lonlat2 = reg2.bounding_lonlat + assert bounding_lonlat2[0] is None + + reg3 = LuneSphericalSkyRegion(SkyCoord(3 * u.deg, -20 * u.deg), + SkyCoord(178 * u.deg, -20 * u.deg), + meta=self.meta, visual=self.visual) + bounding_lonlat3 = reg3.bounding_lonlat + assert bounding_lonlat3[0] is None + + def test_discretize_boundary(self): + polylune = self.reg.discretize_boundary(n_points=100) + assert isinstance(polylune, PolygonSphericalSkyRegion) + assert len(polylune.vertices) == 200 diff --git a/regions/shapes/tests/test_polygon.py b/regions/shapes/tests/test_polygon.py index a9f98f06..8f053bfb 100644 --- a/regions/shapes/tests/test_polygon.py +++ b/regions/shapes/tests/test_polygon.py @@ -5,18 +5,21 @@ import astropy.units as u import numpy as np import pytest -from astropy.coordinates import SkyCoord +from astropy.coordinates import Latitude, Longitude, SkyCoord from astropy.tests.helper import assert_quantity_allclose from numpy.testing import assert_allclose, assert_equal from regions._utils.examples import make_example_dataset from regions._utils.optional_deps import HAS_MATPLOTLIB from regions.core import PixCoord, RegionBoundingBox, RegionMeta, RegionVisual +from regions.shapes.circle import CircleSphericalSkyRegion from regions.shapes.polygon import (PolygonPixelRegion, PolygonSkyRegion, + PolygonSphericalSkyRegion, RegularPolygonPixelRegion) from regions.shapes.tests.test_common import (BaseTestPixelRegion, - BaseTestSkyRegion) -from regions.tests.helpers import make_simple_wcs + BaseTestSkyRegion, + BaseTestSphericalSkyRegion) +from regions.tests.helpers import assert_skycoord_allclose, make_simple_wcs @pytest.fixture(scope='session', name='wcs') @@ -70,6 +73,21 @@ def test_pix_sky_roundtrip(self): assert reg_new.meta['text'] != self.reg.meta['text'] assert reg_new.visual['color'] != self.reg.visual['color'] + def test_to_spherical_sky(self, wcs): + polysphsky = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=False) + assert isinstance(polysphsky, PolygonSphericalSkyRegion) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=True) + + def test_to_spherical_sky_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_spherical_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + def test_bounding_box(self): bbox = self.reg.bounding_box assert bbox == RegionBoundingBox(ixmin=1, ixmax=4, iymin=1, iymax=5) @@ -173,6 +191,20 @@ def test_transformation(self, wcs): assert_quantity_allclose(poly.vertices.data.lat, self.reg.vertices.data.lat, atol=1e-3 * u.deg) + polysphsky = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=False) + assert isinstance(polysphsky, PolygonSphericalSkyRegion) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_spherical_sky(wcs, + include_boundary_distortions=True) + + def test_to_spherical_sky_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_spherical_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + def test_contains(self, wcs): position = SkyCoord([1, 3.25] * u.deg, [2, 3.75] * u.deg) # 1,2 is outside, 3.25,3.75 should be inside the triangle... @@ -186,7 +218,7 @@ def test_eq(self): assert reg != self.reg -class TestRegionPolygonPixelRegion(BaseTestPixelRegion): +class TestRegularPolygonPixelRegion(BaseTestPixelRegion): meta = RegionMeta({'text': 'test'}) visual = RegionVisual({'color': 'blue'}) reg = RegularPolygonPixelRegion(PixCoord(50, 50), 8, 20, angle=25 * u.deg, @@ -281,3 +313,118 @@ def test_to_polygon(self): assert polyreg.meta == self.meta assert polyreg.visual == self.visual assert polyreg.origin == PixCoord(0, 0) + + +class TestPolygonSphericalSkyRegion(BaseTestSphericalSkyRegion): + inside = [(3.1 * u.deg, 3.5 * u.deg)] + outside = [(3 * u.deg, 0 * u.deg)] + meta = RegionMeta({'text': 'test'}) + visual = RegionVisual({'color': 'blue'}) + reg = PolygonSphericalSkyRegion(SkyCoord([3, 4, 3] * u.deg, + [3, 4, 4] * u.deg), + meta=meta, visual=visual) + + expected_repr = (')>') + expected_str = ('Region: PolygonSphericalSkyRegion\nvertices: ') + + def test_copy(self): + reg = self.reg.copy() + assert_allclose(reg.vertices.ra.deg, [3, 4, 3]) + assert_allclose(reg.vertices.dec.deg, [3, 4, 4]) + assert reg.meta == self.meta + assert reg.visual == self.visual + + def test_transformation(self, wcs): + polypix = self.reg.to_pixel(wcs) + assert isinstance(polypix, PolygonPixelRegion) + assert_allclose(polypix.vertices.x, + [11.18799228, 10.97633218, 11.02403192]) + assert_allclose(polypix.vertices.y, + [1.99948607, 2.0390009, 2.07707564]) + + polysky = self.reg.to_sky(wcs) + assert isinstance(polysky, PolygonSkyRegion) + assert_allclose(polysky.vertices.ra.deg, [3, 4, 3]) + assert_allclose(polysky.vertices.dec.deg, [3, 4, 4]) + + polysky2 = polysky.to_spherical_sky(wcs) + + assert_quantity_allclose(self.reg.vertices.ra.deg, + polysky2.vertices.ra.deg) + assert_quantity_allclose(self.reg.vertices.dec.deg, + polysky2.vertices.dec.deg) + + polysky3 = self.reg.to_sky(wcs, + include_boundary_distortions=True, + discretize_kwargs={'n_points': 10}) + assert isinstance(polysky3, PolygonSkyRegion) + assert len(polysky3.vertices) == 30 + + polypix2 = self.reg.to_pixel(wcs, + include_boundary_distortions=True, + discretize_kwargs={'n_points': 10}) + assert isinstance(polypix2, PolygonPixelRegion) + assert len(polypix2.vertices) == 30 + + def test_transformation_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_pixel(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + def test_frame_transformation(self): + reg = self.reg.transform_to('galactic') + assert_skycoord_allclose(reg.centroid_mindist, + self.reg.centroid_mindist.transform_to('galactic')) + assert_skycoord_allclose(reg.vertices, + self.reg.vertices.transform_to('galactic')) + assert isinstance(reg, PolygonSphericalSkyRegion) + assert reg.frame.name == 'galactic' + assert reg != self.reg + + def test_eq(self): + reg = self.reg.copy() + assert reg == self.reg + reg.vertices = SkyCoord([3, 5, 3] * u.deg, [3, 4, 4] * u.deg) + assert reg != self.reg + + def test_bounding_circle(self): + skycoord = SkyCoord(3.3333295725289807 * u.deg, + 3.6666666666666665 * u.deg, + frame='icrs') + reg = CircleSphericalSkyRegion(skycoord, + 0.7451014340169484 * u.deg) + + bc = self.reg.bounding_circle + assert bc == reg + + def test_bounding_lonlat(self): + bounding_lonlat = self.reg.bounding_lonlat + + assert_quantity_allclose(bounding_lonlat[0], + Longitude([3 * u.deg, + 4 * u.deg])) + + # GC connecting top two vertices "bluges up" above + # the latitude value = 4 deg for those two points + assert_quantity_allclose(bounding_lonlat[1], + Latitude([3 * u.deg, + 4.00015182 * u.deg])) + + def test_centroid(self): + # Test default polygon has mindist centroid outside region, + # so falls back to average lon/lat centroid definition: + assert self.reg.centroid == self.reg.centroid_avg + + # Test mindist centroid: + reg = PolygonSphericalSkyRegion(SkyCoord([3, 5, 3] * u.deg, + [3, 4, 5] * u.deg)) + + assert reg.centroid == reg.centroid_mindist diff --git a/regions/shapes/tests/test_range.py b/regions/shapes/tests/test_range.py new file mode 100644 index 00000000..86d9ede1 --- /dev/null +++ b/regions/shapes/tests/test_range.py @@ -0,0 +1,506 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import astropy.units as u +import pytest +from astropy.coordinates import Angle, Latitude, Longitude, SkyCoord +from astropy.io import fits +from astropy.tests.helper import assert_quantity_allclose +from astropy.utils.data import get_pkg_data_filename +from astropy.wcs import WCS +from numpy.testing import assert_allclose, assert_equal + +from regions.core import RegionMeta, RegionVisual +from regions.core.compound import CompoundSphericalSkyRegion +from regions.shapes.circle import CircleSphericalSkyRegion +from regions.shapes.lune import LuneSphericalSkyRegion +from regions.shapes.polygon import (PolygonPixelRegion, PolygonSkyRegion, + PolygonSphericalSkyRegion) +from regions.shapes.range import RangeSphericalSkyRegion +from regions.shapes.tests.test_common import BaseTestSphericalSkyRegion + + +@pytest.fixture(scope='session', name='wcs') +def wcs_fixture(): + filename = get_pkg_data_filename('data/example_header.fits') + header = fits.getheader(filename) + return WCS(header) + + +class TestRangeSphericalSkyRegion(BaseTestSphericalSkyRegion): + inside = [(5 * u.deg, 0 * u.deg)] + outside = [(75 * u.deg, 10 * u.deg)] + meta = RegionMeta({'text': 'test'}) + visual = RegionVisual({'color': 'blue'}) + reg = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[-4, 4] * u.deg, + frame='icrs', + meta=meta, visual=visual) + + expected_repr = ('') + expected_str = ('Region: RangeSphericalSkyRegion\n' + 'frame: icrs\n' + 'longitude_range: [ 0. 10.] deg\n' + 'latitude_range: [-4. 4.] deg') + + def test_no_range_set(self): + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(frame='icrs') + estr = 'A range for at least one of longitude and latitude must be set' + assert estr in str(excinfo.value) + + def test_invalid_lon_range(self): + # Test input types: + for lon_range in [u.Quantity([0, 10], u.m / u.s), [0 * u.m / u.s, 30 * u.m / u.s]]: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(longitude_range=lon_range, frame='icrs') + estr = 'must have angular units' + assert estr in str(excinfo.value) + + for lon_range in [u.Quantity([0], u.deg), Angle([2 * u.radian])]: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(longitude_range=lon_range, frame='icrs') + estr = 'must be length 2' + assert estr in str(excinfo.value) + + for lon_range in [[0 * u.deg], [0, 10], 'string']: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(longitude_range=lon_range, frame='icrs') + estr = 'must be an angle of length 2' + assert estr in str(excinfo.value) + + # Test valid boundary values: + for lon_range in [[-30, 30] * u.deg, [20, 560] * u.deg]: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(longitude_range=lon_range, frame='icrs') + estr = 'Longitude values must be within [0, 360] degrees or equivalent!' + assert estr in str(excinfo.value) + + def test_invalid_lat_range(self): + # Test input types: + for lat_range in [u.Quantity([0, 10], u.m / u.s), [0 * u.m / u.s, 30 * u.m / u.s]]: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(latitude_range=lat_range, frame='icrs') + estr = 'must have angular units' + assert estr in str(excinfo.value) + + for lat_range in [u.Quantity([0], u.deg), Angle([2 * u.radian])]: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(latitude_range=lat_range, frame='icrs') + estr = 'must be length 2' + assert estr in str(excinfo.value) + + for lat_range in [[0 * u.deg], [0, 10], 'string']: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(latitude_range=lat_range, frame='icrs') + estr = 'must be an angle of length 2' + assert estr in str(excinfo.value) + + # Test valid boundary values: + for lat_range in [[-120, -30] * u.deg, [20, 100] * u.deg]: + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(latitude_range=lat_range, frame='icrs') + estr = 'Latitude values must be within [-90, 90] degrees or equivalent!' + assert estr in str(excinfo.value) + + def test_invalid_lon_bounds_input(self): + invalid_bounds = CircleSphericalSkyRegion(SkyCoord(5 * u.deg, 0 * u.deg), + 0.2 * u.deg), + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(_longitude_bounds=invalid_bounds, + frame='icrs') + estr = 'Invalid direct longitude bounds input!' + assert estr in str(excinfo.value) + + def test_invalid_lat_bounds_input(self): + invalid_bounds = LuneSphericalSkyRegion(SkyCoord(3 * u.deg, 0 * u.deg), + SkyCoord(178 * u.deg, 0 * u.deg)) + with pytest.raises(ValueError) as excinfo: + _ = RangeSphericalSkyRegion(_latitude_bounds=invalid_bounds, + frame='icrs') + estr = 'Invalid direct latitude bounds input!' + assert estr in str(excinfo.value) + + def test_lon_only(self): + reg = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + frame='icrs') + assert reg._bound_nverts == 2 + assert reg.latitude_bounds is None + + coos = self.inside + self.outside + lon, lat = zip(*coos) + skycoord = SkyCoord(list(lon), list(lat)) + actual = reg.contains(skycoord) + assert_equal(actual[:len(self.inside)], True) + assert_equal(actual[len(self.inside):], False) + + assert reg.contains(skycoord[0]) + + reg_transf = reg.transform_to('galactic') + skycoord_transf = skycoord.transform_to('galactic') + actual = reg_transf.contains(skycoord_transf) + assert_equal(actual[:len(self.inside)], True) + assert_equal(actual[len(self.inside):], False) + + # Manual "no latitude constraints" + reg2 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[-90, 90] * u.deg, + frame='icrs') + assert reg2._bound_nverts == 2 + actual = reg2.contains(skycoord) + assert_equal(actual[:len(self.inside)], True) + assert_equal(actual[len(self.inside):], False) + + assert reg2.contains(skycoord[0]) + + def test_lat_only(self): + reg = RangeSphericalSkyRegion(latitude_range=[-4, 4] * u.deg, + frame='icrs') + assert reg._bound_nverts == 0 + assert reg.longitude_bounds is None + + coos = self.inside + self.outside + lon, lat = zip(*coos) + skycoord = SkyCoord(list(lon), list(lat)) + actual = reg.contains(skycoord) + assert_equal(actual[:len(self.inside)], True) + assert_equal(actual[len(self.inside):], False) + + assert reg.contains(skycoord[0]) + + reg_transf = reg.transform_to('galactic') + skycoord_transf = skycoord.transform_to('galactic') + actual = reg_transf.contains(skycoord_transf) + assert_equal(actual[:len(self.inside)], True) + assert_equal(actual[len(self.inside):], False) + + def test_lon_over_meridian(self): + reg = RangeSphericalSkyRegion(longitude_range=[358, 2] * u.deg, + latitude_range=[-4, 4] * u.deg, + frame='icrs') + assert_allclose(reg.centroid.ra.deg, 0) + assert_allclose(reg.centroid.dec.deg, 0) + + def test_lat_over_poles(self): + reg = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[80, -80] * u.deg, + frame='icrs') + assert isinstance(reg.latitude_bounds, CompoundSphericalSkyRegion) + assert reg._bound_nverts == 6 + + reg2 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[90, -80] * u.deg, + frame='icrs') + assert isinstance(reg2.latitude_bounds, CircleSphericalSkyRegion) + assert reg2.latitude_bounds.center.dec.deg == -90 + + reg3 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[80, -90] * u.deg, + frame='icrs') + assert isinstance(reg3.latitude_bounds, CircleSphericalSkyRegion) + assert reg3.latitude_bounds.center.dec.deg == 90 + + def test_centroid(self): + assert_allclose(self.reg.centroid.ra.deg, 5) + assert_allclose(self.reg.centroid.dec.deg, 0) + + # Long and narrow range: + reg = RangeSphericalSkyRegion(longitude_range=[0, 350] * u.deg, + latitude_range=[45, 50] * u.deg, + frame='icrs') + assert reg.centroid == reg.centroid_avg + + reg2 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + frame='icrs') + assert_allclose(reg2.centroid.ra.deg, 5) + assert_allclose(reg2.centroid.dec.deg, 0) + + reg3 = RangeSphericalSkyRegion(latitude_range=[-4, 4] * u.deg, + frame='icrs') + assert_allclose(reg3.centroid.ra.deg, 0) + assert_allclose(reg3.centroid.dec.deg, 90) + + def test_centroid_avg(self): + assert_allclose(self.reg.centroid_avg.ra.deg, 5) + assert_allclose(self.reg.centroid_avg.dec.deg, 0) + + # Very wide range: + reg = RangeSphericalSkyRegion(longitude_range=[200, 100] * u.deg, + latitude_range=[-50, 80] * u.deg, + frame='icrs') + assert reg.centroid_avg == SkyCoord(330 * u.deg, 15 * u.deg) + + reg2 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + frame='icrs') + assert_allclose(reg2.centroid_avg.ra.deg, 5) + assert_allclose(reg2.centroid_avg.dec.deg, 0) + + reg3 = RangeSphericalSkyRegion(latitude_range=[-4, 4] * u.deg, + frame='icrs') + assert_allclose(reg3.centroid_avg.ra.deg, 0) + assert_allclose(reg3.centroid_avg.dec.deg, 90) + + def test_vertices(self): + assert_allclose(self.reg.vertices.ra.deg, [0, 10, 10, 0]) + assert_allclose(self.reg.vertices.dec.deg, [-4, -4, 4, 4]) + + def test_copy(self): + reg = self.reg.copy() + assert_allclose(reg.centroid.ra.deg, 5) + assert_allclose(reg.centroid.dec.deg, 0) + assert reg.meta == self.meta + assert reg.visual == self.visual + + def test_transformation(self, wcs): + + # Test error for no boundary distortions: + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_pixel(wcs) + estr = 'Invalid parameter: `include_boundary_distortions=False`!' + assert estr in str(excinfo.value) + + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_sky(wcs) + estr = 'Invalid parameter: `include_boundary_distortions=False`!' + assert estr in str(excinfo.value) + + # Test boundary distortion transformations: + polypix2 = self.reg.to_pixel(wcs, + include_boundary_distortions=True, + discretize_kwargs={'n_points': 2}) + assert isinstance(polypix2, PolygonPixelRegion) + assert len(polypix2.vertices) == 8 + + assert_allclose(polypix2.vertices.x, + [-4861.625867, -4717.364168, + -4536.821565, -5079.772466, + -5685.546225, -5734.531503, + -5771.862922, -5302.114084]) + assert_allclose(polypix2.vertices.y, + [-2724.76778, -2909.927597, + -3091.161999, -3184.698562, + -3236.051353, -3037.141783, + -2837.972245, -2798.470494]) + + polysky2 = self.reg.to_sky(wcs, + include_boundary_distortions=True, + discretize_kwargs={'n_points': 2}) + assert isinstance(polysky2, PolygonSkyRegion) + assert len(polysky2.vertices) == 8 + + # WCS is Galactic: + assert_allclose(polysky2.vertices.l.deg, + [99.222517, 96.337283, + 92.726431, 103.585449, + 115.700925, 116.68063, + 117.427258, 108.032282]) + assert_allclose(polysky2.vertices.b.deg, + [-56.485356, -60.188552, + -63.81324, -65.683971, + -66.711027, -62.732836, + -58.749445, -57.95941]) + + def test_transformation_over_poles(self, wcs): + # Transformation for over-poles case: + reg = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[80, -80] * u.deg, + frame='icrs') + + with pytest.raises(NotImplementedError): + _ = reg.to_pixel(wcs, + include_boundary_distortions=True, + discretize_kwargs={'n_points': 10}) + + def test_transformation_no_wcs(self): + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_sky(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + with pytest.raises(ValueError) as excinfo: + _ = self.reg.to_pixel(include_boundary_distortions=True) + estr = "'wcs' must be set if 'include_boundary_distortions'=True" + assert estr in str(excinfo.value) + + def test_frame_transformation(self): + reg2 = self.reg.transform_to('galactic') + reg2_cent = reg2.centroid + transf_reg_cent = self.reg.centroid.transform_to('galactic') + + assert isinstance(reg2, RangeSphericalSkyRegion) + assert not reg2._is_original_frame + assert transf_reg_cent.frame.name == reg2_cent.frame.name + assert_quantity_allclose(reg2_cent.l, transf_reg_cent.l) + assert_quantity_allclose(reg2_cent.b, transf_reg_cent.b) + assert reg2.frame.name == 'galactic' + + assert reg2 != self.reg + + reg3 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[80, 90] * u.deg, + frame='icrs').transform_to('galactic') + assert reg3 != reg2 + + # Formally won't round trip, because the original lon/lat range info + # is not stored: + reg = reg2.transform_to('icrs') + assert reg != self.reg + + def test_copy_frame_transformation(self): + reg2 = self.reg.transform_to('galactic') + reg3 = reg2.copy() + + assert reg2 == reg3 + + def test_repr_str_frame_transformation(self): + expected_repr_transf = (', ' + 'center_gc2=)>,\n' + 'latitude_bounds=, ' + 'inner_radius=86.0 deg, outer_radius=94.0 deg)>\n' + ')>') + expected_str_transf = ('Region: RangeSphericalSkyRegion\n' + 'frame: galactic\n' + 'longitude_bounds: , ' + 'center_gc2=)>\n' + 'latitude_bounds: , ' + 'inner_radius=86.0 deg, outer_radius=94.0 deg)>') + + reg2 = self.reg.transform_to('galactic') + + assert repr(reg2) == expected_repr_transf + assert str(reg2) == expected_str_transf + + def test_eq(self): + reg = self.reg.copy() + assert reg == self.reg + reg.longitude_range = [0, 4] * u.deg + assert reg != self.reg + + reg2 = CircleSphericalSkyRegion(SkyCoord(5 * u.deg, 0 * u.deg), + 6.3999492595405645 * u.deg) + assert self.reg != reg2 + + def test_bounding_circle(self): + skycoord = SkyCoord(5 * u.deg, 0 * u.deg) + reg = CircleSphericalSkyRegion(skycoord, + 6.3999492595405645 * u.deg) + + bc = self.reg.bounding_circle + assert_quantity_allclose(bc.radius, reg.radius) + assert_quantity_allclose(bc.center.ra, reg.center.ra) + assert_quantity_allclose(bc.center.dec, reg.center.dec) + + reg2 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + frame='icrs') + assert reg2.bounding_circle == reg2.longitude_bounds.bounding_circle + + reg3 = RangeSphericalSkyRegion(latitude_range=[-4, 4] * u.deg, + frame='icrs') + assert reg3.bounding_circle == reg3.latitude_bounds.bounding_circle + + def test_bounding_circle_over_poles(self): + # Wrap over poles: + reg = RangeSphericalSkyRegion(longitude_range=[0, 350] * u.deg, + latitude_range=[80, -80] * u.deg, + frame='icrs') + + with pytest.raises(NotImplementedError): + _ = reg.bounding_circle + + def test_bounding_lonlat(self): + bounding_lonlat = self.reg.bounding_lonlat + + assert_quantity_allclose(bounding_lonlat[0], + Longitude([0. * u.deg, + 10. * u.deg])) + + assert_quantity_allclose(bounding_lonlat[1], + Latitude([-4 * u.deg, + 4 * u.deg])) + + # Transformed: + reg_tr = self.reg.transform_to('galactic') + bll_transf = reg_tr.bounding_lonlat + assert_quantity_allclose(bll_transf[0], + Longitude([92.7264313, 117.42725845] * u.deg)) + assert_quantity_allclose(bll_transf[1], + Latitude([-66.71102705, -56.48535559] * u.deg)) + + # Touching poles: should return longitude bounds + reg2 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[80, 90] * u.deg, + frame='icrs') + bounding_lonlat2 = reg2.bounding_lonlat + assert bounding_lonlat2[0] is not None + + reg3 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + latitude_range=[-90, -80] * u.deg, + frame='icrs') + bounding_lonlat3 = reg3.bounding_lonlat + assert bounding_lonlat3[0] is not None + + # Only lon or lat ranges: + reg4 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + frame='icrs') + bll = reg4.bounding_lonlat + bound_bll = reg4.longitude_bounds.bounding_lonlat + assert_quantity_allclose(bll[0], bound_bll[0]) + assert_quantity_allclose(bll[1], bound_bll[1]) + + reg5 = RangeSphericalSkyRegion(latitude_range=[-4, 4] * u.deg, + frame='icrs') + bll = reg5.bounding_lonlat + bound_bll = reg5.latitude_bounds.bounding_lonlat + assert bll[0] is None + assert_quantity_allclose(bll[1], bound_bll[1]) + + def test_bounding_lonlat_over_poles(self): + # Wrap over poles: + reg = RangeSphericalSkyRegion(longitude_range=[0, 350] * u.deg, + latitude_range=[80, -80] * u.deg, + frame='icrs') + + with pytest.raises(NotImplementedError): + _ = reg.bounding_lonlat + + def test_discretize_boundary(self): + polyrange = self.reg.discretize_boundary(n_points=100) + assert isinstance(polyrange, PolygonSphericalSkyRegion) + assert len(polyrange.vertices) == 400 + + reg2 = RangeSphericalSkyRegion(longitude_range=[0, 10] * u.deg, + frame='icrs') + polyrange2 = reg2.discretize_boundary(n_points=100) + assert polyrange2 == reg2.longitude_bounds.discretize_boundary(n_points=100) + assert len(polyrange2.vertices) == 200 + + reg3 = RangeSphericalSkyRegion(latitude_range=[-4, 4] * u.deg, + frame='icrs') + polyrange3 = reg3.discretize_boundary(n_points=100) + assert polyrange3 == reg3.latitude_bounds.discretize_boundary(n_points=100) + assert len(polyrange3.region1.vertices) == 100 + assert len(polyrange3.region2.vertices) == 100 + + def test_discretize_over_poles(self): + # Wrap over poles: + reg = RangeSphericalSkyRegion(longitude_range=[0, 350] * u.deg, + latitude_range=[80, -80] * u.deg, + frame='icrs') + + with pytest.raises(NotImplementedError): + _ = reg.discretize_boundary() diff --git a/regions/shapes/tests/test_whole_sky.py b/regions/shapes/tests/test_whole_sky.py new file mode 100644 index 00000000..9646673d --- /dev/null +++ b/regions/shapes/tests/test_whole_sky.py @@ -0,0 +1,64 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst + +import astropy.units as u +import pytest +from astropy.io import fits +from astropy.utils.data import get_pkg_data_filename +from astropy.wcs import WCS + +from regions.core import RegionMeta, RegionVisual +from regions.shapes.tests.test_common import BaseTestSphericalSkyRegion +from regions.shapes.whole_sky import WholeSphericalSkyRegion + + +@pytest.fixture(scope='session', name='wcs') +def wcs_fixture(): + filename = get_pkg_data_filename('data/example_header.fits') + header = fits.getheader(filename) + return WCS(header) + + +class TestWholeSphericalSkyRegion(BaseTestSphericalSkyRegion): + meta = RegionMeta({'text': 'test'}) + visual = RegionVisual({'color': 'blue'}) + reg = WholeSphericalSkyRegion(meta=meta, visual=visual) + inside = [(3 * u.deg, 4 * u.deg), (1 * u.deg, 3 * u.deg)] + outside = [] + expected_repr = ('') + expected_str = ('Region: WholeSphericalSkyRegion') + + def test_del(self): + # test that Region shape attributes (descriptors) cannot be + # deleted) + if len(self.reg._params) > 0: + with pytest.raises(AttributeError): + delattr(self.reg, self.reg._params[0]) + + def test_copy(self): + reg = self.reg.copy() + assert isinstance(reg, WholeSphericalSkyRegion) + assert reg.meta == self.meta + assert reg.visual == self.visual + + def test_transformation(self, wcs): + with pytest.raises(NotImplementedError): + _ = self.reg.to_sky(wcs) + + with pytest.raises(NotImplementedError): + _ = self.reg.to_pixel(wcs) + + def test_frame_transformation(self): + reg2 = self.reg.transform_to('galactic') + assert isinstance(reg2, WholeSphericalSkyRegion) + + def test_eq(self): + reg = self.reg.copy() + assert reg == self.reg + + def test_bounding_circle(self): + bounding_circle = self.reg.bounding_circle + assert bounding_circle is None + + def test_bounding_lonlat(self): + bounding_lonlat = self.reg.bounding_lonlat + assert bounding_lonlat is None diff --git a/regions/shapes/text.py b/regions/shapes/text.py index 45d503f6..05bcfddd 100644 --- a/regions/shapes/text.py +++ b/regions/shapes/text.py @@ -74,6 +74,10 @@ def to_sky(self, wcs): return TextSkyRegion(center, self.text, meta=self.meta.copy(), visual=visual.copy()) + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError + def as_artist(self, origin=(0, 0), **kwargs): """ Return a matplotlib Text object for this region @@ -143,3 +147,7 @@ def to_pixel(self, wcs): return TextPixelRegion(center, self.text, meta=self.meta.copy(), visual=visual.copy()) + + def to_spherical_sky(self, wcs=None, include_boundary_distortions=False, + discretize_kwargs=None): + raise NotImplementedError diff --git a/regions/shapes/whole_sky.py b/regions/shapes/whole_sky.py new file mode 100644 index 00000000..8e754de8 --- /dev/null +++ b/regions/shapes/whole_sky.py @@ -0,0 +1,72 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +""" +This module defines a whole-sky spherical sky region, for use with +compound region logic and over-the-pole logic. +""" + +import numpy as np + +from regions.core.core import SphericalSkyRegion +from regions.core.metadata import RegionMeta, RegionVisual + + +class WholeSphericalSkyRegion(SphericalSkyRegion): + """ + Spherical region representing the whole sky. + + Implemented to handle compound region logic, particularly with + ranges and over-the-pole pole logic. + + Parameters + ---------- + meta : `~regions.RegionMeta` or `dict`, optional + A dictionary that stores the meta attributes of the region. + visual : `~regions.RegionVisual` or `dict`, optional + A dictionary that stores the visual meta attributes of the + region. + """ + + def __init__(self, meta=None, visual=None): + self.meta = meta or RegionMeta() + self.visual = visual or RegionVisual() + + def contains(self, coord): + if coord.isscalar: + return True + else: + return np.ones(coord.shape, dtype=bool) + + @property + def bounding_circle(self): + # Not defined + return None + + @property + def bounding_lonlat(self): + # Not defined + return None + + def transform_to(self, frame): + return self.__class__() + + def discretize_boundary(self, n_points=100): + # Not defined + raise NotImplementedError('Not defined.') + + def to_sky( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None + ): + # Not defined + raise NotImplementedError('Not defined.') + + def to_pixel( + self, + wcs=None, + include_boundary_distortions=False, + discretize_kwargs=None + ): + # Not defined + raise NotImplementedError('Not defined.')