PIG 10 - Regions

• Author: Christoph Deil, Axel Donath, Régis Terrier

• Created: May 3, 2019

• Accepted: Jun 17, 2019

• Status: accepted

• Discussion: GH 2129

Abstract

We propose to use astropy-regions to handle spatial sky and pixel regions throughout Gammapy. We already use astropy-regions, so why a PIG? There are a few decisions we need to make concerning where to use sky and where to use pixel regions and where to support both. This affects the algorithms used (e.g. for rotated region background estimation) and the API (where a WCS or exclusion map is needed in functions and methods working with regions). Also astropy-regions was started after Gammapy - so we have some code to work with sky cones and boxes that should partially be removed, partially refactored to use astropy-regions.

The scope of this PIG is small and limited to spatial sky and pixel regions. The work to implement it is small, we propose to complete it for Gammapy v0.13 in July 2019.

The question of more general data subspace selections that include energy, time or phase regions and selections, or general n-dimensional regions, or whether to introduce a field of view (FOV) coordinate frame and special FOV regions is not addressed here. Also some things like the clean-up of gammapy.image which contains some region-related code is left as future work.

Introduction

This long introduction describes the history and design of astropy-regions and regions in Gammapy. It is background information, you are welcome to skip ahead to the Proposal and Task list in this PIG below.

Spatial regions are used a lot in gamma-ray astronomy (and other domains of astronomy as well). Often a users chooses a region of interest, and then selects the observations or events or pixels for analysis. Besides “on regions” or “regions of interest”, there are a few use cases to create regions within Gammapy, e.g. to restrict analysis to a maximum field of view offset, or to compute “off” regions for background estimation for a given “on” region.

In early versions of Gammapy we started to develop gammapy.regions, but soon realised that the functionality really is very generic, and also needed by many other astronomers (radio, optical, X-ray, …). Astropy and photutils was being developed, there was pyregion. At Python in Astronomy 2015 we discussed and decided to create a new Astropy-affiliated regions package from scratch and have been developing it since: astropy-regions. The goal is to merge it as astropy.regions into the Astropy core once it’s in a certain sense complete, for now it is an extra package that has been a dependency of Gammapy since Gammapy v0.5 in 2016. This is OK, the situation for astropy-regions is the same as for reproject and soon astropy-healpix (see GH 1167), they are extra packages and dependencies of Gammapy. We should continue to contribute to these three projects until they are complete and help get them in Astropy core some day.

The astropy-regions package contains sky and pixel regions with methods to filter positions and to make masks, ways to transform between sky and pixel region for a given WCS, support for plotting with matplotlib and I/O (string serialisation and parsing) for the DS9 region format. As described below, there are a few things still missing in astropy-regions that we need in Gammapy, thus the work described in this PIG is partly for astropy-regions and partly in Gammapy.

The places where regions appear in the Gammapy code and API are very few and overall region handling is pretty simple. Gammapy is mostly a binned analysis framework, so mainly regions are used to create pixel masks, and then the analysis is mask-based, not involving region objects any more (see section below: Images and masks). The second place where regions appear is to select objects of interest (events, observations, sources) for a given region (see section below: Sky and pixel regions). Finally, there is the use case of rotated regions for background estimation in 1D spectral analysis (see section below: Rotated regions).

However, it turns out that working with regions is surprisingly tricky, because there are different ways to define the basic semantics and computations for “contains” of a region, how to transform between pixel and sky regions and how to visualise them. The origin of the issue is that to make an image for analysis or visualisation, one has to choose a projection that maps sky to pixel coordinates, and there are many different projections, and they all have distortions. For example, the Wikipedia page on Aitoff projection contains an illustration of how sky circles map to pixel regions that are to a good approximation pixel circles in some parts of the image, but in others have complex shapes and different sizes, and at the edge even split in two disjoint pixel regions on the left and right part of the image.

In astropy-regions, after a lot of discussions, we have agreed to keep it simple, and to follow the lead of DS9 and to handle sky and pixel regions and their transforms in the same way that DS9 does. This has the advantage that we can use the DS9 string format to define and exchange regions (see section below: Region arguments), and we can use DS9 to view them. As can be seen in the CircleSkyRegion.to_pixel method implementation below, the definition used to map a sky circle to a pixel circle region is an approximation where the center is transformed according to the given wcs, and the radius is transformed using the local pixel scale at the center point for the given wcs.

class CircleSkyRegion(SkyRegion):
    def to_pixel(self, wcs):
        center, scale, _ = skycoord_to_pixel_scale_angle(self.center, wcs)
        radius = self.radius.to('deg').value * scale
        return CirclePixelRegion(center, radius, self.meta, self.visual)

So a sky circle is converted to a pixel circle. This works well for most use cases where one has an image with an extent and WCS of roughly constant pixel scale, e.g. a local TAN projection at a source location and 10 deg image size. In general the strategy in Gammapy will be to accept both pixel or sky regions as inputs, but to internally always call to_pixel first if a sky region is passed, and to internally always work with pixel regions.

This approximation fails for all-sky analysis and image parts where the projection has a shear or is varying like the sky circles that are distorted in the image on the Wikipedia page on Aitoff projection. One possible future improvement to astropy-regions would be to implement sky contains directly for CircleSkyRegion using sky distance, and for PolygonSkyRegion using the assumption that polygon edges follow great circles, and for other shapes to implement to_polygon methods that approximate any region as a PolygonSkyRegion (using any number of points and accuracy the user wants), and then transforming the PolygonSkyRegion to PolygonPixelRegion, as defined already by simply transforming the vertex points. This is feasible and useful for circles, ellipses and polygons (one would have to decide what to do with image projection edges like at longitude 180 deg in the Aitoff example, whether to split into disjoint pixel region parts), but is difficult for other shapes such as wedges or already for rectangles, because their shape only becomes well-defined for a given projection. These things are not in scope for this PIG or needed for Gammapy for now, Gammapy simply doesn’t have all-sky analysis support built in yet, although already now users could write Python scripts to do it using a combination of HEALPix and WCS maps and Gammapy.

Now that we know how astropy-regions works and how we plan to handle pixel and sky regions in Gammapy, let’s move on to the detailed proposal and task list outlining additions and changes to make.

Proposal

This section contains a detailed proposal concerning region-related code in Gammapy. The Task list at the end of this PIG references these descriptions.

Region arguments

We propose to add support for DS9 region strings (in addition to region objects) in all functions in Gammapy that take a region argument. This can be achieved via a make_region helper function:

>>> from regions import DS9Parser
>>> def make_region(region):
...     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()
...     else:
...         return region
...
>>> make_region("image;circle(42,43,3)")
<CirclePixelRegion(center=PixCoord(x=41.0, y=42.0), radius=3.0)>
>>> make_region("galactic;circle(42,43,3)")
<CircleSkyRegion(center=<SkyCoord (Galactic): (l, b) in deg
    (42., 43.)>, radius=3.0 deg)>

This is convenient, a simple region="galactic;circle(42,43,3)" is much shorter than having to remember and type this equivalent code:

>>> from astropy.coordinates import SkyCoord, Angle
>>> from regions import CircleSkyRegion
>>> center = SkyCoord(42, 43, unit="deg", frame="galactic")
>>> radius = Angle(3, "deg")
>>> region = CircleSkyRegion(center, radius)
>>> region
<CircleSkyRegion(center=<SkyCoord (Galactic): (l, b) in deg
    (42., 43.)>, radius=3.0 deg)>

There is precedence for this pattern in Gammapy, we usually pass Quantity and Angle arguments through via e.g. angle = Angle(angle) in Gammapy functions, and often call them with strings angle="3 deg".

We propose to add a second function make_pixel_region which in addition to the string case handling also does the sky region case handling, and gives good error messages on invalid inputs.

from regions import DS9Parser, SkyRegion, PixelRegion

def make_pixel_region(region, wcs=None):
    if isinstance(region, str):
        region = make_region(region)

    if isinstance(region, SkyRegion):
        if wcs is None:
            raise ValueError("Sky regions not supported here.")
        return region.to_pixel(wcs)
    elif isinstance(region, PixelRegion):
        return region
    else:
        raise TypeError("What is this? Giving up.")

Since in Gammapy we implement almost all algorithms using pixel regions, we could then call this helper function on each Gammapy method or function on the first line. E.g. gammapy.maps.WcsGeom.region_mask would start like this:

def region_mask(self, region):
    region = make_pixel_region(region)
    # do computation, knowing region is a PixelRegion object

We suggest to start by implementing these two helper functions in gammapy/utils/regions.py and to use it internally, but not expose it as part of the public API via the Gammapy docs. Once this is done, we would open an issue in astropy-regions asking whether inclusion there is welcome. This might or might not be the case, because make_region is basically a one-liner, and is DS9 format specific, whereas astropy-regions supports other region formats as well. Based on the outcome of this discussion, we would then either move gammapy/utils/regions.py to astropy-regions, or expose it as part of the public Gammapy API in the docs, since users might want to call them directly in their scripts or analysis codes that build on Gammapy.

We note that this means that we commit to support DS9 region strings as part of the Gammapy API. Many users will use this feature, so changing to something else would break user scripts. If we add this now, there are still a few months to gather user feedback, and to possibly add as similar make_sky_coord helper function throughout Gammapy and to consider how the format there compares with the regions format. There is also the HEALPix region string format in the HEALPix map spec and support in gammapy.maps; ideally this format wouldn’t exist, but avoiding it probably isn’t possible, since it includes the DISK_INC and HPX_PIXEL regions cannot be represented in the DS9 or any other existing region format. So at this time, there are no plans to improve or unify that region format.

We propose to change the current gammapy.maps.WcsGeom.region_mask input (and everywhere else) from a list of regions, to a single region. The same effect of a list of regions (meaning “union”) can be achieved by creating a Compound region. That should be as efficient computationally, and more flexible (e.g. one could also do “intersection” in the region definition). It might be useful to add a factory function in astropy-regions that makes the creation of such compound regions a bit simpler.

Sky and pixel regions

As already explained in the Introduction above, there are projection-related issues when it comes to defining the basic operations for sky regions like contains, to_pixel or plot. Currently contains and plot methods simply aren’t defined for SkyRegion, and for to_pixel only the local pixel scale approximation method is implemented.

This could be improved in astropy-regions:

• Implement contains for PolygonSkyRegion (see astropy-regions GH 46)

• Is is possible to define RectangleSkyRegion.contains without a WCS (i.e. to compute corners from data members)?

• Add SkyRegion.plot (see astropy-regions GH 221)?

• Add generic SkyRegion.contains using a local TAN or CAR projection with reference point at the region center, to allow use without a WCS?

• Add SkyRegion.to_polygon, try to be consistent with plot and contains.

Having improved support for sky regions would make it convenient to do some things without requiring a WCS. Especially for circle and polygon, it can be useful to filter positions from an observation or event list, and it is possible and common to do this. In Gammapy we currently have select_sky_box, select_sky_circle and skycoord_from_table which do this in an ad-hoc way. They are called from gammapy.data.EventList, gammapy.data.ObservationTable and gammapy.catalog.SourceCatalog which each are a wrapper or subclass of an astropy.table.Table.

We propose to replace those functions and methods with a select_region method, which for now requires a WCS to be passed on call.

In the future, if astropy-healpix sky region methods improve, passing a WCS might no longer be needed. Depending on the details of how SkyRegion.contains is defined, this change might or might not be backward compatible for Gammapy users. So ideally, this work should be done soon, but it is better discussed in the astropy-healpix issue tracker, not this Gammapy PIG.

Rotated regions

Rotated regions are used for background estimation in 1D spectral analysis in gamma-ray astronomy (see [Berge2007]). The rotation is done around the pointing position of the IACT telescope array. This method is also called “reflected regions”, but really that’s a bad name, reflecting would only give one region on the other side; what is done is rotating the regions.

In Gammapy, rotated regions are implemented in gammapy/background/reflected.py. This is one of the most complex pieces of code in Gammapy, because it was started before astropy.coordinates or astropy-regions existed and then adapted many times, never achieving a clean design and implementation until now. The problem was that astropy-region objects didn’t support a rotate operation, and also their representation is such that it’s not clear how to compute and represent e.g. a rotated RectangleSkyRegion (see discussion of ambiguity in SkyRegion.contains above). Also, the current code tried to handle the case with and without having an exclusion mask and a WCS.

A recent attempt to improve our rotated regions code is GH 2089, which shows that it is difficult to to in a clean way without having a Region.rotate method, because if-else and specific case handling is needed for the operation.

We propose to add a PixelRegion.rotate(center, angle) method for all regions in astropy-regions (or all that can be represented again as a pixel region object after rotating), and to only handle the case of pixel regions and a given exclusion mask and WCS in reflected.py. For convenience, a dummy WCS as a local TAN projection to the pointing position is used, if the user doesn’t need of have an exclusion mask, as is already done in the current code. This is a common case e.g. for isolated point sources like AGN, where the rotated regions method is most commonly used.

One complication is that in the current code, to have fewer off test regions, the algorithm takes a large step arctan(2*radius/offset) at the start and after placing an off region. This requires a “radius”, which is is not available for regions of a complex shape. We propose that this optimisation is only left for the circle case (which already has a specialised fast algorithm anyways using a distance mask), and for other regions for now the normal thing is done to test all positions in small angle increments.

Images and masks

Often regions are used to analyse images, and usually this happens via binary masks (with False for “outside” and True for “inside”), or in some cases integer label images (with 0 for “outside”, 1 for “in region 1”, 2 for “in region 2” and N for “in region N”). This is used in Gammapy, but also in scipy.ndimage, scikit-image, photutils and other packages.

We leave the description of the details how we work with images and masks in Gammapy as future work. This could be done in another PIG, or just directly by improving the Gammapy code and documentation.

There are open issues already discussing when to use bool or int arrays, and a lot of code in gammapy.image (especially gammapy/image/measure.py and gammapy/image/profile.py) should either be removed from Gammapy, replaced by documentation examples, or improved to in some cases used regions. This is not high priority and a separate task, since this functionality isn’t used anywhere else in Gammapy.

Outlook

This PIG only describes some short-term improvements to spatial regions. Long-term, we might want to develop a more general solution for regions that also includes non-spatial regions, e.g. spectral regions or time regions.

Just to mention what we have at the moment: spectral “regions” are handled via energy_bounds or (energy_min, energy_max) in a not very consistent way. We could look towards the Astropy specutils package and the specutils.SpectralRegion object as an example how to implement spectral regions. Currently there is plan to integrate astropy-regions and specutils, astropy-regions GH 228 has been opened for discussion.

Supporting an energy-dependent maximum field of view radius selection, or a field-of-view dependent safe energy offset might be a motivation to couple spatial and spectral region, possibly supported via gammapy.maps.

For time regions we currently often have (time_min, time_max) in the API, and we have the gammapy.data.GTI class to represent multiple intervals. It’s not clear if this needs to be changed or wrapped somehow into a “time region”, or if there are any use cases to couple time with the spatial and spectral regions, e.g. for a time-dependent safe energy threshold or user analysis option.

There is an issue GH 1805 “Introduce RegionGeom and RegionNDMap?” with first thoughts in this direction. Discussion or working out a solution for that is beyond the scope of this PIG.

In the HESS software, there is a region class that wraps a mask, and that is sometimes useful, specifically because in the HESS software exclusion regions are region objects, sometimes compound with a mask region and other regions like circles. It’s clearly an option, easy to implement as a SkyRegion sub-class. But also causes problems, especially when it comes to serialisation. Thoughts? I’m OK to not add this for now, and just use Map objects for exclusion masks, and require that users create that before using Gammapy, in case they e.g. have a list of circles or whatever they want to use as exclusion region.

Task list

We propose to implement this PIG though a series of small pull requests. Some are dependent on others, but for many also the order doesn’t matter and overall the amount of work is limited.

• Add region copy method, needed for region rotate (see astropy-regions GH 266)

• Add region rotate method (see astropy-regions GH 265 and Rotated regions above)

• Release astropy-regions v0.4, then bump Gammapy dependency to that version.

• Improve rotated region code and tests and docs in Gammapy. This could mean a rebase and continuation of GH 2089, or closing that and re-starting from master, copying over the useful parts or cherry-picking useful commits from that branch.

Other tasks, mostly independent from the previous ones, can happen in any order:

• Add gammapy/utils/regions.py with two helper functions. See Region arguments above.

• Use the helper functions throughout Gammapy. Should be 5-10 cases, could be several PRs. See Region arguments above.

• Replace sky circle and box select with generic region select (see Sky and pixel regions, resolves GH 1172)

Beyond those core tasks that really should be done for Gammapy v0.13, further improvements in astropy-healpix and Gammapy should be done, but are lower priority:

• Implement SkyRegion.contains(point) (see Sky and pixel regions above)

• Implement SkyRegion.plot (see Sky and pixel regions above)

• Review and improve region-related tests, especially for complicated edge lon = 0 / 360 and the poles.

• Review and improve region-related end-user docs. Possibly add a dedicated tutorial notebook.

• “Support region_mask for multi-resolution maps” (GH 1715)

• “HPX up/down sample issue with partial skymaps” (GH 1445)

• “Introduce RegionGeom and RegionNDMap?” (GH 1805)

Alternatives

As mentioned in Sky and pixel regions and Rotated regions, we could use sky regions more instead of pixel regions, and restrict analysis to sky regions where contains or rotate is well defined, namely circles and polygons, and to polygonise other sky regions using some approximation, before passing them to Gammapy.

In the future we might want a more general solution, and there’s the question whether a mask-backed region should be added (see Outlook)

Decision

The PIG was discussed in GH 2129 and the weekly Gammapy developer calls. A final review announced on the Gammapy and CC mailing list did not yield any additional comments. Therefore the PIG was accepted on June 17, 2019.