.. include:: ../../references.txt .. _pig-010: **************** 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)") >>> make_region("galactic;circle(42,43,3)") , 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 , 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. .. These explicit URLs to Gammapy classes are to avoid possible future breakage of .. the links in the PIG if those classes are removed or changed: .. _gammapy.data.GTI: https://docs.gammapy.org/0.11/api/gammapy.data.GTI.html .. _gammapy.maps.WcsGeom.region_mask: https://docs.gammapy.org/0.11/api/gammapy.maps.WcsGeom.html#gammapy.maps.WcsGeom.region_mask .. _gammapy.data.EventList: https://docs.gammapy.org/0.11/api/gammapy.data.EventList.html .. _gammapy.data.ObservationTable: https://docs.gammapy.org/0.11/api/gammapy.data.ObservationTable.html .. _gammapy.catalog.SourceCatalog: https://docs.gammapy.org/0.11/api/gammapy.catalog.SourceCatalog.html .. _select_sky_box: https://docs.gammapy.org/0.11/api/gammapy.catalog.select_sky_box.html .. _select_sky_circle: https://docs.gammapy.org/0.11/api/gammapy.catalog.select_sky_circle.html .. _skycoord_from_table: https://docs.gammapy.org/0.11/api/gammapy.catalog.skycoord_from_table.html .. _gammapy.image: https://docs.gammapy.org/0.11/image/index.html#reference-api .. _Wikipedia page on Aitoff projection: https://en.wikipedia.org/wiki/Aitoff_projection .. _HEALPix region string: https://gamma-astro-data-formats.readthedocs.io/en/latest/skymaps/healpix/index.html#healpix-region-string .. _photutils: https://photutils.readthedocs.io .. _Python in Astronomy 2015: http://openastronomy.org/pyastro/2015/ .. _astropy-regions: https://astropy-regions.readthedocs.io .. _pyregion: https://pyregion.readthedocs.io .. _DS9 region format: http://ds9.si.edu/doc/ref/region.html .. _reproject: https://reproject.readthedocs.io .. _astropy-healpix: https://astropy-healpix.readthedocs.io .. _spherical_geometry: https://github.com/spacetelescope/spherical_geometry .. _regions.PolygonSkyRegion: https://astropy-regions.readthedocs.io/en/latest/api/regions.PolygonSkyRegion.html .. _specutils: https://specutils.readthedocs.io .. _specutils.SpectralRegion: https://specutils.readthedocs.io/en/latest/spectral_regions.html .. _Compound region: https://astropy-regions.readthedocs.io/en/latest/compound.html .. _GH 2129: https://github.com/gammapy/gammapy/pull/2129 .. _GH 2089: https://github.com/gammapy/gammapy/pull/2089 .. _GH 1172: https://github.com/gammapy/gammapy/issues/1172 .. _GH 1167: https://github.com/gammapy/gammapy/pull/1167 .. _GH 1715: https://github.com/gammapy/gammapy/issues/1715 .. _GH 2068: https://github.com/gammapy/gammapy/issues/2068 .. _GH 1445: https://github.com/gammapy/gammapy/issues/1445 .. _GH 1805: https://github.com/gammapy/gammapy/issues/1805 .. _astropy-regions GH 46: https://github.com/astropy/regions/issues/46 .. _astropy-regions GH 91: https://github.com/astropy/regions/pull/91 .. _astropy-regions GH 217: https://github.com/astropy/regions/issues/217 .. _astropy-regions GH 221: https://github.com/astropy/regions/pull/221 .. _astropy-regions GH 228: https://github.com/astropy/regions/issues/228 .. _astropy-regions GH 265: https://github.com/astropy/regions/issues/265 .. _astropy-regions GH 266: https://github.com/astropy/regions/issues/266