# Licensed under a 3-clause BSD style license - see LICENSE.rst
import copy
import logging
from functools import lru_cache
import numpy as np
from astropy import units as u
from astropy.coordinates import Angle, SkyCoord
from astropy.io import fits
from astropy.table import QTable, Table
from astropy.utils import lazyproperty
from astropy.visualization.wcsaxes import WCSAxes
from astropy.wcs.utils import (
proj_plane_pixel_area,
proj_plane_pixel_scales,
wcs_to_celestial_frame,
)
from regions import (
CompoundSkyRegion,
PixCoord,
PointSkyRegion,
RectanglePixelRegion,
Regions,
SkyRegion,
)
import matplotlib.pyplot as plt
from gammapy.utils.regions import (
compound_region_center,
compound_region_to_regions,
regions_to_compound_region,
)
from gammapy.visualization.utils import ARTIST_TO_LINE_PROPERTIES
from ..axes import MapAxes
from ..coord import MapCoord
from ..core import Map
from ..geom import Geom, pix_tuple_to_idx
from ..utils import _check_width
from ..wcs import WcsGeom
__all__ = ["RegionGeom"]
[docs]
class RegionGeom(Geom):
"""Map geometry representing a region on the sky.
The spatial component of the geometry is made up of a
single pixel with an arbitrary shape and size. It can
also have any number of non-spatial dimensions. This class
represents the geometry for the `RegionNDMap` maps.
Parameters
----------
region : `~regions.SkyRegion`
Region object.
axes : list of `MapAxis`
Non-spatial data axes.
wcs : `~astropy.wcs.WCS`
Optional WCS object to project the region if needed.
binsz_wcs : `~astropy.units.Quantity`
Angular bin size of the underlying `~WcsGeom` used to evaluate
quantities in the region. Default is "0.1 deg". This default
value is adequate for the majority of use cases. If a WCS object
is provided, the input of ``binsz_wcs`` is overridden.
"""
is_regular = True
is_allsky = False
is_hpx = False
is_region = True
_slice_spatial_axes = slice(0, 2)
_slice_non_spatial_axes = slice(2, None)
projection = "TAN"
def __init__(self, region, axes=None, wcs=None, binsz_wcs="0.1 deg"):
self._region = region
self._axes = MapAxes.from_default(axes, n_spatial_axes=2)
self._binsz_wcs = u.Quantity(binsz_wcs)
if wcs is None and region is not None:
if isinstance(region, CompoundSkyRegion):
self._center = compound_region_center(region)
else:
self._center = region.center
wcs = WcsGeom.create(
binsz=binsz_wcs,
skydir=self._center,
proj=self.projection,
frame=self._center.frame.name,
).wcs
self._wcs = wcs
# TODO : can we get the width before defining the wcs ?
wcs = WcsGeom.create(
binsz=binsz_wcs,
width=tuple(self.width),
skydir=self._center,
proj=self.projection,
frame=self._center.frame.name,
).wcs
self._wcs = wcs
self.ndim = len(self.data_shape)
# define cached methods
self.get_wcs_coord_and_weights = lru_cache()(self.get_wcs_coord_and_weights)
def __setstate__(self, state):
for key, value in state.items():
if key in ["get_wcs_coord_and_weights"]:
state[key] = lru_cache()(value)
self.__dict__ = state
@property
def frame(self):
"""Coordinate system, either Galactic ("galactic") or Equatorial ("icrs")."""
if self.region is None:
return "icrs"
try:
return self.region.center.frame.name
except AttributeError:
return wcs_to_celestial_frame(self.wcs).name
@property
def binsz_wcs(self):
"""Angular bin size of the underlying `~WcsGeom`.
Returns
-------
binsz_wcs: `~astropy.coordinates.Angle`
Angular bin size.
"""
return Angle(proj_plane_pixel_scales(self.wcs), unit="deg")
@lazyproperty
def _rectangle_bbox(self):
if self.region is None:
raise ValueError("Region definition required.")
regions = compound_region_to_regions(self.region)
regions_pix = [_.to_pixel(self.wcs) for _ in regions]
bbox = regions_pix[0].bounding_box
for region_pix in regions_pix[1:]:
bbox = bbox.union(region_pix.bounding_box)
try:
rectangle_pix = bbox.to_region()
except ValueError:
rectangle_pix = RectanglePixelRegion(
center=PixCoord(*bbox.center[::-1]), width=1, height=1
)
return rectangle_pix.to_sky(self.wcs)
@property
def width(self):
"""Width of bounding box of the region.
Returns
-------
width : `~astropy.units.Quantity`
Dimensions of the region in both spatial dimensions.
Units: ``deg``
"""
rectangle = self._rectangle_bbox
return u.Quantity([rectangle.width.to("deg"), rectangle.height.to("deg")])
@property
def region(self):
"""The spatial component of the region geometry as a `~regions.SkyRegion`."""
return self._region
@property
def is_all_point_sky_regions(self):
"""Whether regions are all point regions."""
regions = compound_region_to_regions(self.region)
return np.all([isinstance(_, PointSkyRegion) for _ in regions])
@property
def axes(self):
"""List of non-spatial axes."""
return self._axes
@property
def axes_names(self):
"""All axes names."""
return ["lon", "lat"] + self.axes.names
@property
def wcs(self):
"""WCS projection object."""
return self._wcs
@property
def center_coord(self):
"""Center coordinate of the region as a `astropy.coordinates.SkyCoord`."""
return self.pix_to_coord(self.center_pix)
@property
def center_pix(self):
"""Pixel values corresponding to the center of the region."""
return tuple((np.array(self.data_shape) - 1.0) / 2)[::-1]
@lazyproperty
def center_skydir(self):
"""Sky coordinate of the center of the region."""
if self.region is None:
return SkyCoord(np.nan * u.deg, np.nan * u.deg)
return self._rectangle_bbox.center
@property
def npix(self):
"""Number of spatial pixels."""
return ([1], [1])
[docs]
def contains(self, coords):
"""Check if a given map coordinate is contained in the region.
Requires the `.region` attribute to be set.
For `PointSkyRegion` the method always returns True.
Parameters
----------
coords : tuple, dict, `MapCoord` or `~astropy.coordinates.SkyCoord`
Object containing coordinate arrays we wish to check for inclusion
in the region.
Returns
-------
mask : `~numpy.ndarray`
Boolean array with the same shape as the input that indicates
which coordinates are inside the region.
"""
if self.region is None:
raise ValueError("Region definition required.")
coords = MapCoord.create(coords, frame=self.frame, axis_names=self.axes.names)
if self.is_all_point_sky_regions:
return np.ones(coords.skycoord.shape, dtype=bool)
return self.region.contains(coords.skycoord, self.wcs)
[docs]
def contains_wcs_pix(self, pix):
"""Check if a given WCS pixel coordinate is contained in the region.
For `PointSkyRegion` the method always returns True.
Parameters
----------
pix : tuple
Tuple of pixel coordinates.
Returns
-------
containment : `~numpy.ndarray`
Boolean array.
"""
if self.is_all_point_sky_regions:
return np.ones(pix[0].shape, dtype=bool)
region_pix = self.region.to_pixel(self.wcs)
return region_pix.contains(PixCoord(pix[0], pix[1]))
[docs]
def separation(self, position):
"""Angular distance between the center of the region and the given position.
Parameters
----------
position : `astropy.coordinates.SkyCoord`
Sky coordinate we want the angular distance to.
Returns
-------
sep : `~astropy.coordinates.Angle`
The on-sky separation between the given coordinate and the region center.
"""
return self.center_skydir.separation(position)
@property
def data_shape(self):
"""Shape of the `~numpy.ndarray` matching this geometry."""
return self._shape[::-1]
@property
def data_shape_axes(self):
"""Shape of data of the non-spatial axes and unit spatial axes."""
return self.axes.shape[::-1] + (1, 1)
@property
def _shape(self):
"""Number of bins in each dimension.
The spatial dimension is always (1, 1), as a
`RegionGeom` is not pixelized further.
"""
return tuple((1, 1) + self.axes.shape)
[docs]
def get_coord(self, mode="center", frame=None, sparse=False, axis_name=None):
"""Get map coordinates from the geometry.
Parameters
----------
mode : {'center', 'edges'}, optional
Get center or edge coordinates for the non-spatial axes.
Default is "center".
frame : str or `~astropy.coordinates.Frame`, optional
Coordinate frame. Default is None.
sparse : bool, optional
Compute sparse coordinates. Default is False.
axis_name : str, optional
If mode = "edges", the edges will be returned for this axis only.
Default is None.
Returns
-------
coord : `~MapCoord`
Map coordinate object.
"""
if mode == "edges" and axis_name is None:
raise ValueError("Mode 'edges' requires axis name")
coords = self.axes.get_coord(mode=mode, axis_name=axis_name)
coords["skycoord"] = self.center_skydir.reshape((1, 1))
if frame is None:
frame = self.frame
coords = MapCoord.create(coords, frame=self.frame).to_frame(frame)
if not sparse:
coords = coords.broadcasted
return coords
def _pad_spatial(self, pad_width):
raise NotImplementedError("Spatial padding of `RegionGeom` not supported")
[docs]
def crop(self):
raise NotImplementedError("Cropping of `RegionGeom` not supported")
[docs]
def solid_angle(self):
"""Get solid angle of the region.
Returns
-------
angle : `~astropy.units.Quantity`
Solid angle of the region in steradians.
Units: ``sr``
"""
if self.region is None:
raise ValueError("Region definition required.")
# compound regions do not implement area()
# so we use the mask representation and estimate the area
# from the pixels in the mask using oversampling
if isinstance(self.region, CompoundSkyRegion):
# oversample by a factor of ten
oversampling = 10.0
wcs = self.to_binsz_wcs(self.binsz_wcs / oversampling).wcs
pixel_region = self.region.to_pixel(wcs)
mask = pixel_region.to_mask()
area = np.count_nonzero(mask) / oversampling**2
else:
# all other types of regions should implement area
area = self.region.to_pixel(self.wcs).area
solid_angle = area * proj_plane_pixel_area(self.wcs) * u.deg**2
return solid_angle.to("sr")
[docs]
def bin_volume(self):
"""If the `RegionGeom` has a non-spatial axis, it returns the volume of the region.
If not, it returns the solid angle size.
Returns
-------
volume : `~astropy.units.Quantity`
Volume of the region.
"""
bin_volume = self.solid_angle() * np.ones(self.data_shape)
for idx, ax in enumerate(self.axes):
shape = self.ndim * [1]
shape[-(idx + 3)] = -1
bin_volume = bin_volume * ax.bin_width.reshape(tuple(shape))
return bin_volume
[docs]
def to_wcs_geom(self, width_min=None):
"""Get the minimal equivalent geometry which contains the region.
Parameters
----------
width_min : `~astropy.quantity.Quantity`, optional
Minimum width for the resulting geometry. Can be a single number or two,
for different minimum widths in each spatial dimension.
Default is None.
Returns
-------
wcs_geom : `~WcsGeom`
A WCS geometry object.
"""
if width_min is not None:
width = np.max(
[self.width.to_value("deg"), _check_width(width_min)], axis=0
)
else:
width = self.width
wcs_geom_region = WcsGeom(wcs=self.wcs)
wcs_geom = wcs_geom_region.cutout(position=self.center_skydir, width=width)
wcs_geom = wcs_geom.to_cube(self.axes)
return wcs_geom
[docs]
def to_binsz_wcs(self, binsz):
"""Change the bin size of the underlying WCS geometry.
Parameters
----------
binsz : float, str or `~astropy.quantity.Quantity`
Bin size.
Returns
-------
region : `~RegionGeom`
Region geometry with the same axes and region as the input,
but different WCS pixelization.
"""
new_geom = RegionGeom(self.region, axes=self.axes, binsz_wcs=binsz)
return new_geom
[docs]
def get_wcs_coord_and_weights(self, factor=10):
"""Get the array of spatial coordinates and corresponding weights.
The coordinates are the center of a pixel that intersects the region and
the weights that represent which fraction of the pixel is contained
in the region.
Parameters
----------
factor : int, optional
Oversampling factor to compute the weights.
Default is 10.
Returns
-------
region_coord : `~MapCoord`
MapCoord object with the coordinates inside
the region.
weights : `~numpy.ndarray`
Weights representing the fraction of each pixel
contained in the region.
"""
wcs_geom = self.to_wcs_geom()
weights = wcs_geom.to_image().region_weights(
regions=[self.region], oversampling_factor=factor
)
mask = weights.data > 0
weights = weights.data[mask]
# Get coordinates
coords = wcs_geom.get_coord(sparse=True).apply_mask(mask)
return coords, weights
[docs]
def to_binsz(self, binsz):
"""Return self."""
return self
[docs]
def to_cube(self, axes):
"""Append non-spatial axes to create a higher-dimensional geometry.
Returns
-------
region : `~RegionGeom`
Region geometry with the added axes.
"""
axes = copy.deepcopy(self.axes) + axes
return self._init_copy(region=self.region, wcs=self.wcs, axes=axes)
[docs]
def to_image(self):
"""Remove non-spatial axes to create a 2D region.
Returns
-------
region : `~RegionGeom`
Region geometry without any non-spatial axes.
"""
return self._init_copy(region=self.region, wcs=self.wcs, axes=None)
[docs]
def upsample(self, factor, axis_name=None):
"""Upsample a non-spatial dimension of the region by a given factor.
Returns
-------
region : `~RegionGeom`
Region geometry with the upsampled axis.
"""
axes = self.axes.upsample(factor=factor, axis_name=axis_name)
return self._init_copy(region=self.region, wcs=self.wcs, axes=axes)
[docs]
def downsample(self, factor, axis_name):
"""Downsample a non-spatial dimension of the region by a given factor.
Returns
-------
region : `~RegionGeom`
Region geometry with the downsampled axis.
"""
axes = self.axes.downsample(factor=factor, axis_name=axis_name)
return self._init_copy(region=self.region, wcs=self.wcs, axes=axes)
[docs]
def pix_to_coord(self, pix):
lon = np.where(
(-0.5 < pix[0]) & (pix[0] < 0.5),
self.center_skydir.data.lon,
np.nan * u.deg,
)
lat = np.where(
(-0.5 < pix[1]) & (pix[1] < 0.5),
self.center_skydir.data.lat,
np.nan * u.deg,
)
coords = (lon, lat)
for p, ax in zip(pix[self._slice_non_spatial_axes], self.axes):
coords += (ax.pix_to_coord(p),)
return coords
[docs]
def pix_to_idx(self, pix, clip=False):
idxs = list(pix_tuple_to_idx(pix))
for i, idx in enumerate(idxs[self._slice_non_spatial_axes]):
if clip:
np.clip(idx, 0, self.axes[i].nbin - 1, out=idx)
else:
np.putmask(idx, (idx < 0) | (idx >= self.axes[i].nbin), -1)
return tuple(idxs)
[docs]
def coord_to_pix(self, coords):
# inherited docstring
if isinstance(coords, tuple) and len(coords) == len(self.axes):
skydir = self.center_skydir.transform_to(self.frame)
coords = (skydir.data.lon, skydir.data.lat) + coords
elif isinstance(coords, dict):
valid_keys = ["lon", "lat", "skycoord"]
if not any([_ in coords for _ in valid_keys]):
coords.setdefault("skycoord", self.center_skydir)
coords = MapCoord.create(coords, frame=self.frame, axis_names=self.axes.names)
if self.region is None:
pix = (0, 0)
else:
in_region = self.contains(coords.skycoord)
x = np.zeros(coords.skycoord.shape)
x[~in_region] = np.nan
y = np.zeros(coords.skycoord.shape)
y[~in_region] = np.nan
pix = (x, y)
pix += self.axes.coord_to_pix(coords)
return pix
[docs]
def get_idx(self):
idxs = [np.arange(n, dtype=float) for n in self.data_shape[::-1]]
return np.meshgrid(*idxs[::-1], indexing="ij")[::-1]
def _make_bands_cols(self):
return []
[docs]
@classmethod
def create(cls, region, **kwargs):
"""Create region geometry.
The input region can be passed in the form of a ds9 string and will be parsed
internally by `~regions.Regions.parse`. See:
* https://astropy-regions.readthedocs.io/en/stable/region_io.html
* http://ds9.si.edu/doc/ref/region.html
Parameters
----------
region : str or `~regions.SkyRegion`
Region definition.
**kwargs : dict
Keyword arguments passed to `RegionGeom.__init__`.
Returns
-------
geom : `RegionGeom`
Region geometry.
"""
return cls.from_regions(regions=region, **kwargs)
def __str__(self):
axes = ["lon", "lat"] + [_.name for _ in self.axes]
try:
frame = self.center_skydir.frame.name
lon = self.center_skydir.data.lon.deg
lat = self.center_skydir.data.lat.deg
except AttributeError:
frame, lon, lat = "", np.nan, np.nan
return (
f"{self.__class__.__name__}\n\n"
f"\tregion : {self.region.__class__.__name__}\n"
f"\taxes : {axes}\n"
f"\tshape : {self.data_shape[::-1]}\n"
f"\tndim : {self.ndim}\n"
f"\tframe : {frame}\n"
f"\tcenter : {lon:.1f} deg, {lat:.1f} deg\n"
)
[docs]
def is_allclose(self, other, rtol_axes=1e-6, atol_axes=1e-6):
"""Compare two data IRFs for equivalency.
Parameters
----------
other : `RegionGeom`
Region geometry to compare against.
rtol_axes : float, optional
Relative tolerance for the axes comparison.
Default is 1e-6.
atol_axes : float, optional
Relative tolerance for the axes comparison.
Default is 1e-6.
Returns
-------
is_allclose : bool
Whether the geometry is all close.
"""
if not isinstance(other, self.__class__):
return TypeError(f"Cannot compare {type(self)} and {type(other)}")
if self.data_shape != other.data_shape:
return False
axes_eq = self.axes.is_allclose(other.axes, rtol=rtol_axes, atol=atol_axes)
# TODO: compare regions based on masks...
regions_eq = True
return axes_eq and regions_eq
def __eq__(self, other):
if not isinstance(other, self.__class__):
return False
return self.is_allclose(other=other)
def _to_region_table(self):
"""Export region to a FITS region table."""
if self.region is None:
raise ValueError("Region definition required.")
region_list = compound_region_to_regions(self.region)
pixel_region_list = []
for reg in region_list:
pixel_region_list.append(reg.to_pixel(self.wcs))
table = Regions(pixel_region_list).serialize(format="fits")
header = WcsGeom(wcs=self.wcs).to_header()
table.meta.update(header)
return table
[docs]
def to_hdulist(self, format="ogip", hdu_bands=None, hdu_region=None):
"""Convert geometry to HDU list.
Parameters
----------
format : {"ogip", "gadf", "ogip-sherpa"}
HDU format. Default is "ogip".
hdu_bands : str, optional
Name or index of the HDU with the BANDS table.
Default is None.
hdu_region : str, optional
Name or index of the HDU with the region table.
Not used for the "gadf" format.
Default is None.
Returns
-------
hdulist : `~astropy.io.fits.HDUList`
HDU list.
"""
if hdu_bands is None:
hdu_bands = "HDU_BANDS"
if hdu_region is None:
hdu_region = "HDU_REGION"
if format != "gadf":
hdu_region = "REGION"
hdulist = fits.HDUList()
hdulist.append(self.axes.to_table_hdu(hdu_bands=hdu_bands, format=format))
# region HDU
if self.region:
region_table = self._to_region_table()
region_hdu = fits.BinTableHDU(region_table, name=hdu_region)
hdulist.append(region_hdu)
return hdulist
[docs]
@classmethod
def from_regions(cls, regions, **kwargs):
"""Create region geometry from list of regions.
The regions are combined with union to a compound region.
Parameters
----------
regions : list of `~regions.SkyRegion` or str
Regions.
**kwargs: dict
Keyword arguments forwarded to `RegionGeom`.
Returns
-------
geom : `RegionGeom`
Region map geometry.
"""
if isinstance(regions, str):
regions = Regions.parse(data=regions, format="ds9")
elif isinstance(regions, SkyRegion):
regions = [regions]
elif isinstance(regions, SkyCoord):
regions = [PointSkyRegion(center=regions)]
elif isinstance(regions, list) and len(regions) == 0:
regions = None
if regions:
regions = regions_to_compound_region(regions)
return cls(region=regions, **kwargs)
[docs]
@classmethod
def from_hdulist(cls, hdulist, format="ogip", hdu=None):
"""Read region table and convert it to region list.
Parameters
----------
hdulist : `~astropy.io.fits.HDUList`
HDU list.
format : {"ogip", "ogip-arf", "gadf"}
HDU format. Default is "ogip".
hdu : str, optional
Name of the HDU. Default is None.
Returns
-------
geom : `RegionGeom`
Region map geometry.
"""
region_hdu = "REGION"
if format == "gadf" and hdu:
region_hdu = f"{hdu}_{region_hdu}"
if region_hdu in hdulist:
try:
region_table = QTable.read(hdulist[region_hdu])
regions_pix = Regions.parse(data=region_table, format="fits")
except TypeError:
region_table = Table.read(hdulist[region_hdu])
regions_pix = Regions.parse(data=region_table, format="fits")
wcs = WcsGeom.from_header(region_table.meta).wcs
regions = []
for region_pix in regions_pix:
# TODO: remove workaround once regions issue with fits serialization is sorted out
# see https://github.com/astropy/regions/issues/400
# requires update regions to >=v0.6
region_pix.meta["include"] = True
regions.append(region_pix.to_sky(wcs))
region = regions_to_compound_region(regions)
else:
region, wcs = None, None
if format == "ogip":
hdu_bands = "EBOUNDS"
elif format == "ogip-arf":
hdu_bands = "SPECRESP"
elif format == "gadf":
hdu_bands = hdu + "_BANDS"
else:
raise ValueError(f"Unknown format {format}")
axes = MapAxes.from_table_hdu(hdulist[hdu_bands], format=format)
return cls(region=region, wcs=wcs, axes=axes)
[docs]
def union(self, other):
"""Stack a `RegionGeom` by making the union."""
if not self == other:
raise ValueError("Can only make union if extra axes are equivalent.")
if other.region:
if self.region:
self._region = self.region.union(other.region)
else:
self._region = other.region
[docs]
def plot_region(self, ax=None, kwargs_point=None, path_effect=None, **kwargs):
"""Plot region in the sky.
Parameters
----------
ax : `~astropy.visualization.WCSAxes`, optional
Axes to plot on. If no axes are given,
the region is shown using the minimal
equivalent WCS geometry.
Default is None.
kwargs_point : dict, optional
Keyword arguments passed to `~matplotlib.lines.Line2D` for plotting
of point sources. Default is None.
path_effect : `~matplotlib.patheffects.PathEffect`, optional
Path effect applied to artists and lines.
Default is None.
**kwargs : dict
Keyword arguments forwarded to `~regions.PixelRegion.as_artist`.
Returns
-------
ax : `~astropy.visualization.WCSAxes`
Axes to plot on.
"""
if self.region:
kwargs_point = kwargs_point or {}
if ax is None:
ax = plt.gca()
if not isinstance(ax, WCSAxes):
ax.remove()
wcs_geom = self.to_wcs_geom()
m = Map.from_geom(geom=wcs_geom.to_image())
ax = m.plot(add_cbar=False, vmin=-1, vmax=0)
kwargs.setdefault("facecolor", "None")
kwargs.setdefault("edgecolor", "tab:blue")
kwargs_point.setdefault("marker", "*")
for key, value in kwargs.items():
key_point = ARTIST_TO_LINE_PROPERTIES.get(key, None)
if key_point:
kwargs_point[key_point] = value
for region in compound_region_to_regions(self.region):
region_pix = region.to_pixel(wcs=ax.wcs)
if isinstance(region, PointSkyRegion):
artist = region_pix.as_artist(**kwargs_point)
else:
artist = region_pix.as_artist(**kwargs)
if path_effect:
artist.add_path_effect(path_effect)
ax.add_artist(artist)
return ax
else:
logging.info("Region definition required.")
