Repository: gammapy/gammapy
Pull request #5715 (2025-03-04 08:33:30)

Modified Polygon region for use as ON region
<!-- These are hidden commments. You can leave them or delete them. -->

<!-- Don't be shy! Github issues or pull requests are welcome any time. -->
<!-- But if you first want to chat, use the Gammapy Slack. -->
<!-- See https://gammapy.org/contact.html -->

**Description**
<!-- What is the motivation for this pull request? -->
<!-- Briefly: what changes are done here? -->
<!-- If this is to address a Github issue, mention its number and use the specific keywords to create a link -->
<!-- See docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue -->
<!-- Example: This PR is to resolve #42 -->

This is a minimal reimplementation of the polygon region class to allow it to be used as a signal extraction region in analysis. The significant addition to this was the addition of a center class member which is missing in the regions implementation of this class. This centre is calculated as the geometric centre of the points when the class is created. The Polygon object was creasted as a sky region and a pixel region to match the regions class behaviour and the centre is propagated between them.

**Dear reviewer**
<!-- Let the reviewer and Gammapy team know what you want: -->
<!-- * Is this ready for review? Or is it work in progress and you want some feedback? -->
<!-- * Do you want to go through review here? Or if someone should just finish this up and merge it? -->
<!-- Do you have any specific questions, e.g. about API or implementation? -->
<!-- Do you include a test executing the new code you're adding (to make sure it runs)? -->
<!-- Do you include some documentation? Is it needed? -->
I created this as a draft to check that I was taking the correct approach here for implementation, also it is not clear to me how best to implement tests for this class.

QRemy: Thanks @ParsonsRD, maybe it should be added to regions rather than gammapy, opinions @adonath @registerrier ?
adonath: @QRemy I have discussed this with @ParsonsRD shortly in private and the problem with adding this to `regions` right now, is that there is no single agreed upon definition of what e.g. the center is, than there is the question of spherical polygons, which algorithm to use for "contains" etc. I think it is fine to provide a prototype / reference implementation in Gammapy and then strat discussions and finally move upstream. This is true for most functionality in `gammapy/utils/regions.py`. At the same time it provides a real benefit for Gammapy users already.
bkhelifi: I am wondering whether the exact definition of this region centre is very critical. How thought about that? This implementation looks promising to me...
In the mid-run, should we also propose a feature request to astropy?
adonath: > I am wondering whether the exact definition of this region centre is very critical. How thought about that? This implementation looks promising to me...

I think it is mostly a matter of convention. E.g. I have introduce a helper method to compute the "center" of a set of points via it's geometric median: https://github.com/gammapy/gammapy/blob/main/gammapy/utils/regions.py#L38

Not sure whether this give the same results as what @ParsonsRD introduces in this PR. We do not necessarily have to strive for consistency, but be should be ware that there are different possible approaches. 

registerrier: Thanks @ParsonsRD . This is interesting. To proceed this would require tests to be implemented.
Can we compare the center found in simple cases with expected values?
ParsonsRD: Indeed this was what I was thinking, to create a region that is a equivalent to another implemented region shape and check we get the same result. Do you know if there is an example of a test like this that I could try to reproduce?
QRemy: Added AI generated unit tests.

Patch:
diff --git a/gammapy/utils/regions.py b/gammapy/utils/regions.py
index f68d8e388..92cca0d83 100644
--- a/gammapy/utils/regions.py
+++ b/gammapy/utils/regions.py
@@ -12,6 +12,7 @@ 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
 import numpy as np
 from scipy.optimize import Bounds, minimize
@@ -24,8 +25,14 @@ from regions import (
     EllipseSkyRegion,
     RectangleSkyRegion,
     Regions,
+    PolygonSkyRegion,
+    PolygonPixelRegion,
 )
 
+from regions.core.pixcoord import PixCoord
+from regions.core.metadata import RegionMeta, RegionVisual
+from regions._utils.wcs_helpers import pixel_scale_angle_at_skycoord
+
 __all__ = [
     "compound_region_to_regions",
     "make_concentric_annulus_sky_regions",
@@ -137,6 +144,44 @@ def regions_to_compound_region(regions):
     return region_union
 
 
+def get_centroid(vertices):
+    """Compute centroid of a polygon. Implicitly assumes a flat
+    cartesian projection, will probably break for very large polygons.
+
+    Code comes from:
+    https://stackoverflow.com/questions/75699024/finding-the-centroid-of-a-polygon-in-python
+
+    Parameters
+    ----------
+    vertices : `~astropy.coordinates.SkyCoord`
+        List of vertices.
+
+    Returns
+    -------
+    centroid : `~astropy.coordinates.SkyCoord`
+        Centroid of the polygon.
+    """
+    polygon = []
+    for vertex in vertices:
+        polygon.append((vertex.ra.degree, vertex.dec.degree))
+    polygon = np.array(polygon)
+
+    # Same polygon, but with vertices cycled around. Now the polygon
+    # decomposes into triangles of the form origin-polygon[i]-polygon2[i]
+    polygon2 = np.roll(polygon, -1, axis=0)
+
+    # Compute signed area of each triangle
+    signed_areas = 0.5 * np.cross(polygon, polygon2)
+
+    # Compute centroid of each triangle
+    centroids = (polygon + polygon2) / 3.0
+
+    # Get average of those centroids, weighted by the signed areas.
+    centroid = np.average(centroids, axis=0, weights=signed_areas)
+
+    return SkyCoord(centroid[0] * u.deg, centroid[1] * u.deg, frame=vertices.frame)
+
+
 class SphericalCircleSkyRegion(CircleSkyRegion):
     """Spherical circle sky region.
 
@@ -157,6 +202,109 @@ class SphericalCircleSkyRegion(CircleSkyRegion):
         return separation < self.radius
 
 
+class PolygonPointsSkyRegion(PolygonSkyRegion):
+    """Polygon sky region defined by a list of points."""
+
+    def __init__(self, vertices, meta=None, visual=None):
+        """Create a polygon sky region.
+
+        Parameters
+        ----------
+        vertices : `~astropy.coordinates.SkyCoord`
+            List of vertices.
+        meta : `~regions.RegionMeta`, optional
+            Region meta data.
+        visual : `~regions.RegionVisual`, optional
+            Region visual meta data.
+        """
+        self.vertices = vertices
+        self.meta = meta or RegionMeta()
+        self.center = get_centroid(vertices)
+        self.visual = visual or RegionVisual()
+
+    def to_pixel(self, wcs):
+        """Convert to pixel region."""
+        x, y = wcs.world_to_pixel(self.vertices)
+        center, _, _ = pixel_scale_angle_at_skycoord(self.center, wcs)
+
+        vertices_pix = PixCoord(x, y)
+        return PolygonPointsPixelRegion(
+            vertices_pix,
+            center=center,
+            meta=self.meta.copy(),
+            visual=self.visual.copy(),
+        )
+
+
+class PolygonPointsPixelRegion(PolygonPixelRegion):
+    """Polygon pixel region defined by a list of points."""
+
+    def __init__(self, vertices, center=None, meta=None, visual=None, origin=None):
+        """Create a polygon pixel region.
+
+        Parameters
+        ----------
+        vertices : `~regions.PixCoord`
+            List of vertices.
+        center : `~regions.PixCoord`, optional
+            Center of the region.
+        meta : `~regions.RegionMeta`, optional
+            Region meta data.
+        visual : `~regions.RegionVisual`, optional
+            Region visual meta data.
+        origin : `~regions.PixCoord`, optional
+            Origin of the region. Default is `PixCoord(0, 0)`
+        """
+
+        if origin is None:
+            origin = PixCoord(0, 0)
+
+        self._vertices = vertices
+        self.meta = meta or RegionMeta()
+        self.visual = visual or RegionVisual()
+        self.origin = origin
+        self.vertices = vertices + origin
+        self.center = center
+
+    def to_sky(self, wcs):
+        """Convert to sky region.
+
+        Parameters
+        ----------
+        wcs : `~astropy.wcs.WCS`
+            WCS transformation object.
+
+        """
+        vertices_sky = wcs.pixel_to_world(self.vertices.x, self.vertices.y)
+
+        return PolygonPointsSkyRegion(
+            vertices=vertices_sky, meta=self.meta.copy(), visual=self.visual.copy()
+        )
+
+    def rotate(self, center, angle):
+        """
+        Rotate the region.
+
+        Positive ``angle`` corresponds to counter-clockwise rotation.
+
+        Parameters
+        ----------
+        center : `~regions.PixCoord`
+            The rotation center point.
+        angle : `~astropy.coordinates.Angle`
+            The rotation angle.
+
+        Returns
+        -------
+        region : `PolygonPixelRegion`
+            The rotated region (which is an independent copy).
+        """
+        vertices = self.vertices.rotate(center, angle)
+        center = self.center.rotate(center, angle)
+
+        return self.copy(vertices=vertices, center=center)
+
+
 def make_orthogonal_rectangle_sky_regions(start_pos, end_pos, wcs, height, nbin=1):
     """Utility returning an array of regions to make orthogonal projections.
 


Test patch:
diff --git a/gammapy/utils/tests/test_regions.py b/gammapy/utils/tests/test_regions.py
index 14faafe40..b6467e585 100644
--- a/gammapy/utils/tests/test_regions.py
+++ b/gammapy/utils/tests/test_regions.py
@@ -7,16 +7,30 @@ in https://astropy-regions.readthedocs.io that we rely on in Gammapy.
 That package is still work in progress and not fully developed and
 stable, so need to establish a bit what works and what doesn't.
 """
+
+import numpy as np
 from numpy.testing import assert_allclose, assert_equal
 import astropy.units as u
-from astropy.coordinates import SkyCoord
-from regions import CircleSkyRegion, EllipseSkyRegion, Regions
+from astropy.coordinates import SkyCoord, Angle
+from regions import (
+    CircleSkyRegion,
+    EllipseSkyRegion,
+    Regions,
+    RegionMeta,
+    RegionVisual,
+    PixCoord,
+)
+from astropy.wcs import WCS
+
 from gammapy.utils.regions import (
     SphericalCircleSkyRegion,
     compound_region_center,
     region_circle_to_ellipse,
     region_to_frame,
     regions_to_compound_region,
+    get_centroid,
+    PolygonPointsSkyRegion,
+    PolygonPointsPixelRegion,
 )
 
 
@@ -102,3 +116,71 @@ def test_region_circle_to_ellipse():
     assert_allclose(region_new.height, region.radius, rtol=1e-3)
     assert_allclose(region_new.width, region.radius, rtol=1e-3)
     assert_allclose(region_new.angle, 0.0 * u.deg, rtol=1e-3)
+
+
+def test_get_centroid():
+    vertices = SkyCoord(
+        [
+            (0, 0),
+            (0, 1),
+            (1, 1),
+            (1, 0),
+        ],
+        unit="deg",
+    )
+
+    expected_centroid = SkyCoord(0.5 * u.deg, 0.5 * u.deg)
+
+    calculated_centroid = get_centroid(vertices)
+
+    assert_allclose(calculated_centroid.ra.deg, expected_centroid.ra.deg, rtol=1e-7)
+    assert_allclose(calculated_centroid.dec.deg, expected_centroid.dec.deg, rtol=1e-7)
+
+
+def test_polygon_points_sky_region_init():
+    vertices = SkyCoord([(0, 0), (0, 1), (1, 1), (1, 0)], unit="deg")
+    region = PolygonPointsSkyRegion(vertices)
+    assert np.all(region.vertices == vertices)
+    assert isinstance(region.meta, RegionMeta)
+    assert isinstance(region.visual, RegionVisual)
+    expected_centroid = SkyCoord(0.5 * u.deg, 0.5 * u.deg)
+    assert_allclose(region.center.ra.deg, expected_centroid.ra.deg, rtol=1e-7)
+    assert_allclose(region.center.dec.deg, expected_centroid.dec.deg, rtol=1e-7)
+
+
+def test_polygon_points_sky_region_to_pixel_to_sky():
+    vertices = SkyCoord([(0, 0), (0, 1), (1, 1), (1, 0)], unit="deg")
+    region = PolygonPointsSkyRegion(vertices)
+    wcs = WCS(naxis=2)
+    wcs.wcs.crpix = [0, 0]
+    wcs.wcs.cdelt = [1, 1]
+    wcs.wcs.crval = [0, 0]
+    wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
+    pixel_region = region.to_pixel(wcs)
+    assert isinstance(pixel_region, PolygonPointsPixelRegion)
+
+    region_new = pixel_region.to_sky(wcs)
+
+    assert isinstance(region_new, PolygonPointsSkyRegion)
+    assert_allclose(region_new.vertices.ra[0], 0 * u.deg, atol=1e-7)
+    assert_allclose(region_new.vertices.dec[0], 0 * u.deg, atol=1e-7)
+
+
+def test_polygon_points_pixel_region_init():
+    vertices = PixCoord([0, 0], [1, 1])
+    region = PolygonPointsPixelRegion(vertices)
+    assert region.vertices.x[0] == 0
+    assert region.vertices.y[0] == 1
+    assert isinstance(region.meta, RegionMeta)
+    assert isinstance(region.visual, RegionVisual)
+
+
+def test_polygon_points_pixel_region_rotate():
+    vertices = PixCoord([0, 1], [1, 0])
+    center = PixCoord(0.5, 0.5)
+    region = PolygonPointsPixelRegion(vertices, center=center)
+    angle = Angle(90, unit="deg")
+    rotated_region = region.rotate(center, angle)
+    expected_vertices = PixCoord([0, 1], [0, 1])
+    assert_allclose(rotated_region.vertices.x, expected_vertices.x, atol=1e-7)
+    assert_allclose(rotated_region.vertices.y, expected_vertices.y, atol=1e-7)