Source code for gammapy.utils.regions

# Licensed under a 3-clause BSD style license - see LICENSE.rst
"""Regions helper functions.

Throughout Gammapy, we use `regions` to represent and work with regions.

https://astropy-regions.readthedocs.io

The functions ``make_region`` and ``make_pixel_region`` should be used
throughout Gammapy in all functions that take ``region`` objects as input.
They do conversion to a standard form, and some validation.

We might add in other conveniences and features here, e.g. sky coord contains
without a WCS (see "sky and pixel regions" in PIG 10), or some HEALPix integration.

TODO: before Gammapy v1.0, discuss what to do about ``gammapy.utils.regions``.
Options: keep as-is, hide from the docs, or to remove it completely
(if the functionality is available in ``astropy-regions`` directly.
"""
import operator
from regions import (
    CircleSkyRegion,
    CompoundSkyRegion,
    DS9Parser,
    PixelRegion,
    Region,
    SkyRegion,
)

__all__ = ["make_region", "make_pixel_region"]


[docs]def make_region(region): """Make region object (`regions.Region`). See also: * `gammapy.utils.regions.make_pixel_region` * https://astropy-regions.readthedocs.io/en/latest/ds9.html * http://ds9.si.edu/doc/ref/region.html Parameters ---------- region : `regions.Region` or str Region object or DS9 string representation Examples -------- If a region object in DS9 string format is given, the corresponding region object is created. Note that in the DS9 format "image" or "physical" coordinates start at 1, whereas `regions.PixCoord` starts at 0 (as does Python, Numpy, Astropy, Gammapy, ...). >>> from gammapy.utils.regions import make_region >>> make_region("image;circle(10,20,3)") <CirclePixelRegion(PixCoord(x=9.0, y=19.0), radius=3.0)> >>> make_region("galactic;circle(10,20,3)") <CircleSkyRegion(<SkyCoord (Galactic): (l, b) in deg (10., 20.)>, radius=3.0 deg)> If a region object is passed in, it is returned unchanged: >>> region = make_region("image;circle(10,20,3)") >>> region2 = make_region(region) >>> region is region2 True """ if isinstance(region, str): # This is basic and works for simple regions # It could be extended to cover more things, # like e.g. compound regions, exclusion regions, .... return DS9Parser(region).shapes[0].to_region() elif isinstance(region, Region): return region else: raise TypeError(f"Invalid type: {region!r}")
[docs]def make_pixel_region(region, wcs=None): """Make pixel region object (`regions.PixelRegion`). See also: `gammapy.utils.regions.make_region` Parameters ---------- region : `regions.Region` or str Region object or DS9 string representation wcs : `astropy.wcs.WCS` WCS Examples -------- >>> from gammapy.maps import WcsGeom >>> from gammapy.utils.regions import make_pixel_region >>> wcs = WcsGeom.create().wcs >>> region = make_pixel_region("galactic;circle(10,20,3)", wcs) >>> region <CirclePixelRegion(PixCoord(x=570.9301128316974, y=159.935542455567), radius=6.061376992149382)> """ if isinstance(region, str): region = make_region(region) if isinstance(region, SkyRegion): if wcs is None: raise ValueError("Need wcs to convert to pixel region") return region.to_pixel(wcs) elif isinstance(region, PixelRegion): return region else: raise TypeError(f"Invalid type: {region!r}")
def compound_region_to_list(region): """Create list of regions from compound regions. Parameters ---------- regions : `~regions.CompoundSkyRegion` or `~regions.SkyRegion` Compound sky region Returns ------- regions : list of `~regions.SkyRegion` List of regions. """ regions = [] if isinstance(region, CompoundSkyRegion): if region.operator is operator.or_: regions_1 = compound_region_to_list(region.region1) regions.extend(regions_1) regions_2 = compound_region_to_list(region.region2) regions.extend(regions_2) else: raise ValueError("Only union operator supported") else: return [region] return regions def list_to_compound_region(regions): """Create compound region from list of regions, by creating the union. Parameters ---------- regions : list of `~regions.SkyRegion` List of regions. Returns ------- compound : `~regions.CompoundSkyRegion` Compound sky region """ region_union = regions[0] for region in regions[1:]: region_union = region_union.union(region) return region_union class SphericalCircleSkyRegion(CircleSkyRegion): """Spherical circle sky region. TODO: is this separate class a good idea? - If yes, we could move it to the ``regions`` package? - If no, we should implement some other solution. Probably the alternative is to add extra methods to the ``CircleSkyRegion`` class and have that support both planar approximation and spherical case? Or we go with the approach to always make a TAN WCS and not have true cone select at all? """ def contains(self, skycoord, wcs=None): """Defined by spherical distance.""" separation = self.center.separation(skycoord) return separation < self.radius