Support rotation for cutouts

[astrofrog]

Background

As discussed on the tag-up, it would be nice to support specifying a position angle when creating cutouts. I had a look at the existing code:

https://github.com/spacetelescope/mosviz/blob/master/mosviz/utils/cutout_tools.py#L22

and I think it would be easy to integrate.

@pllim - I decided to just write this up for now - would you like to implement what I am suggesting below, or would you like me to do it? (either way is fine, though I don't have much example data to try it on). If you decide to implement it, just let me know if you run into any issues!

Specifying the angle

There are different possibilities for specifying the position angle for the cutouts. I think that if the user doesn't specify anything, it should default to what it does right now, which is to line up the cutout with the pixels.

If the user wants to make sure the cutouts are lined up with some position angle, then we could consider either allowing them to specify a column in the catalog, or as an extra keyword argument to make_cutouts. If the latter, then we should think whether we want to just allow a scalar, or also an array of values (one value per source). I think either way we should use Astropy Quantities to make sure there is no ambiguity in terms of units.

We should make sure there is a difference between angle=None (line up with the pixels) and angle=0 (line up with celestial north).

Making the cutout

If the user doesn't specify a PA, then we should just keep the code exactly as-is for now, because we don't want to do any reprojection.

Let's now assume a PA is specified. The easiest thing to do is then to take the following information:

    * ``'id'`` - ID string; no unit necessary.
    * ``'ra'`` - RA (e.g., in degrees).
    * ``'dec'`` - DEC (e.g., in degrees).
    * ``'cutout_x_size'`` - Cutout width (e.g., in arcsec).
    * ``'cutout_y_size'`` - Cutout height (e.g., in arcsec).
    * ``'spatial_pixel_scale'`` - Pixel scale (e.g., in arcsec/pix).

and create a header from this which does not include any rotation (so we specify CDELT1 and CDELT2 manually), and then if the user specified a non-zero angle, we add CROTA2. For CTYPE I think we could just use RA---TAN and DEC--TAN for simplicity (the projection probably doesn't matter for cutouts)

Once we have the header, then calling reproject is easy:

from reproject import reproject_interp
array, footprint = reproject_interp(imagename, cutout_header)

and array is your new cutout array! :tada:

Advanced options

There are several things one can do to then allow the quality of the cutout to be adjusted. By default it will use bilinear interpolation, but the interpolation order can be adjusted using the order= keyword argument. We could consider adding a similar option to make_cutouts to give this kind of flexibility to the user. We could also allow the user to specify that they want to use flux-conserving reprojection and then use this reproject function:

http://reproject.readthedocs.io/en/stable/api/reproject.reproject_exact.html#reproject.reproject_exact

Future

Future work could include implementing something like this into Astropy core or the reproject package. In addition, we could make it possible for people to specify which north to use (equatorial, galactic, etc.)

[pllim]

Given that we're short on time at the sprint, I would propose some corner cutting -- In the input catalog, when slit_pa column is given and a new function keyword flag that says "yes please rotate", then do the reproject_interp(...). Otherwise, same as before.

So, just to be clear... We are reprojecting just the cutout, right? What I have in mind...

try:
    cutout = cutcls(position, size=(y_pix, x_pix))
except ...:
    blahblah
...

if some_condition_met:
    array = reproject_interp(cutout.data, cutout.wcs.to_header())  # Can I pass in non-HDU here?
    hdr = ???  # How do I construct WCS of the reprojected cutout?
else:
    array = cutout.data
    hdr = cutout.wcs.to_header()

hdu = fits.PrimaryHDU(array)
hdu.header.update(hdr)

blahblah

EDIT: Corrected example.

[larrybradley]

@pllim If you rotate/reproject just the cutout then you will literally cut corners. :smiley: At the very least, you would need to make a cutout that is sqrt(2) times larger than the longest cutout side length to ensure no corners are cut after rotation.

[astrofrog]

@pllim - actually I was thinking it would make more sense to actually just reproject using the original image, which won't incur a loss in performance, because all pixels outside the cutout footprint will be ignored. The reason I suggest that is because otherwise when you rotate the cutout, there will be some empty areas. It's far easier to simply construct the header of the cutout you want including rotation, and reproject from the original image/mosaic. In pseudo-code:

if some_condition_met:
    wcs = WCS(naxis=2)
    wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']
    wcs.wcs.... # basically just set these using the central RA/Dec, pixel scales, rotation optionally
    hdr = wcs.to_header()
    array = reproject_interp(original_image, hdr)
else:
    cutout = cutcls(position, size=(y_pix, x_pix))
    array = cutout.data
    hdr = cutout.wcs.to_header()

hdu = fits.PrimaryHDU(array)
hdu.header.update(hdr)

Does this make sense?

[astrofrog]

@pllim If you rotate/reproject just the cutout then you will literally cut corners. 😃 At the very least, you would need to make a cutout that is sqrt(2) times larger than the longest cutout side length to ensure no corners are cut after rotation.

Indeed, and it's far easier to just instead do what I'm suggesting and let reproject handle it all. Using the original image for reprojection doesn't mean all the pixels in there will actually be reprojected, so it shouldn't be slow.

[astrofrog]

I'm happy to provide real/non-pseudo code if the above is confusing btw!

[larrybradley]

I agree, @astrofrog.

[pllim]

Ahaha... I love puns!

Yes, what you both said make sense. I was just confused. I'll come up with a PR shortly, so @astrofrog can concentrate on the Viz's (MOSViz, SpecViz, AstroViz, etc).

[pllim]

@larrybradley , are you bored? Not too late to join the sprint! :stuck_out_tongue_winking_eye:

[larrybradley]

@pllim I think you've got it covered! 😈

[pllim]

Is this the right way to build that reproject WCS?

wcs = ...  # From master image

if some_condition_met:
    cutout_wcs = WCS(naxis=2)
    cutout_wcs.wcs.ctype = ['RA---TAN', 'DEC--TAN']
    cutout_wcs.wcs.crval = [position.ra.deg, position.dec.deg]
    cutout_wcs.wcs.crpix = [x_pix * 0.5, y_pix * 0.5]
    cutout_wcs.wcs.cdelt = wcs.wcs.cdelt
    cutout_wcs.wcs.pc = wcs.wcs.pc

[astrofrog]

Is x_pix defined as the size in pixels? If so, I think you might want to be setting crpix to (x_pix - 1) * 0.5 - that is, the center of a 2-pixel image is crpix=0.5?

I don't think you should set pc because you will preserve the rotation from the original image. So I think you want to set crota[1] to the rotation angle requested.

Should cdelt be taken from spatial_pixel_scale instead of the original image scale?

[pllim]

I think I got the basic functionality implemented in #15 but I am keeping this open for now in case we want to revisit the extra features @astrofrog mentioned.

Just to clarify... Positive angle rotates counter-clockwise, right? Because that is what it is currently doing.

[astrofrog]

In terms of position angle convention, I'm not 100% sure - I think we may actually want it to go the other way - a position angle of 45 is a line to the upper left, so then if we said that was the new north, it would actually rotate the image 45 clockwise.

[astrofrog]

@pllim - by the way, you may want to try and hard-code order=2 or order=3 when doing the reprojection, since that may provide a better 'default' quality.

[pllim]

@astrofrog , @larrybradley -- What are the plans to have Cutout2D eventually support rotation? It can probably just use reproject as an optional dependency (as we use SciPy in other parts of Astropy)? Should I open an issue on Astropy GitHub about this?

[larrybradley]

IMHO, Cutout2D shouldn't be responsible for rotation. That should be in a separate general-purpose function/class.

[astrofrog]

RotatedCutout2D? :)

[larrybradley]

I was thinking more general (not necessarily cutout related), like RotateImage.

[astrofrog]

@larrybradley - there are benefits to be able to get just a rotated cutout - rotating the whole image before calling Cutout2D is going to be slow. I guess I don't see why Cutout2D couldn't include a position angle option?

[larrybradley]

One could rotate after making the cutout (assuming sizes are adjusted). :)

For rotated cutouts, would the the size be automatically expanded to prevent the corners from being cut? And is the input cutout "size" the size before rotation or after rotation? Also, you would need to add keyword options for interpolation schemes (e.g. order) too. I'm not strongly opposed to it, but my initial gut reaction was based on the philosophy that "a function should do one thing (and do it well)" and a general purpose image-rotation function/class is probably useful in other contexts.

[hcferguson]

We certainly want simple way to rotate images, and for the use case where you are making lots of cutouts it will be more efficient to rotate the whole image once and make the cutouts. But I think it might be pretty common to be asking for cutouts just a few at a time (e.g. from GUI tools) in which case having a function that cuts out the bigger size, rotates it, and cuts out the smaller size might be needed for performance. I think users will expect the cutout to be rectangular and filled after rotation, so one would need to make a larger cutout, rotate and then cut again. I'm agnostic about whether this should be part of Cutout2D or a special RotatedCutout2D.

[astrofrog]

@larrybradley - the size would refer to the one after rotation - if you use reproject then you don't need to worry about the corner cutting, with the new fixes to reproject you can basically just reproject from the original large image with minimal memory usage (you don't have to do it as two steps, with the larger cutout etc).

While I agree there's also a place for a general rotation routine, i think there is also room for a function that takes position, size, rotation, resolution, and makes cutouts. I'm fine if that's separate from Cutout2D but I think it does deserve to be a function/class that's not just a general 'rotate' one.

[pllim]

@astrofrog , I implemented your suggestions in #17. Thanks!

[pllim]

@larrybradley , @astrofrog , @hcferguson -- Perhaps we need a rotating APE for Astropy? :smiley:

[astrofrog]

A more general class could also be called ReprojectedCutout2D to make it clear that this is not just about rotation, but also removing distortions, etc.

[pllim]

I think this is no longer an issue for cutout tools...