# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""
Tools for converting between `regions.Region` and Aperture objects and
between `shapely.Polygon` and `regions.PolygonRegion` objects.
"""
import astropy.units as u
import numpy as np
from photutils.aperture.circle import (CircularAnnulus, CircularAperture,
SkyCircularAnnulus, SkyCircularAperture)
from photutils.aperture.core import Aperture
from photutils.aperture.ellipse import (EllipticalAnnulus, EllipticalAperture,
SkyEllipticalAnnulus,
SkyEllipticalAperture)
from photutils.aperture.polygon import PolygonAperture, SkyPolygonAperture
from photutils.aperture.rectangle import (RectangularAnnulus,
RectangularAperture,
SkyRectangularAnnulus,
SkyRectangularAperture)
__all__ = ['aperture_to_region', 'region_to_aperture']
__doctest_requires__ = {'region_to_aperture': ['regions >= 0.12.dev'],
'aperture_to_region': ['regions >= 0.12.dev'],
'_scalar_aperture_to_region': ['regions >= 0.12.dev'],
'_shapely_polygon_to_region': ['regions >= 0.12.dev',
'shapely']}
[docs]
def region_to_aperture(region):
"""
Convert a given `regions.Region` object to an
`~photutils.aperture.Aperture` object.
Parameters
----------
region : `regions.Region`
A supported `regions.Region` object.
Returns
-------
aperture : `~photutils.aperture.Aperture`
An equivalent ``photutils`` aperture.
Raises
------
`TypeError`
The given `regions.Region` object is not supported.
Notes
-----
The ellipse ``width`` and ``height`` region parameters represent
the full extent of the shapes and thus are divided by 2 when
converting to elliptical aperture objects, which are defined using
the semi-major (``a``) and semi-minor (``b``) axes. The ``width``
and ``height`` parameters are mapped to the semi-major (``a``) and
semi-minor (``b``) axes parameters, respectively, of the elliptical
apertures.
The region ``angle`` for sky-based regions is defined as the angle
of the ``width`` axis relative to WCS longitude axis (PA=90).
However, the sky-based apertures define the ``theta`` as the
position angle of the semimajor axis relative to the North celestial
pole (PA=0). Therefore, for sky-based regions the region ``angle``
is converted to the aperture ``theta`` parameter by subtracting 90
degrees.
.. |rarr| unicode:: U+0279E .. RIGHTWARDS ARROW
The following `regions.Region` objects are supported, shown with
their equivalent `~photutils.aperture.Aperture` object:
* `~regions.CirclePixelRegion` |rarr|
`~photutils.aperture.CircularAperture`
* `~regions.CircleSkyRegion` |rarr|
`~photutils.aperture.SkyCircularAperture`
* `~regions.EllipsePixelRegion` |rarr|
`~photutils.aperture.EllipticalAperture`
* `~regions.EllipseSkyRegion` |rarr|
`~photutils.aperture.SkyEllipticalAperture`
* `~regions.RectanglePixelRegion` |rarr|
`~photutils.aperture.RectangularAperture`
* `~regions.RectangleSkyRegion` |rarr|
`~photutils.aperture.SkyRectangularAperture`
* `~regions.CircleAnnulusPixelRegion` |rarr|
`~photutils.aperture.CircularAnnulus`
* `~regions.CircleAnnulusSkyRegion` |rarr|
`~photutils.aperture.SkyCircularAnnulus`
* `~regions.EllipseAnnulusPixelRegion` |rarr|
`~photutils.aperture.EllipticalAnnulus`
* `~regions.EllipseAnnulusSkyRegion` |rarr|
`~photutils.aperture.SkyEllipticalAnnulus`
* `~regions.RectangleAnnulusPixelRegion` |rarr|
`~photutils.aperture.RectangularAnnulus`
* `~regions.RectangleAnnulusSkyRegion` |rarr|
`~photutils.aperture.SkyRectangularAnnulus`
* `~regions.PolygonPixelRegion` |rarr|
`~photutils.aperture.PolygonAperture`
* `~regions.PolygonSkyRegion` |rarr|
`~photutils.aperture.SkyPolygonAperture`
Examples
--------
>>> from regions import CirclePixelRegion, PixCoord
>>> from photutils.aperture import region_to_aperture
>>> region = CirclePixelRegion(center=PixCoord(x=10, y=20), radius=5)
>>> aperture = region_to_aperture(region)
>>> aperture
<CircularAperture([10., 20.], r=5.0)>
"""
from regions import (CircleAnnulusPixelRegion, CircleAnnulusSkyRegion,
CirclePixelRegion, CircleSkyRegion,
EllipseAnnulusPixelRegion, EllipseAnnulusSkyRegion,
EllipsePixelRegion, EllipseSkyRegion,
PolygonPixelRegion, PolygonSkyRegion,
RectangleAnnulusPixelRegion,
RectangleAnnulusSkyRegion, RectanglePixelRegion,
RectangleSkyRegion, Region)
if not isinstance(region, Region):
msg = 'Input region must be a Region object'
raise TypeError(msg)
if isinstance(region, CirclePixelRegion):
aperture = CircularAperture(region.center.xy, region.radius)
elif isinstance(region, CircleSkyRegion):
aperture = SkyCircularAperture(region.center, region.radius)
elif isinstance(region, EllipsePixelRegion):
aperture = EllipticalAperture(
region.center.xy, region.width * 0.5, region.height * 0.5,
theta=region.angle)
elif isinstance(region, EllipseSkyRegion):
aperture = SkyEllipticalAperture(
region.center, region.width * 0.5, region.height * 0.5,
theta=(region.angle - (90 * u.deg)))
elif isinstance(region, RectanglePixelRegion):
aperture = RectangularAperture(
region.center.xy, region.width, region.height,
theta=region.angle)
elif isinstance(region, RectangleSkyRegion):
aperture = SkyRectangularAperture(
region.center, region.width, region.height,
theta=(region.angle - (90 * u.deg)))
elif isinstance(region, CircleAnnulusPixelRegion):
aperture = CircularAnnulus(
region.center.xy, region.inner_radius, region.outer_radius)
elif isinstance(region, CircleAnnulusSkyRegion):
aperture = SkyCircularAnnulus(
region.center, region.inner_radius, region.outer_radius)
elif isinstance(region, EllipseAnnulusPixelRegion):
aperture = EllipticalAnnulus(
region.center.xy,
region.inner_width * 0.5, region.outer_width * 0.5,
region.outer_height * 0.5, b_in=region.inner_height * 0.5,
theta=region.angle)
elif isinstance(region, EllipseAnnulusSkyRegion):
aperture = SkyEllipticalAnnulus(
region.center,
region.inner_width * 0.5, region.outer_width * 0.5,
region.outer_height * 0.5, b_in=region.inner_height * 0.5,
theta=(region.angle - (90 * u.deg)))
elif isinstance(region, RectangleAnnulusPixelRegion):
aperture = RectangularAnnulus(
region.center.xy,
region.inner_width, region.outer_width,
region.outer_height, h_in=region.inner_height,
theta=region.angle)
elif isinstance(region, RectangleAnnulusSkyRegion):
aperture = SkyRectangularAnnulus(
region.center,
region.inner_width, region.outer_width,
region.outer_height, h_in=region.inner_height,
theta=(region.angle - (90 * u.deg)))
elif isinstance(region, PolygonPixelRegion):
verts = np.stack([region.vertices.x, region.vertices.y], axis=1)
aperture = PolygonAperture.from_vertices(verts)
elif isinstance(region, PolygonSkyRegion):
aperture = SkyPolygonAperture.from_vertices(region.vertices)
else:
msg = (f'Cannot convert {region.__class__.__name__!r} to an '
'Aperture object')
raise TypeError(msg)
return aperture
[docs]
def aperture_to_region(aperture):
"""
Convert a given `~photutils.aperture.Aperture` object to a
`regions.Region` or `regions.Regions` object.
Because a `regions.Region` object can only have one position, a
`regions.Regions` object will be returned if the input ``aperture``
has more than one position. Otherwise, a `regions.Region` object
will be returned.
Parameters
----------
aperture : `~photutils.aperture.Aperture`
An `~photutils.aperture.Aperture` object to convert.
Returns
-------
region : `regions.Region` or `regions.Regions`
An equivalent `regions.Region` object. If the input ``aperture``
has more than one position then a `regions.Regions` will be
returned.
Notes
-----
The elliptical aperture ``a`` and ``b`` parameters represent the
semi-major and semi-minor axes, respectively. The ``a`` and ``b``
parameters are mapped to the ellipse ``width`` and ``height`` region
parameters, respectively, by multiplying by 2 because they represent
the full extent of the ellipse.
The region ``angle`` for sky-based regions is defined as the angle
of the ``width`` axis relative to WCS longitude axis (PA=90).
However, the sky-based apertures define the ``theta`` as the
position angle of the semimajor axis relative to the North celestial
pole (PA=0). Therefore, for sky-based apertures the ``theta``
parameter is converted to the region ``angle`` by adding 90 degrees.
.. |rarr| unicode:: U+0279E .. RIGHTWARDS ARROW
The following `~photutils.aperture.Aperture` objects are supported,
shown with their equivalent `regions.Region` object:
* `~photutils.aperture.CircularAperture` |rarr|
`~regions.CirclePixelRegion`
* `~photutils.aperture.SkyCircularAperture` |rarr|
`~regions.CircleSkyRegion`
* `~photutils.aperture.EllipticalAperture` |rarr|
`~regions.EllipsePixelRegion`
* `~photutils.aperture.SkyEllipticalAperture` |rarr|
`~regions.EllipseSkyRegion`
* `~photutils.aperture.RectangularAperture` |rarr|
`~regions.RectanglePixelRegion`
* `~photutils.aperture.SkyRectangularAperture` |rarr|
`~regions.RectangleSkyRegion`
* `~photutils.aperture.CircularAnnulus` |rarr|
`~regions.CircleAnnulusPixelRegion`
* `~photutils.aperture.SkyCircularAnnulus` |rarr|
`~regions.CircleAnnulusSkyRegion`
* `~photutils.aperture.EllipticalAnnulus` |rarr|
`~regions.EllipseAnnulusPixelRegion`
* `~photutils.aperture.SkyEllipticalAnnulus` |rarr|
`~regions.EllipseAnnulusSkyRegion`
* `~photutils.aperture.RectangularAnnulus` |rarr|
`~regions.RectangleAnnulusPixelRegion`
* `~photutils.aperture.SkyRectangularAnnulus` |rarr|
`~regions.RectangleAnnulusSkyRegion`
* `~photutils.aperture.PolygonAperture` |rarr|
`~regions.PolygonPixelRegion`
* `~photutils.aperture.SkyPolygonAperture` |rarr|
`~regions.PolygonSkyRegion`
Examples
--------
>>> from photutils.aperture import CircularAperture, aperture_to_region
>>> aperture = CircularAperture((10, 20), r=5)
>>> region = aperture_to_region(aperture)
>>> region
<CirclePixelRegion(center=PixCoord(x=10.0, y=20.0), radius=5.0)>
>>> aperture = CircularAperture(((10, 20), (30, 40)), r=5)
>>> region = aperture_to_region(aperture)
>>> region
<Regions([
<CirclePixelRegion(center=PixCoord(x=10.0, y=20.0), radius=5.0)>,
<CirclePixelRegion(center=PixCoord(x=30.0, y=40.0), radius=5.0)>
])>
"""
from regions import Regions
if not isinstance(aperture, Aperture):
msg = 'Input aperture must be an Aperture object'
raise TypeError(msg)
if aperture.shape == ():
return _scalar_aperture_to_region(aperture)
# Multiple aperture positions return a Regions object
regs = [_scalar_aperture_to_region(aper) for aper in aperture]
return Regions(regs)
def _scalar_aperture_to_region(aperture):
"""
Convert a given scalar `~photutils.aperture.Aperture` object to a
`regions.Region` object.
Parameters
----------
aperture : `~photutils.aperture.Aperture`
An `~photutils.aperture.Aperture` object to convert. The
``aperture`` must have a single position (scalar).
Returns
-------
region : `regions.Region` or `regions.Regions`
An equivalent `regions.Region` object.
"""
from regions import (CircleAnnulusPixelRegion, CircleAnnulusSkyRegion,
CirclePixelRegion, CircleSkyRegion,
EllipseAnnulusPixelRegion, EllipseAnnulusSkyRegion,
EllipsePixelRegion, EllipseSkyRegion, PixCoord,
PolygonPixelRegion, PolygonSkyRegion,
RectangleAnnulusPixelRegion,
RectangleAnnulusSkyRegion, RectanglePixelRegion,
RectangleSkyRegion)
if aperture.shape != ():
msg = 'Only scalar (single-position) apertures are supported'
raise ValueError(msg)
if isinstance(aperture, CircularAperture):
region = CirclePixelRegion(PixCoord(*aperture.positions), aperture.r)
elif isinstance(aperture, SkyCircularAperture):
region = CircleSkyRegion(aperture.positions, aperture.r)
elif isinstance(aperture, EllipticalAperture):
region = EllipsePixelRegion(
PixCoord(*aperture.positions), aperture.a * 2, aperture.b * 2,
angle=aperture.theta)
elif isinstance(aperture, SkyEllipticalAperture):
region = EllipseSkyRegion(
aperture.positions, aperture.a * 2, aperture.b * 2,
angle=(aperture.theta + (90 * u.deg)))
elif isinstance(aperture, RectangularAperture):
region = RectanglePixelRegion(
PixCoord(*aperture.positions), aperture.w, aperture.h,
angle=aperture.theta)
elif isinstance(aperture, SkyRectangularAperture):
region = RectangleSkyRegion(
aperture.positions, aperture.w, aperture.h,
angle=(aperture.theta + (90 * u.deg)))
elif isinstance(aperture, CircularAnnulus):
region = CircleAnnulusPixelRegion(
PixCoord(*aperture.positions), aperture.r_in, aperture.r_out)
elif isinstance(aperture, SkyCircularAnnulus):
region = CircleAnnulusSkyRegion(
aperture.positions, aperture.r_in, aperture.r_out)
elif isinstance(aperture, EllipticalAnnulus):
region = EllipseAnnulusPixelRegion(
PixCoord(*aperture.positions),
aperture.a_in * 2, aperture.a_out * 2,
aperture.b_in * 2, aperture.b_out * 2,
angle=aperture.theta)
elif isinstance(aperture, SkyEllipticalAnnulus):
region = EllipseAnnulusSkyRegion(
aperture.positions,
aperture.a_in * 2, aperture.a_out * 2,
aperture.b_in * 2, aperture.b_out * 2,
angle=(aperture.theta + (90 * u.deg)))
elif isinstance(aperture, RectangularAnnulus):
region = RectangleAnnulusPixelRegion(
PixCoord(*aperture.positions),
aperture.w_in, aperture.w_out,
aperture.h_in, aperture.h_out,
angle=aperture.theta)
elif isinstance(aperture, SkyRectangularAnnulus):
region = RectangleAnnulusSkyRegion(
aperture.positions,
aperture.w_in, aperture.w_out,
aperture.h_in, aperture.h_out,
angle=(aperture.theta + (90 * u.deg)))
elif isinstance(aperture, PolygonAperture):
verts = aperture.vertices
region = PolygonPixelRegion(
PixCoord(x=verts[:, 0], y=verts[:, 1]))
elif isinstance(aperture, SkyPolygonAperture):
region = PolygonSkyRegion(aperture.vertices)
else:
msg = 'Cannot convert input aperture to a Region object'
raise TypeError(msg)
return region
def _shapely_polygon_to_region(polygon, *, label=None, visual_kwargs=None):
"""
Convert a `shapely.Polygon` object to a `regions.PolygonPixelRegion`
object.
Parameters
----------
polygon : `shapely.Polygon` or `shapely.MultiPolygon`
A Shapely Polygon or MultiPolygon object.
label : str or `None`, optional
A label for the region. If provided, it will be stored in the
meta attribute of the returned `regions.PolygonPixelRegion`
objects.
visual_kwargs : dict or `None`, optional
A dictionary of visual keyword arguments to pass to
`regions.RegionVisual`. If provided, the visual attributes will
be set on the returned region(s).
Returns
-------
result : list of `regions.PolygonPixelRegion` or `regions.Regions`
If the polygon is a `shapely.Polygon`, then a
`regions.PolygonPixelRegion` object is returned. If the polygon
is a `shapely.MultiPolygon`, then a `regions.Regions` object is
returned containing one or more `regions.PolygonPixelRegion`
objects.
Notes
-----
The `regions.PolygonPixelRegion` object does not include the
last Shapely vertex, which is the same as the first vertex. The
`regions.PolygonPixelRegion` does not need to include the last
vertex to close the polygon.
Examples
--------
>>> from shapely import Polygon
>>> from photutils.aperture.converters import _shapely_polygon_to_region
>>> polygon = Polygon([(1, 1), (3, 1), (2, 4), (1, 2)])
>>> region = _shapely_polygon_to_region(polygon)
>>> region
<PolygonPixelRegion(vertices=PixCoord(x=[1.0 3.0 2.0 1.0],
y=[1.0 1.0 4.0 2.0]))>
"""
from regions import PixCoord, PolygonPixelRegion, Regions, RegionVisual
from shapely import MultiPolygon, Polygon
meta = {'label': label} if label is not None else None
visual = RegionVisual(visual_kwargs) if visual_kwargs else None
if isinstance(polygon, Polygon):
x, y = np.transpose(polygon.exterior.coords[:-1])
return PolygonPixelRegion(vertices=PixCoord(x=x, y=y), meta=meta,
visual=visual)
if isinstance(polygon, MultiPolygon):
geoms = []
for poly in polygon.geoms:
x, y = np.transpose(poly.exterior.coords[:-1])
geoms.append(PolygonPixelRegion(vertices=PixCoord(x=x, y=y),
meta=meta, visual=visual))
return Regions(geoms)
msg = 'Input must be a Polygon or MultiPolygon object'
raise TypeError(msg)
