sunpy / sunpy

SunPy - Python for Solar Physics
http://www.sunpy.org
BSD 2-Clause "Simplified" License
919 stars 591 forks source link

Inconsistent results from `_parse_submap_input` when using coordinates versus pixel values #6706

Open wtbarnes opened 1 year ago

wtbarnes commented 1 year ago

Describe the bug

When parsing possible inputs to submap, the function _parse_submap_inputs is used to get the pixel coordinates of all corners of the inputs, regardless of whether the inputs are coordinates or quantities, specified in widths and heights, etc.

I would expect then that, for a given set of bottom left and top right corners, whether specified in pixel or world coordinates, that they should return the same set of pixel coordinates. In particular, for the bottom left and top right corners of a map, I would expect that this function should always return ((0,0), (N-1,0), (N-1,N-1),(0,N-1)), for an image of dimensions (N,N), where the ordering of the result here is (bottom left, bottom_right, top_right, top_left). This is indeed the case for pixel coordinates but is not the case for those same coordinates when specified in terms of their world coordinates.

Interestingly, in the world coordinate input case, some of the coordinates are off by several pixels.

To Reproduce

import sunpy.data.sample
import sunpy.map
import astropy.units as u

m = sunpy.map.Map(sunpy.data.sample.AIA_171_IMAGE)

# Get pixel coordinates of corners when specified in world coordinates
corners_coord = m._parse_submap_input(m.bottom_left_coord, m.top_right_coord, None, None)
print('Corners from coordinates')
print(corners_coord)

# Convert coordinates to pixel space and pass them through the same function
blc_pix = m.world_to_pixel(m.bottom_left_coord)*u.pix
trc_pix = m.world_to_pixel(m.top_right_coord)*u.pix
corners_pix = m._parse_submap_input(blc_pix,trc_pix,None, None)
print('Corners from pixels')
print(corners_pix)

gives the result,

Corners from coordinates
(<Quantity [-7.96376298e-11, -9.09494702e-12] pix>, <Quantity [  -2.47269596, 1020.51525106] pix>, <Quantity [1023., 1023.] pix>, <Quantity [1025.47277787,    2.48474914] pix>)
Corners from pixels
(<Quantity [-7.96376298e-11, -9.09494702e-12] pix>, <Quantity [ 1.02300000e+03, -9.09494702e-12] pix>, <Quantity [1023., 1023.] pix>, <Quantity [-7.96376298e-11,  1.02300000e+03] pix>)

System Details

ayshih commented 1 year ago

That AIA map has a slight rotation, which is why a rectangle defined by corners in world coordinates will be slightly rotated relative to the rectangle defined by those same corners in pixel coordinates.

ayshih commented 1 year ago

Marginally related, but I'll xref #4896

ayshih commented 1 year ago

Here's an illustration where I've exaggerated the rotation for clarity. For the same bottom-left and top-right points (green dots), the rectangle defined when those points are provided in pixel coordinates (magenta rectangle) is different from the rectangle defined when those points are provided in world coordinates (blue solid rectangle). Furthermore, in the latter case, submap() crops in pixel space with the desire to include the region defined by the blue rectangle, so it actually returns the region defined by the blue dashed rectangle.

download (75)

wtbarnes commented 1 year ago

Thanks, that picture makes it very clear what's going on and why my assumption of always returning the pixel coordinates of the corners of the image is wrong.

Calling .rotate on the map in my snippet above and then computing the corners gives,

Corners from coordinates
(<Quantity [-5.45696821e-11,  3.41060513e-12] pix>, <Quantity [4.12734574e-05, 1.02700000e+03] pix>, <Quantity [1027., 1027.] pix>, <Quantity [ 1.02700004e+03, -1.49316293e-09] pix>)
Corners from pixels
(<Quantity [-5.45696821e-11,  3.41060513e-12] pix>, <Quantity [1.02700000e+03, 3.41060513e-12] pix>, <Quantity [1027., 1027.] pix>, <Quantity [-5.45696821e-11,  1.02700000e+03] pix>)

The values are as expected, but the ordering is reversed.

wtbarnes commented 1 year ago

Ok I think I know what is going on is that the ordering of the corners returned is not consistent between the two input types. When inputting world coordinates, the ordering of the corners is: bottom left, top left, top right, bottom right. When inputting pixel coordinates, the ordering of the corners is: bottom left, bottom right, top right, top left

The latter is because of an error in the way the top left and bottom right corners are defined in pixel space: https://github.com/sunpy/sunpy/blob/ae2a93f7ed438b1e5a3e7a49231541c246b7c80b/sunpy/map/mapbase.py#L1976-L1977

However, in terms of actually creating a submap, this does not matter because when working out the bounds, the min and max bounds for the x and y pixel coordinates are the only thing that actually matters. Thus, any differences in the ordering of the corners are not important.

dstansby commented 1 year ago

I'm slightly confused by the example given at the top here, since it uses private API. Would it be possible to re-write using public API to make clearer what the potential user-facing issue is here?

wtbarnes commented 1 year ago

I'm slightly confused by the example given at the top here, since it uses private API. Would it be possible to re-write using public API to make clearer what the potential user-facing issue is here?

Albert cleared up what I thought was the user facing issue. The issue with the ordering of the results from the private method is certainly not a user-facing one, but this inconsistency makes it a difficult API for a developer to program against (as I was trying to do in #6707). Additionally, the only reason it doesn't introduce a fairly major bug is because we take the extrema of all the corners when constructing the cutout.

If this is going to be "fixed", my commit in #6707 should be pulled out into its own PR.