Repository navigation
Postage stamps from WCS #3594
Description
I see a good amount of duplicate code on GitHub for making "Postage Stamps" from a larger file. A GitHub search for "astropy postage stamp WCS" provides 15 hits. A search for "postage stamp" within the astropy repository yielded nothing, so I think it's safe to say that this functionality does not exist specifically in astropy. Am I wrong? Making "postage stamps" is a common use case in astronomy.
I think it would make sense to add this functionality directly into astropy, probably in the wcs module. I know it's lame to make a feature request rather than a pull request, but when it comes to wcs there is a lot of nuance that is probably over my head. Like I said, there are many existing examples of this functionality, so those would be a great starting point. I would be happy to collaborate with someone on this feature, if it is deemed worthy of integration into astropy.
Activity
astrofrog commentedon Mar 13, 2015
Just a quick note - maybe this can just be a simple mash-up of
astropy.coordinates,astropy.wcsandWCSAxes.gully commentedon Mar 13, 2015
Just a thought- this might be better suited for APLpy too
embray commentedon Mar 13, 2015
Yeah I have a feeling this kind of specific application doesn't belong in the Astropy core, but I'm willing to be wrong on that.
patti commentedon Mar 25, 2015
I agree. I expected to find a function to do this in Astropy- as it seems many other have- before writing my own. I don't see why it wouldn't belong given that it's just a straightforward application of the existing functionality. The actual plotting can be left up to the user if that's the concern. I'll work on this.
astrofrog commentedon Mar 25, 2015
I wonder if this should maybe live in
nddataand be a utility function that would make use of extract_array and basically take an array, WCS, SkyCoord, and then return arrays returned byextract_arrays?The plotting can indeed be left to the user as suggested by @patti, but the function to go from (data, wcs, coordinates) -> (subarrays, subwcs) should live in astropy core.
@patti - feel free to give it a shot!
cdeil commentedon Mar 25, 2015
What exactly is a postage stamp?
Is it something like (or exactly)
extract_arrayor the proposedcropfunction here:astropy/imageutils#4
Note that if this postage stamp function should use
extract_arrayit might not be a good idea to start this PR now, as long as #3333 is still open. (I haven't followed the discussion, but it seems to be an endless discussion because different people have different opinions what this function should do exactly).I'm all for having nice functions to cut out a series of images for a specified list of catalog positions and width / height and maybe even to plot them on a grid of panels in astropy core ... it just seems a bit hard to get consensus how to implement this and get something merged, so @patti , this is my warning that it might be a little frustrating if this is your first PR for Astropy and then you have to re-write it five times before it's merged.
astrofrog commentedon Mar 25, 2015
Right, there are two ingredients here - the ability to extract cutouts, and making actual postage stamps plots. I think that what astropy should provide is definitely the former, and maybe in the future the latter.
@cdeil is right that it might be worth waiting a few days until
extract_arrayis fixed before attacking this, but then I think it would be handy to have a cutout/crop function in the core. While we can go more complicated as suggested in astropy/imageutils#4 and have bounding boxes etc. I think there is definitely room for a function that just takesSkyCoordand radius since it's a common case.patti commentedon Mar 25, 2015
I haven't used extract_array, but the issue is relevant so it's probably worth holding off. The way I've done this is just to take the RA, Dec position and some angular width as input to get the angular bounds and then use WCS.all_world2pix to get the pixel bounds of the sub-array. Nothing too fancy. Return the subarray, pixel bounds, and perhaps RA,Dec maps.
astrofrog commentedon Mar 25, 2015
@patti - as a side note, there are two functions in
astropy.wcs.utilsthat can help with this, which areskycoord_to_pixelandpixel_to_skycoord:https://github.com/astropy/astropy/blob/master/astropy/wcs/utils.py
These will automatically convert the coordinates to the right frame for the WCS (e.g. Galactic coordinates will be converted to FK5 if WCS is in FK5). This could help generalize the cutout code.
embray commentedon Mar 25, 2015
👍 Sounds good so far.
pllim commentedon Mar 25, 2015
Sounds like a nice addition to
imutilstoo. What do you think, @larrybradley ?larrybradley commentedon Mar 25, 2015
Actually, the
Cutoutclass that I wrote forimutils(https://github.com/spacetelescope/imutils/blob/master/imutils/core.py#L333) already handles WCS, at least in the sense that the central position can be aSkyCoord(which is converted to pixel space withskycoord_to_pixel). The cutout shape is always in pixels though (one could have an angular size, but that becomes trickier if there are geometric distortions). I haven't yet updatedCutoutto return an updated WCS for the cutout data, but will soon. That essentially requires shifting the CRPIX_i to keep the tangent point the same (which is also the reference point for geometric distortion corrections).If there is a desire to put this in the Astropy core (and it seems like there is), I'll submit a PR with my
Cutoutclass.larrybradley commentedon Mar 25, 2015
I should mention that I too want to use
extract_arrayinCutoutwhen #3333 is merged.gully commentedon Mar 25, 2015
👍
eteq commentedon May 29, 2015
@larrybradley - so are you planning to extract
Cutoutinto the core? It sounds from this thread like that's much desired, so I think it makes sense to create a PR to do that.astrofrog commentedon May 31, 2015
+1!
larrybradley commentedon Jun 5, 2015
I've created a PR for
Cutoutin #3823. I still need to add tests though.mcara commentedon Jun 11, 2015
I would like to bring to your consideration the following issues:
astrofrog commentedon Jun 11, 2015
@mcara - I think it makes sense to keep CRPIX be the same world position, because it's not just the distortions that need to be recomputed but the spherical projection will be incorrect so the image would have to be reprojected (except in the rare case of the
-CARprojection).larrybradley commentedon Jun 11, 2015
@mcara I don't think it's a problem to have CRPIX outside of the cutout image. The advantage of shifting CRPIX (and not CRVAL) is that the existing SIP coefficients will work fine and the spherical projection remains the same (as @astrofrog mentioned).
Good point on the tabular distortions. Unfortunately, I don't think there's any good way to propagate them and apply them to the cutout image (thankfully these distortions are usually small). If they are present, then I'll issue a warning and remove them from the cutout wcs object. Otherwise, the tabular distortions would be applied incorrectly (wrong pixels) to the cutout image. Thanks.
cdeil commentedon Aug 4, 2015
@larrybradley's
Cutout2Dclass pull request #3823 was merged. The high-level docs are here.What about this feature request issue? Can it be closed? Or is there something left to do?
mwcraig commentedon Aug 4, 2015
Closing, since as @cdeil points out the functionality is in #3823 (which, in the PR comments, mentions this issue but didn't use the magic word "closes" :))