# 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.core.pixcoord import PixCoord
from regions.shapes.circle import CircleSphericalSkyRegion
from regions.shapes.polygon import (PolygonPixelRegion,
PolygonSphericalSkyRegion)
__all__ = ['LuneSphericalSkyRegion']
[docs]
class LuneSphericalSkyRegion(SphericalSkyRegion):
"""
Class for a spherical "lune", the intersection between two 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
[docs]
def contains(self, coord):
return self._compound_region.contains(coord)
[docs]
def discretize_boundary(self, n_points=100):
bound_verts = discretize_all_edge_boundaries(
self.vertices, self._edge_circs, n_points,
)[::-1] # inverted for lune
# TODO: properly check for CW order even for a nverts=2 spherical polygon (lune).
return PolygonSphericalSkyRegion(bound_verts, meta=self.meta.copy(),
visual=self.visual.copy())
[docs]
def to_polygon(self, n_points=100):
"""
Return a `~regions.PolygonSphericalSkyRegion` that approximates this
lune.
Parameters
----------
n_points : int, optional
The number of polygon vertices. Default is 100.
Returns
-------
polygon : `~regions.PolygonSphericalSkyRegion`
A polygon region approximating the lune.
"""
return self.discretize_boundary(n_points=n_points)
[docs]
def to_sky(
self,
wcs=None,
include_boundary_distortions=False,
n_points=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.',
)
self._validate_planar_spherical_transform(wcs, include_boundary_distortions)
# 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,
n_points=n_points,
).to_sky(wcs)
[docs]
def to_pixel(
self,
wcs=None,
include_boundary_distortions=False,
n_points=None,
):
if not include_boundary_distortions:
raise ValueError(
'Invalid parameter: `include_boundary_distortions=False`!\n'
'Transforming lune to planar pixel region is only possible by '
'including boundary distortions, as there is no analogous pixel region.',
)
self._validate_planar_spherical_transform(wcs, include_boundary_distortions)
disc_kwargs = {} if n_points is None else {'n_points': n_points}
# Requires spherical to planar projection (from WCS) and discretization
disc_bound = self.discretize_boundary(**disc_kwargs)
verts = wcs.world_to_pixel(disc_bound.vertices)
return PolygonPixelRegion(
PixCoord(*verts), meta=self.meta.copy(),
visual=self.visual.copy(),
)