extract_array can give results of different size - #3333

cc @adonath @mhvk

Cleanest may be to always return a MaskedArray, then the user can decide to do what they want to do if any elements are masked. Indeed, if the small array falls completely outside the large array, one could return a fully masked array. If we make this an option, then I'd suggest that for the normal case the code raises an exception if the requested region does not completely fall within the large array.

I think extract_array should always return an array of the requested shape, so
+1 to @mhvk's suggestion for returning a MaskedArray and raising an exception if there is no overlap.

+1

My suggestion was a bit more complicated, allowing one to ask for either masked or regular array, but always returning a masked array and raising an exception if fully outside is probably easier and thus better.

This turns out to be trickier than I thought, but I'm working on it.

OK, I pushed code. I suspect that this will not make it into 1.0 and thus I have not written anything into the CHANGES.rst yet.

travis-ci failed because of this Sphinx issue:
https://travis-ci.org/astropy/astropy/jobs/48295532#L4253

I'm not convinced this change is the right approach. Maybe it would be better to have functions that compute nd slices and then the user can decide what to do (e.g. change in place, or copy with or without padding) ... maybe using other helper functions.

I don't have a concrete suggestion at the moment, but my suggestion would be to think about this and discuss some more, unless most people think it's important to make this change before the 1.0 release.

@adonath You're using extract_array in gammapy/gammapy#234 ... what do you think?

@cdeil The function to compute ndslices is overlap_arrays. If a user wants just a slice, he/she can call that. All that extract_array does is to return large_array[slices] (padded with masked entries in my branch to ensure that the return array has the size requested).

I should point out that in the course of this PR, I found that the current implementation no only returns arrays or different sizes (which I consider a bug, but others may consider that a feature) it also gives clearly wrong results when called with invalid (negative) inputs (examples with current master, are fixed with the modified implementation in this PR):

>>> from astropy.nddata.utils import overlap_slices, extract_array
>>> l = np.arange(5)
>>> extract_array(l, (2,), (-5,))
array([0, 1])     # correct answer would be empty or (per discussion above) an exception
>>> extract_array(l, (2,), (-3,))
array([0, 1, 2, 3])   # incorrect and much longer than requested
>>> # It works again is the requested position is close to valid
>>> extract_array(l, (2,), (-2,))
array([], dtype=int32)

The problem is independent of the number of dimensions, so I chose 1d arrays in these example, but it is a problem for 2 d arrays, too:

>>> extract_array(la, (2,3), (-5,-4))
array([[0, 1, 2, 3]])

Ah, OK ... I didn't know overlap_slices was part of the public API.

Still ... aren't there use cases where one does want to get smaller cutouts at the edged (and not have those padded or masked)? I think that's what we need in gammapy/gammapy#234 .

Would it be bad to add an option to let the user control the behaviour at the large array edges where the small array doesn't completely fit?

When you use the current version of extract_array you get a cutout of unknown size. Let's say you have some X-ray (or gammy ray) event with a peak in pixel (x, y) and not you want to extract a 5*5 pixel island around it.
Let's say you determine the centroid of that and find it's at (0.5, 2.2) of the small array. How do you know what coordinate that is in the big array? Your code might try (x+0.5, y+2.2), but that is wrong at the edges.
Instead the code in gammypy would have to do:

#event on the left side of array:
if x < 5/2:
    x_event = x - (5 - small_array.shape[0]) / 2 + 0.5
else:
    x_event = x + 0.5
# event on left side of array:
if y < 5/2:
    y_event = y - (5 - small_array.shape[1]) / 2 + 2.2
else:
    y_event = y + 2.2

(Note that the example above does not include the correct rounding for even/odd array sizes and that in practice 5, 0.5 and 2.2 would be variable names.) Does that not look a little more confusing than it should be?

I'm not opposed to offering a keyword to switch that behaviour (that's what I suggested at the top of this thread) but I've been convinced that it would be easier to support one and only one way to do it.

@cdeil For gammapy/gammapy#234 it doesn't make a difference because the boundary handling is done outside the extract_array function. In general, I guess, that's what most users will do anyway, because it depends a lot on the application what you want to do with the extracted 'incomplete' arrays at the boundaries.

@hamogu Sorry for not commenting earlier. I'm not in favor of returning a MaskedArray, because it could change the type of the array. When you extract a smaller array out of a larger one you would expect it to have exactly the same type as the larger one.

I would keep the default behaviour as it is. That leaves is to the user, what to do with those 'incomplete' extracted arrays, because it think there are surely cases were you just don't care and discard them, as it was done in the photutils psf photometry prototype . But I would add an optional boundary parameter for convenience. It could even be the same options as in astropy.convolution.convolve i.e. 'none', 'fill', 'wrap' and 'extend' or possibly 'reflect'. But that might be overkill... and padding with NaN is sufficient. If the smaller array is completely outside the boundaries an error is raised.

Regarding the discussion re: maskedarray one can also return an array with NaN values outside the boundary. That could actually be a keyword argument if needed.

@hamogu - on the negative indices: yes, this is indeed inconsistent. In principle, there is no reason not to be able to give a negative offset -- meaning, as usual, counting from the far end of the array -- but the addition of the size should not be allowed to become positive. Still, the effort would likely not be worth the pain, so it may be easiest to just disallow negative indices.

What does a negative position mean?

I interpreted that as the position in a geometical sense (The array as indexes 0,1,2,3, ... ,so "-1" would be on to the left of "0"), but I now realize that this can also be interpreted as the "index of the position" (negative numbers simply mean it's counted from the right of the array.)

Which one was the intention of the original authors of this code?

If they intended the meaning "array indices" then I take everything back that I said and close this PR. In that case, I just suggest the simple addition of "Index of the Position of the small array's center, with respect to the large array." to the docstring.

@hamogu - for what it is worth, I think your interpretation is much more logical than mine: let's treat pos as a position, not an index, i.e., something that might be returned by a conversion from RA, Dec.

@adonath's point about keeping the array type the same is well-taken (especially since MaskedArray deals so poorly with subclasses). At the same time, I do like the idea that one always gets an array with the right shape back. Filling with nan would seem the best compromise.

I'd second @astrofrog's suggestion that the fill value be a keyword parameter (e.g. with default=np.nan). The user could then later create a MaskedArray for the boundary cases if needed.

@hamogu My intention was position measured in pixel coordinates, rather than index, I just didn't care about negative positions. But as you suggested it should be handled in some way, because I think e.g. wcs_world2pix also returns negative pixel coordinates, when the world position is outside the image.

Just to underline my comment above: I don't think extract_array should always return the given shape, because the information that the array was extracted from the boundary region can be useful for how the array data is processed later. E.g. as I've written above, you might just want to discard this cases or apply a renormalisation. My suggestion would be the same as @astrofrog, but additonally having the boundary handling optional with a parameter boundary='fill' and 'none'. When the option is set to'none' the behaviour could be either as it is now (returning 'incomplete' arrays) or maybe better raising some SmallerArrayIsNotFullyContained error. boundary='fill' could be the default...

Just thinking out loud here, but perhaps a better approach would be to make this a class (e.g. DataCutout?) that returns both a cutout array (which would be truncated at the edges) and a bunch of positional attributes (e.g. center, bounding box, slice tuple, etc.) relative to both the large and small images/data?

I see this more as a low-level function that could be used to build a classes like @larrybradley suggested on top of it. If we are talking e.g. about an image, that class could also handle changes to the WCS etc.
On the other hand, in some classes (e.g. photometry) I probably want to use this operation in an inner loop and I really only want one return value (and as other comments point out above, even the padding with nans might already be too much for some use cases), so making this into a class might be more complicated than necessary.

The Travis failure us just a few sphinx warnings which I will fix when I add the changelog in the end.
Before that, I want to make sure that we agree that this is the way it should be implemented.

@larrybradley: For odd arrays we use the middle element as the center already, which, I believe, is exactly what you suggest. (I looked at both the scipy documentation and the scipy source code and it is not clear to me what happens for even sized arrays).

@astrofrog : If you want to know exactly what the coordinates of the input position are in terms of your new array, use mode='partial'. The position you put in will be exactly the center of the extracted array, with two complications:

• if you put in an float position instead of an int, then it might not be at the center of the pixel, but it's in the same pixel (e.g. if you put in 3.2, the center of the new array will have the same coordinates as 3.0 in the original array)
• for even sized arrays it's rounded up, i.e. the "center" of a (2,2) array is the pixel (1,1), thus pixel (1,1) will have the same coordinates as the position you put in.

However, I understand that this is not necessarily sufficient. A user will want an option that returns the coordinates of position in the new array for even and odd shaped arrays and for any mode. Unfortunately, this is somewhat non-trivial in the presence of edges (and precisely this is the reason why I originally implemented the mode that we now call "partial", where it's easy). I think I will have to redesign overlap_slices to calculate these numbers. Sounds simple, but is a little tedious to get it right for even shaped arrays and rounding, so it won't happen today and maybe not even this weekend.

I'm sorry that now it's me to hold this up. I know that by now there are at least two other PRs that wait for this. Please bear with me.

@hamogu Yes, that's right. I was just trying to clarify that the generic filter origin=0 means the "middle", not index=0 (origin was probably not the best name).

BTW, for even-sized arrays, scipy appears to use (shape / 2) as the "center" index (as you have done here I think).

Travis is still running, but this is a heads up that I added a keyword to return the position in the coordinate system of the new array. Please review.

@hamogu -- only a couple minor comments, thanks for adding the option to return position! If you don't have time to make these changes please speak up. They are not critical to merging this.

@astrofrog -- position is now returned, haven't tried it extensively myself, but the test suite is pretty thorough. Unless you have concerns beyond my line comments, this looks good to merge; those line comments are not critical, can be fixed later.

Any more comments on my longest ever running PR?

But it was created only in this year!

This might not get merged until the python in astronomy workshop is done.

Perhaps we could open a PR for a new unit? 1 hans = 1 comment/sloc or 1 hans = 1 day-since-pr-opened/sloc and any PR that reaches 0.5 hans is automatically merged? 😁

@hamogu: @mwcraig and myself just looked through the comments and agree that this is good to merge. However, it can't be automatically merged. Can you rebase and ping us? Thanks for all your effort!

I get a merge conflict in every single of 19 commits. That's why I don't like PR's that are open this long (and if I remember right, this is about the forth time that I rebase).
I will close this PR and start from scratch from a fresh master checkout and just copy and paste the final changes in a text editor. That also cleans up the commit history, because this entire PR will only be a single commit.
(I understand that technically I could do that on the same branch, but I want to keep this branch in my repro in case anything goes wrong. )

@mwcraig @wkerzendorf : Please see #3727 which is this PR, with copy & paste of the changes into a new branch based of the current master.

I get a merge conflict in every single of 19 commits.

Something about that doesn't sound right. In the future I could probably try rebasing for you.