Postage stamp cutouts with WCS - #3823

This PR is ready for review.

Thanks! Before I review this, I was wondering whether it would make sense to also add a section in the actual documentation? (and minor comment: needs a changelog entry)

👍 on adding actual documentation in docs. As a user that does not really know what this is supposed to do, I would want to see example with screenshots (or plt.imshow() results) in the actual doc.

@astrofrog I've added the changelog entry and I will add a section in the documentation.

I added a new documentation page for nddata.utils, which includes the Cutout documentation.

This is ready for review again.

@patti -- finally took a quick look at #3909 -- you have been busy!

Is the issue here, in a nutshell, that a tuple of tuples generalizes nicely to higher dimensions than two but a flat tuple doesn't? That makes sense to me.

@larrybradley -- any objections to merging this in and then @patti doing a PR to make the bounding box change? I don't have any sense of how this fits into, e.g. photutils, but there aren't any issues with backwards compatibility since this is a new feature.

Thanks, @patti. I modified the bbox_original and bbox_cutout to match your convention in Model.bounding_bound.

@mwcraig Hopefully this is ready to merge now (if all tests pass)!

There are 8 failures against numpy dev, but those seem to be unrelated to this PR and they are occurring in at least a couple other PRs, so merging.

Thanks to @larrybradley for putting this together and to everyone who provided feedback.

Thanks @larrybradley and @mwcraig : Cutout2D is a really useful class.
I'm using FITS images with TPV convention and 3rd degree polynomial distortion (PV?_?).
I have two questions:

1. is there a class/method of Cutout2D which returns True/False depending if input (coord,boxsize) overlap with the input image?
   This has been implemented in astLib.WCS.coordsAreInImage(self, RADeg, decDeg), and in def coords_in_image(xc,yc,header,coordsys='celestial'):.
   However, I would like to keep using Cutout2D in astropy, because of the availability of improved wcslib5.9, which is crucial for TPV images.

2. which is the schedule for Cutout2D to support TPV images with PV?_? distortions?
   I can work on testing side.

@ofors For 1, Cutout2D will raise a NoOverlapError is there is no overlap of the cutout with the input image. The behavior of partial overlaps can be controlled with the mode keyword. For example, using mode='strict' will raise a PartialOverlapError if the cutout only partially overlaps the input data. A partial overlap can also be detected with mode='trim'. In that case the cutout image (e.g. cutout.data) will not have the same shape as the requested input shape.

A generic function to determine if a SkyCoord position is in an image (like astLib.WCS.coordsAreInImage) could be useful, but it would belong in astropy.wcs.utils. PRs welcome!

For 2, I've never used TPV images. Cutout2D simply shifts the CRPIX values of the input WCS to create the output WCS. This keeps the tangent point of the projection the same and also works for SIP distortion corrections because CRPIX is the reference pixel for those distortions. If TPV also uses CRPIX as the reference pixel, then I would guess that it should work, but I really don't know anything about TPV. If not, then perhaps someone with TPV experience could make that work in Cutout2D.

@ofors -- nothing to add beyond @larrybradley's comments except to clarify that @larrybradley did all the hard work :)

@larrybradley You're right. I updated my script using Cutout2D and now it's using try & except with mode='partial' or 'trim'. That way, I can tolerate PartialOverlapError exceptions, but I don't cutout the ones throwing NoOverlapError.
I agree a function like coords_in_image would be nice to have (I encourage others to pull a request :). However, right now I managed to live with the exceptions which Cutout2D throws.

As 2, I tried Cutout2D with a TPV P?_? polynomial distorted image and got a 256x256pix cutout.
The original image is 23.8degx16.0deg FoV, 13.17arcsec/pixel.
The problem with the cutout is that when displaying with ds9 is not showing the correct central target coordinates (ra=311.08,dec=-68.09). Instead it is showing (ra=308.54,dec=-63.43).

@larrybradley: "For 2, I've never used TPV images. Cutout2D simply shifts the CRPIX values of the
input WCS to create the output WCS."

does this mean that Cutout2D doesn't rotate the cutout according the rotation matrix CD[1,2][1,2] and it only shifts?
I see @keflavich cutout.py:L115,116 that he is rotating the cutout with CD[1,2][1,2].
Does anybody think this could solve my coordinate offset?

@ofors Neither Cutout2D or @keflavich's cutout.py perform any rotations (why would they?). L115:116 of cutout.py is simply converting (xw, yw) from WCS units into pixel units using the x and y pixel scales (and actually the calculated pixel scales are not correct in cutout.py if CDELT1,2 are not in the header).

Note that Cutout2D does not yet allow shape to be specified in angular units. It's on my "to do" list.

@ofors Cutout2D seems to work fine with TPV images (because like SIP, the reference pixel is also CRPIX). Here's a simple example starting from the FITS file you sent me:

from astropy.io import fits
from astropy.wcs import WCS
from astropy.nddata.utils import Cutout2D
from astropy.wcs.utils import pixel_to_skycoord

f = fits.open('ML0751914_20150701_035326_calib.fits')
data = f[0].data
w = WCS(f[0].header)
cutout = Cutout2D(data, (2000, 2000), (100, 100), wcs=w)

# find the skycoord of a pixel in the cutout image
x_cutout, y_cutout = (5, 10)
skycoord_cutout = pixel_to_skycoord(x_cutout, y_cutout, cutout.wcs)

# now find the skycoord of the corresponding pixel in the original image
x_data, y_data = cutout.to_original_position((x_cutout, y_cutout))
skycoord_original = pixel_to_skycoord(x_data, y_data, w)

print(skycoord_cutout)
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
    (298.10121912, -71.86024496)>

print(skycoord_original)
<SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
    (298.10121912, -71.86024496)>

As you see the SkyCoord match. Perhaps DS9 is giving the wrong coordinates?

I forgot to mention that my example above requires the current development version of astropy (post #3905) to work with FITS headers containing TPV distortions.

Thanks so much @larrybradley for your code snippet.
Does anybody know how to store the cutout in a FITS file with the original image header and the CRPIX[1,2] values of cutout.wcs?
An example would be most than welcome.

@ofors Please see the documentation for creating a FITS file:
http://astropy.readthedocs.org/en/latest/io/fits/#creating-a-new-fits-file

It is trivial to update a FITS header:

>>> header['CRPIX1'] = cutout.wcs.wcs.crpix[0]
>>> header['CRPIX2'] = cutout.wcs.wcs.crpix[1]

Also see the documentation on how to get a FITS header object from a WCS object: http://astropy.readthedocs.org/en/latest/api/astropy.wcs.WCS.html#astropy.wcs.WCS.to_header

Thanks @larrybradley ! That worked fine.
Please let us know when you plan to implement the polynomial distortion of the cutout (not a copy of the original image).

@ofors Cutout2D already handles polynomial distortions (e.g. SIP, TPV). I don't plan on adding anything there.

Oh! I thought the Warning here saying : "The cutout WCS object does not currently handle cases where the input WCS object contains tabular distortions." was referring the polynomial distortions P?_?.

No, that refers to the table lookup distortions (i.e. WCS Paper IV). I'm not planning to handle them anytime soon (if at all).