compute_pixel_area_map#

photutils.utils.compute_pixel_area_map(wcs, shape, *, interpolate=True, step=None)[source]#

Compute the on-sky area of every pixel in an image.

By default, the areas are computed with compute_pixel_areas on a coarse grid of positions and then interpolated onto the full image grid with a bicubic spline. Distortions vary smoothly on scales far larger than the grid spacing, so the interpolation error is negligible compared with the variation in area that the map captures. Set interpolate=False to instead evaluate the area directly at every pixel.

Parameters:
wcsWCS object

A world coordinate system (WCS) transformation that supports the astropy shared interface for WCS (e.g., astropy.wcs.WCS, gwcs.wcs.WCS).

shape2-tuple of int

The (ny, nx) shape of the image.

interpolatebool, optional

Whether to interpolate the areas from a coarse grid. If False, the area is evaluated directly at every pixel with compute_pixel_areas. That is exact but slower (a few seconds for a 4096 x 4096 image). For any image more than a few dozen pixels on a side, interpolation is faster than direct evaluation at every step of 2 or more (about 3 times at a step of 2 and more than 10 times at a step of 4), so a small step is not a reason to disable it. Set interpolate=False only when the exact per-pixel areas are wanted.

stepint or None, optional

The spacing in pixels of the coarse grid on which the areas are evaluated. The grid extends two steps before the first pixel and at least one step beyond the last pixel, so the interpolation never extrapolates and the spline is well defined for any image shape. Steps from about 32 to 128 give equivalent results. Smaller steps only add runtime and larger steps lose accuracy on strongly distorted or very wide fields. The largest allowed step is min(shape) // 8 (and at least 1), so that the coarse grid always has at least eight intervals across the image. If None, the step is 64 or the largest allowed step, whichever is smaller. A larger step than allowed is reduced to the largest allowed step with a warning. A step of 1 samples every pixel, so the area is then evaluated directly at every pixel without interpolation, as if interpolate were False. step must be None if interpolate is False.

Returns:
areasndarray

The 2D array of pixel areas in arcsec2, with the same shape as the image. The area is NaN for pixels where the WCS is undefined (e.g., outside the valid region of the projection).

Notes

If any node of the padded coarse grid lies outside the valid region of the projection, which can happen for an all-sky map even when every image pixel is valid, the areas are instead evaluated directly at every pixel with compute_pixel_areas, as if interpolate were False.

Examples

>>> import numpy as np
>>> from photutils.datasets import make_wcs
>>> from photutils.utils import compute_pixel_area_map
>>> shape = (100, 100)
>>> wcs = make_wcs(shape)
>>> areas = compute_pixel_area_map(wcs, shape)
>>> areas.shape
(100, 100)
>>> print(np.round(areas[50, 50], 6))  # arcsec**2
0.01