Source code for gammapy.maps.region

import copy
from functools import lru_cache
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord, Angle
from astropy.io import fits
from astropy.table import Table
from astropy.wcs.utils import proj_plane_pixel_area, wcs_to_celestial_frame, proj_plane_pixel_scales
from regions import FITSRegionParser, fits_region_objects_to_table, SkyRegion, CompoundSkyRegion, PixCoord
from gammapy.utils.regions import (
    compound_region_to_list,
    list_to_compound_region,
    make_region,
    compound_region_center
)
from gammapy.maps.wcs import _check_width
from .core import MapCoord, Map
from .geom import Geom, pix_tuple_to_idx
from .axes import MapAxes, MapAxis
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 : `float` Angular bin size of the underlying `~WcsGeom` used to evaluate quantities in the region. Default size is 0.01 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 = binsz_wcs if wcs is None and region is not None: if isinstance(region, CompoundSkyRegion): center = compound_region_center(region) else: center = region.center wcs = WcsGeom.create( binsz=binsz_wcs, skydir=center, proj=self.projection, frame=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) @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` """ return Angle(proj_plane_pixel_scales(self.wcs), unit="deg") @property def _rectangle_bbox(self): if self.region is None: raise ValueError("Region definition required.") regions = compound_region_to_list(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) rectangle_pix = bbox.to_region() 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): """`~regions.SkyRegion` object that defines the spatial component of the region geometry""" return self._region @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): """(`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] @property 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. 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 Numpy 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) 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. Parameters ---------- pix : tuple Tuple of pixel coordinates. Returns ------- containment : `~numpy.ndarray` Bool array. """ 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 data array 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'} Get center or edge coordinates for the non-spatial axes. frame : str or `~astropy.coordinates.Frame` Coordinate frame sparse : bool Compute sparse coordinates axis_name : str If mode = "edges", the edges will be returned for this axis only. 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 return MapCoord.create(coords, frame=self.frame).to_frame(frame)
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 sr. Units: ``sr`` """ if self.region is None: raise ValueError("Region definition required.") 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 just retuns 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` Minimal width for the resulting geometry. Can be a single number or two, for different minimum widths in each spatial dimension. 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, npix=self.wcs.array_shape) 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 ---------- binzs : float, string or `~astropy.quantity.Quantity` Returns ------- region : `~RegionGeom` A RegionGeom 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 Oversampling factor to compute the weights Returns ------- region_coord : `~MapCoord` MapCoord object with the coordinates inside the region. weights : `~np.array` 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): """Returns self""" return self
[docs] def to_cube(self, axes): """Append non-spatial axes to create a higher-dimensional geometry. Returns ------- region : `~RegionGeom` RegionGeom with the added axes. """ axes = copy.deepcopy(self.axes) + axes return self._init_copy(axes=axes)
[docs] def to_image(self): """Remove non-spatial axes to create a 2D region. Returns ------- region : `~RegionGeom` RegionGeom without any non-spatial axes. """ return self._init_copy(axes=None)
[docs] def upsample(self, factor, axis_name): """Upsample a non-spatial dimension of the region by a given factor. Returns ------- region : `~RegionGeom` RegionGeom with the upsampled axis. """ axes = self.axes.upsample(factor=factor, axis_name=axis_name) return self._init_copy(axes=axes)
[docs] def downsample(self, factor, axis_name): """Downsample a non-spatial dimension of the region by a given factor. Returns ------- region : `~RegionGeom` RegionGeom with the downsampled axis. """ axes = self.axes.downsample(factor=factor, axis_name=axis_name) return self._init_copy(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): coords = MapCoord.create(coords, frame=self.frame, axis_names=self.axes.names) if self.region is None: pix = (0, 0) else: in_region = self.region.contains(coords.skycoord, wcs=self.wcs) 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. Parameters ---------- region : str or `~regions.SkyRegion` Region axes : list of `MapAxis` Non spatial axes. Returns ------- geom : `RegionGeom` Region geometry """ if isinstance(region, str): region = make_region(region) return cls(region, **kwargs)
def __repr__(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" ) def __eq__(self, other): # check overall shape and axes compatibility if self.data_shape != other.data_shape: return False for axis, otheraxis in zip(self.axes, other.axes): if axis != otheraxis: return False # TODO: compare regions return True def _to_region_table(self): """Export region to a FITS region table.""" if self.region is None: raise ValueError("Region definition required.") # TODO: make this a to_hdulist() method region_list = compound_region_to_list(self.region) pixel_region_list = [] for reg in region_list: pixel_region_list.append(reg.to_pixel(self.wcs)) table = fits_region_objects_to_table(pixel_region_list) header = WcsGeom(wcs=self.wcs, npix=self.wcs.array_shape).to_header() table.meta.update(header) return table
[docs] def to_hdulist(self, format="ogip", hdu_bands=None, hdu_region=None): """Convert geom to hdulist Parameters ---------- format : {"gadf", "ogip", "ogip-sherpa"} HDU format hdu : str Name of the HDU with the map data. 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 geom from list of regions The regions are combined with union to a compound region. Parameters ---------- regions : `~regions.SkyRegion` Regions **kwargs: dict Keyword arguments forwarded to `RegionGeom` Returns ------- geom : `RegionGeom` Region map geometry """ if isinstance(regions, (SkyRegion, str)): regions = [make_region(regions)] region = list_to_compound_region(regions) return cls(region, **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 Returns ------- geom : `RegionGeom` Region map geometry """ region_hdu = "REGION" if format == "gadf" and hdu: region_hdu = hdu + "_" + region_hdu if region_hdu in hdulist: region_table = Table.read(hdulist[region_hdu]) parser = FITSRegionParser(region_table) pix_region = parser.shapes.to_regions() wcs = WcsGeom.from_header(region_table.meta).wcs regions = [] for reg in pix_region: regions.append(reg.to_sky(wcs)) region = list_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): """Plot region in the sky. Parameters ---------- ax : `~astropy.vizualisation.WCSAxes` Axes to plot on. If no axes are given, the region is shown using the minimal equivalent WCS geometry. **kwargs : dict Keyword arguments forwarded to `~regions.PixelRegion.as_artist` """ import matplotlib.pyplot as plt from matplotlib.collections import PatchCollection from astropy.visualization.wcsaxes import WCSAxes if ax is None: ax = plt.gca() if not isinstance(ax, WCSAxes): ax.remove() wcs_geom = self.to_wcs_geom() m = Map.from_geom(wcs_geom.to_image()) fig, ax, cbar = m.plot(add_cbar=False) regions = compound_region_to_list(self.region) artists = [region.to_pixel(wcs=ax.wcs).as_artist() for region in regions] kwargs.setdefault("fc", "None") patches = PatchCollection(artists, **kwargs) ax.add_collection(patches) return ax