Add spherical sky region classes

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
 ---------
 

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")

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: <SkyCoord (ICRS): (ra, dec) in deg
+        (42., 43.)>
+    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])

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
 <https://docs.astropy.org/en/stable/wcs/wcsapi.html>`_ 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
+<https://docs.astropy.org/en/stable/wcs/wcsapi.html>`_, 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
+<https://docs.astropy.org/en/stable/wcs/wcsapi.html>`_ 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: <SkyCoord (ICRS): (ra, dec) in deg
+        (42., 43.)>
+    radius: 3.0 deg
+
+You can access its properties via attributes:
+
+.. code-block:: python
+
+   >>> region.center
+    <SkyCoord (ICRS): (ra, dec) in deg
+        (42., 43.)>
+   >>> region.radius
+   <Quantity 3. deg>
+
+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.