Cutout2D behaviour for non-rectangular pixel sizes

[adonath]

Description

When specifying a non-rectangular pixel size for the WCS provided to the Cutout2D class, the cutout does not return the expected cutout size and data shape.

Steps to Reproduce

from astropy.nddata.utils import Cutout2D
from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales
from astropy.coordinates import SkyCoord
import numpy as np

center = SkyCoord("0d", "0d", frame="icrs")
wcs = celestial_frame_to_wcs(center.frame, projection="CAR")

wcs.wcs.crval = (0, 0)
wcs.wcs.cdelt = (-5, 5)
wcs.wcs.crpix = (36.5, 18.5)

data = np.ones((36, 72))
print(data.shape)

size = [90, 180] * u.deg
c2d = Cutout2D(data=data, position=center, size=size, wcs=wcs)
print(c2d.data.shape)

Prints:

(36, 72)
(18, 36)

So it behaves as expected and returns a data array with half the input shape in all axes.

In case of specifying an un-even pixel size:

from astropy.nddata.utils import Cutout2D
from astropy.wcs.utils import celestial_frame_to_wcs, proj_plane_pixel_scales
from astropy.coordinates import SkyCoord
import numpy as np

center = SkyCoord("0d", "0d", frame="icrs")
wcs = celestial_frame_to_wcs(center.frame, projection="CAR")

wcs.wcs.crval = (0, 0)
wcs.wcs.cdelt = (-10, 5)
wcs.wcs.crpix = (18.5, 18.5)

data = np.ones((36, 36))
print(data.shape)

size = [90, 180] * u.deg
c2d = Cutout2D(data=data, position=center, size=size, wcs=wcs)
print(c2d.data.shape)

It prints:

(36, 36)
(9, 36)

While I would have expected (18, 18). Maybe it seems the lon and lat axes are switched in this case? cc @LauraOlivera

System Details

macOS-10.14.6-x86_64-i386-64bit
Python 3.8.6 | packaged by conda-forge | (default, Jan 25 2021, 23:22:12) 
[Clang 11.0.1 ]
Numpy 1.20.0
astropy 4.1
Scipy 1.6.0
Matplotlib 3.3.4

[adonath]

From a quick look at the code, I think I already found the bug:

• The size argument of the Cutout2D class is given in (lat, lon) order, as documented in the docstring
• When converting the size to pixel in the loop here proj_plane_pixel_scales(wcs) returns the pixel size in (lon, lat) order
• So it seems like the incorrect pixel size is used to convert the size to pixels...lon for lat and the other way around.

If someone confirms the bug, I could possibly implement a fix myself.

[pllim]

Feel free to open an exploratory PR with a test added for this use case. If the CI passes (minus the unrelated failures), then you have proof. Thanks! 😄

[mcara]

I am not sure this is a bug. proj_plane_pixel_scales(wcs) is supposed to return something similar to cdelt (identical to cdelt when PC matrix is identity).

For your WCS, proj_plane_pixel_scales(wcs) recovers the cdelt (in absolute value) that you have specified:

>>> proj_plane_pixel_scales(wcs)
array([10.,  5.])

Let's look at the WCS itself. You said that the pixel scale along the first image coordinate is 10 deg and along the second image coordinate is 5 deg. Given the shape of the extraction box is [90, 180] * u.deg, I do not find it surprising that it gets translated to image space as 90/10=9 and 180/5=36. I may be wrong (I am not expert in carree projection).

[adonath]

@mcara The point is just, that the size parameter is given in (lat, lon) order in the Cutout2D class. So the size of the extraction box is defined as 90 deg in lat and 180 deg in lon, while the pixel size is 10 deg in lon and 5 deg in lat. So in your example this would be 180 / 10 = 18 and 90 / 5 = 18, which is the behaviour I would expect. I hope this makes sense...

[mcara]

the size parameter is given in (lat, lon) order in the Cutout2D class

Could you please provide a reference?

[mcara]

Also, documentation says:

If size is in angular units, the cutout size is converted to pixels using the pixel scales along each axis of the image at the CRPIX location. Projection and other non-linear distortions are not taken into account.

[adonath]

The docstring says as well:

then a square cutout of size will be created. If size has two elements, they should be in (ny, nx) order.

So it is the y axis (lat) first, which is also the behaviour I see in the first (correct) example I provided in my initial post.

[mcara]

Based on my reading of the docs, I think there is an issue either with the documentation or with the code. Indeed, when size is provided in pixels, it is indeed in the (ny, nx) order. It is not clear in what order are axes when specified using "world" units: in the same order as CTYPE, or the same order as when it is specified in pixels (ny, nx). Considering CTYPE can "flip" world axes, this lack of documentation (for the very least) on how to specify size using angular quantities is concerning. If developers intended ang_x = nx * cdelt1 and ang_y = ny * cdelt2, then indeed I would agree with @adonath that this is a bug. If they meant that size should be specified as (ang_x, ang_y) then there is no bug. And if they meant that it should be specified as (lat, lon) (or even (lon, lat)) - that is, the specification of the axis order uses "world coordinates" - then the code is buggy regardless, just because the code seems to ignore that world axes are controlled by CTYPE.

I now believe this is a bug, but I am not yet sure what is the fix: it would be helpful to know original developer's intention.

[pllim]

Hard to tell where this crept in but a quick blame pointed to #4187 . Overall, nddata/utils.py is mostly touched by @larrybradley and @astrofrog .

[mcara]

And since Cutout2D seem to be specifically designed to deal with image extraction boxes in the image space, allowing shape to be specified in angular units seems unnecessary and confusing. I would remove support for this.

For example, specifying a size of (1 deg, 2 deg) - what does it mean? Especially when CD/PC matrix are not diagonal? is it a rhomb? Are these dimensions along world lon/lat lines? Is it "natural" to think 1 deg is size along Y axis and 2 deg is the size along X axis? So, wouldn't it be better to let the user to figure out the correct extraction box in image coordinates?

[mcara]

What's interesting, is that the code is written in such a way that it would allow size specification of the form (1 * u.pix, 2 * u.deg)!

I would remove support for angles altogether, but if they must be used, I would re-design the code to:

1. check that all units are of the same kind;
2. get rid of the for axis, side in enumerate(size) loop - it could be replaced with vectorized code, I think.

[adonath]

@mcara I think specifying the width / height of the cutout box as an angle is a nice additional convenience from the user perspective. For projection regions with few distortions this just works fine and I think Cutout2D documents correctly, how the conversion is handled internally (the green note https://docs.astropy.org/en/stable/api/astropy.nddata.utils.Cutout2D.html#astropy.nddata.utils.Cutout2D), so I would argue the behaviour is probably fine as is.

I certainly agree on the mixed units, I think this is confusing and super rare as a use case and not worth supporting. I could also try to replace the for loop in #11363 with a corresponding Numpy expression.

[larrybradley]

It has been a while, but my memory of this is essentially what is here in the docs:

"If size is in angular units, the cutout size is converted to pixels using the pixel scales along each axis of the image at the CRPIX location. Projection and other non-linear distortions are not taken into account."

So, essentially nx = ang_x / x_pixscale, etc. (I think @mcara said this above).

[adonath]

Thanks @larrybradley! However I think the statement is ambiguous:

"If size is in angular units, the cutout size is converted to pixels using the pixel scales along each axis of the image at the CRPIX location. Projection and other non-linear distortions are not taken into account."

Does "axis of the image" correspond to the axis of the data or the axis of the WCS? I think internally Cutout2D still relies on the numpy convention to map the first axis of the WCS to the 2nd of the numpy array? Or is it not the case?

[larrybradley]

"axis of the image" refers to the axis of data, e.g. (x, y) not WCS axis. In general (lat, lon) (in whichever order) aren't going to be aligned along the (x, y) axes, so I don't know another way to do this other than to use proj_plane_pixel_scales (the data cutouts are always rectangles in (x, y) space). We should update the documentation if that is confusing. Perhaps even explicitly by stating the cutout pixel size will be:

• nx_cutout = x_size (ang) / x_pixscale
• ny_cutout = y_size (ang) / y_pixscale

where the x and y pixel scales are computed using proj_plane_pixel_scales

Would that help?