-
-
Notifications
You must be signed in to change notification settings - Fork 1.7k
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Postage stamp cutouts with WCS #3823
Conversation
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 |
@astrofrog I've added the changelog entry and I will add a section in the documentation. |
I added a new documentation page for 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. |
[ nan 4. 5.]] | ||
""" | ||
|
||
if mode not in ['partial', 'trim', 'strict']: |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Isn't that caught in overlap_slices
and extract_array
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, it is caught now in both overlap_slices
and extract_array
. extract_array
orginally didn't have the 'trim'
option (I added it in this PR) when I wrote that check. I removed the now redundant mode
check in Cutout2D
. Thanks for catching that.
Thanks, @patti. I modified the |
@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. |
Postage stamp cutouts with WCS
Thanks @larrybradley and @mwcraig : Cutout2D is a really useful class.
|
@ofors For 1, A generic function to determine if a For 2, I've never used TPV images. |
@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. As 2, I tried Cutout2D with a TPV P?_? polynomial distorted image and got a 256x256pix cutout.
does this mean that Cutout2D doesn't rotate the cutout according the rotation matrix CD[1,2][1,2] and it only shifts? |
@ofors Neither Note that |
@ofors 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 |
I forgot to mention that my example above requires the current development version of |
Thanks so much @larrybradley for your code snippet. |
@ofors Please see the documentation for creating a 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 |
Thanks @larrybradley ! That worked fine. |
@ofors |
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). |
This PR is an implementation of #3594. I still need to add tests.