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(PixCoord(x=41.0, y=42.0), radius=3.0)>
>>> make_region("galactic;circle(42,43,3)")
<CircleSkyRegion(<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(<SkyCoord (Galactic): (l, b) in deg
(42., 43.)>, radius=3.0 deg)>
There is precendence 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 arount 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)