Introduce SkyRegion.solid_angle

[adonath]

The PixelRegion already defines an .area property. It would be nice to see the equivalent SkyRegion.solid_angle property.

[adonath]

I guess technically this might not be simple to achieve. While of course for some regions once can rely on analytical expressions, one would still need a numerical approximations in other situations. In any case https://mathworld.wolfram.com/LHuiliersTheorem.html is useful...

[larrybradley]

I'm not sure how this would work because the sky-based regions are simple representations that preserve the shape. They aren't shapes in spherical coordinates, but more like cartesian shapes in a tangent plane (a circle in pixel coordinates converts to a circle in sky coordinates based only on the pixel scale). Does a solid angle make sense for non-spherical coordinates?

[adonath]

Thanks for the feedback @larrybradley! Yes, I now remember that methods like SkyRegion.contains() actually require a wcs to handle the projection to the plane. However I think .solid_angle() could be maybe supported in a similar way?

• E.g., compute the pixel area and then multiply with the local approximate wcs bin size, that is used for the sky to pix conversion
• Or maybe compute an array of solid angle per pixel and multiply with PixRegion.to_mask() similar to:

def solid_angle(self, wcs, mode):
    """Solid angle enclosed by the region (`astropy.units.Quantity`)"""
    region_pix = self.to_pix(wcs=wcs)
    mask = region_pix.to_mask(mode="exact")
    
    solid_angle_array = compute_solid_angle_array_from_wcs(wcs, shape=mask.data.shape)
    
    return np.sum(mask.data * solid_angle_array)

[tjgalvin]

I would be very curious to see how this progresses. Is there a currently utility to compute the area of a pixel for some arbitary position anywhere? I know of something in astropy wcs, but that performs this at the reference pixel exclusively. I'd would try to help out but not really sure where to start.

[adonath]

@tjgalvin We have something like this in gammapy.maps:

from gammapy.maps import WcsGeom
from astropy import units as u
from astropy.coordinates import SkyCoord

center = SkyCoord("0d", "0d", frame="galactic")

geom = WcsGeom.create(
    skydir=center,
    binsz=0.1,
    width=[10, 8] * u.deg
)
geom.solid_angle()

Which returns an array of solid angle per pixel, computed using L'Huilllers theorem.