crop_by_coords is broken when data is rotated about ndcube HOT 22 CLOSED

ping @tiagopereira

from ndcube.

Issue 100!!!! 🍾

I am puzzled by this too. Our original idea was to make the API of crop_by_coords like GenericMap.submap. I agree this is something that should be fixed.

from ndcube.

So we don't break the api of crop_by_coords in 1.1 I propose we do the following to maintain parity with submap as much as possible:

I propose the method have the following signature:

def crop_by_coords(bottom_left, interval_widths=None, top_right=None):

Where only one of interval_widths and top_right can be specified. (or neither if the shape of the first argument is (2, N).)

We should then add a pending deprecation warning for interval_widths for 2.0.

Also, we should make sure we use the same internal logic as submap where it calculates the pixel positions of the four corners of the physical coordinates and then takes the maximum, so as @tiagopereira says it takes the maximal crop not the minimal one. This second fix should be backported to 1.0.2.

from ndcube.

[Conversation taken from IRISpy chatroom]

Danny: Hi @Cadair . @BaptistePellorceAstro has started working on #100. We should discuss/lay out what the API is to him/in the issue before he gets too far. I was a little confused about the current sunpy GenericMap.submap API only having bottom left and top right kwargs. It seems to me that this does not result in a unique solution to the desired region of interest.
There is a valid solution aligned with the pixel grid and at an angle to the pixel grid. Surely, we need at least one more corner to make this unique.

@Cadair: how so? why? Are you are specifying a rectangle in physical coordinates? Then calculating a pixel box that encompasses the whole of that rectangle

Danny: You mean both possible rectangles? If the user supplies 4 corners, then the crop_by_coords can determine the minimum pixel-grid-aligned rectangle that includes them all. Then a more advanced version could even mask out regions that in that minimum pixel-grid-aligned rectangle not within the user-provided corners.

from ndcube.

As far as I can tell you only need two corners because you can use those two corners to calculate the position of the other two. i.e. https://github.com/sunpy/sunpy/blob/master/sunpy/map/mapbase.py#L1270

from ndcube.

See this image for two valid solutions to the lower left, upper right API.

Should crop_by_coords return the green region? This should definitely include the region the user wants but likely quite a bit more as well. This will probably seem unintuitive. In order to mask out the excess regions, at least one if not two more corners will need to be supplied.

from ndcube.

Yes the returning the green region is what sunpy Map currently does.

The black box in your diagram is almost irrelevant, if the user want's that, that is what slicing is for.

from ndcube.

OK, that's good to have clear up! :)

That's almost true, except in the case where the user has a region of interest in real world coordinates and understandably doesn't want to work out what it is in pixels before cropping.

There is still an advantage of supplying four corners though. And that is that excess region could be masked as part of crop_by_coords. So the user really on does get what they want.

from ndcube.

of course if the pixel and world grids are aligned the black box is a degenerate case of the green one 😉 I think cropping the black box when you only have the corners of a rotated system in world coords is a pretty edge case. Having to convert to pixel yourself is not unreasonable in this case imo.

it would be non-trivial to calculate the mask, I think you would basically have to compute the world coordinates of all the pixels (get's very expensive at higher dimensions).

from ndcube.

Also, you could calculate the mask with only the two corners, as again you can calculate all for corners of the world box given two of them.

from ndcube.

Ah OK, I'm following how this works now! So you always get a rectangle aligned with the real world coordinate system then find the minimum rectangle aligned with the pixel grid that contains the real world aligned rectangle.

Regarding the mask, we could have a kwarg to enable the user have the excess regions masked out. But because it's expensive computationally it be defaulted to False.

from ndcube.

Can't believe how slow I was on the uptake here ;)

@BaptistePellorceAstro does this make sense to you?

from ndcube.

Yes, I read all of that. For now, I am focusing to solve the last problems from SJI but I also look at crop_by_coords. We have had some discussion with Tiago about that because it's not really clear. This discussion follow what Tiago told me and I am trying to understand the code. I think, as Tiago suggested me, I will recreate it on my SJI class for now and I will try to solve the bug by my own work using Sunpy map as reference.

from ndcube.

Hi @BaptistePellorceAstro. I'm a little confused by your answer. Do you mean that you will create a crop_by_coords on your new SJICube?

from ndcube.

Yes, I don't understand all the code and how I have to modify it. I am trying to do it on my SJICube now but only to test it, I think I will not pull it. That can help me to understand clearly the problem and how to solve it.

from ndcube.

Hi, here the code I wrote to solve the bug in 3D (time, y axis, x axis). Do you have an idea to improve that in general cases @DanRyanIrish ?

def crop_by_coords(self, min_coord_values, interval_widths):
n_dim = len(self.dimensions)

So I create a 2D array to stock by lines the coordinates of each point and by columns the time, y axis, x axis.
coord = [ ]

I input the first point, it will be a reference for creating all others points.
coord.append(min_coord_values)

I input the "last" point, the one with all its coordinates are enhanced by the interval width.
coord.append([min_coord_values[i] + interval_widths[i] for i in range(n_dim)])

I input one of the two other points that I will use to crop the data on the y axis.
coord.append([min_coord_values[0], min_coord_values[1] + interval_widths[1], min_coord_values[2]])

I input the second point used to crop by y axis.
coord.append([min_coord_values[0], min_coord_values[1], min_coord_values[2] + interval_widths[2]])

I will now create some lines to limit these data on the boundaries condition. Do you want a new PR to see my entire code ?

from ndcube.

Hi @BaptistePellorceAstro. Here are a few points. First, a couple about API.
If you read through the above discussion between @Cadair and myself on this issue, you should get a better idea of what the API needs to be, namely:

def crop_by_coords(lower_corner, interval_widths=None, upper_corner=None):

We want to move the API to be more like the SunPy submap API, namely inputting lower_corner and upper_corner, rather than interval_widths. This was the intention, but somehow I screwed this up before we released ndcube. In order that we don't break the old API overnight for people already using NDCube, we must also keep interval_widths for the time being. However, only one of interval_widths and upper_corner can be set. Otherwise, an error should be raised. If interval_widths is set, upper_corner should be calculated from lower_corner and interval_widths. Then the derived upper_corner should be used throughout the rest of crop_by_coords. We also need a deprecation warning to be raised if interval_widths is set by the user, saying support for that API will be removed in version 2.0.

The second API point is regarding the format of the inputs. As with the current implementation of crop_by_coords, lower_corner and upper_corner should be iterables of astropy.units.Quantity objects, one Quantity for each dimension, each with a unit compatible with that dimension. This way, we can have two inputs that describe N-dimensions.

Next, some points about implementation. Once you have lower_corner and upper_corner, you can determine the other corners in all dimensions. For example, in the 2D case with time and wavelength as our dimensions:

>>> import astropy.units as u
>>> lower_corner = [u.Quantity(1, unit=u.s), u.Quantity(1, unit=u.nm)]
>>> upper_corner = [u.Quantity(10, unit=u.s), u.Quantity(10, unit=u.nm)]
>>> all_corners = [[1*u.s, 1*u.nm], [1*u.s, 10*u.nm], [10*u.s, 1*u.nm], [10*u.s, 10*u.nm]]

In the ND case, of course, the calculation of all_corners must be generalised using a for loop or, better yet, list comprehension.
Next, NDCube.world_to_pixel must be used to get the pixel values of all the corners.
Then we need to build a new lower_corner and upper_corner (lower_pixels and upper_pixels) made up of the minimum and maximum pixel in each axis.
Next, these minimum (maximum) pixel values must be rounded down (up) to integer values:

lower_pixels = [int(np.floor(l.value)) for l in lower_pixels]
upper_pixels = [int(np.ceil(u.value)) for u in upper_pixels]

Finally, we can build an item with which to slice the NDCube using the regular slicing API:

item = tuple([slice(lower_pixels[i], upper_pixels[i]+1) for i in range(n_dim)]) # The +1 is so that the upper corner will definitely be in the returned cube since in Python slicing, the interval goes from the lower value up to but NOT including the upper value.
return self[item]

This is a mixture of code and pseudo code, of course. Not all of it will work by just by copy-and-paste. There are many details that need to filled in when going from one step to another. But hopefully it gives you a sense of how you can fix crop_by_coords without having to assume a fixed number of dimensions.

I know this might seem a bit opaque or advanced. So let me know if you have questions. And it's always possible I've overlooked something.

from ndcube.

Thanks for that, I will mind about it tomorrow ! I will mainly mind about creating all corners for n dimensions.
Is time alreay in first position in coordinates ?
If yes, do we really need time dependant corners ?
If yes, what is the use for time dependant corners ?

from ndcube.

Time is a dimension and must be treated as though it were any other dimension. After all, although the WCS for SJI includes time, the WCS for the spectrograph doesn't. So in the SJI case, time (in seconds) is included in the two corner inputs. Thus, the implementation for NDCube is independent of whether the WCS includes time or not.

from ndcube.

Ok, I was expecting this answer, just wanted to check.

from ndcube.

😄

from ndcube.

Here the first version of the pull request #113

from ndcube.