Coordinates are not conserved after cropping fits file

Question

Consider the following code (download test.fits):

from astropy.io import fits
from photutils.utils import cutout_footprint

# Read fits file.
hdulist = fits.open('test.fits')
hdu_data = hdulist[0].data
hdulist.close()

# Some center and box to crop
xc, yc, xbox, ybox = 267., 280., 50., 100.
# Crop image.
hdu_crop = cutout_footprint(hdu_data, (xc, yc), (ybox, xbox))[0]
# Add comment to header
prihdr = hdulist[0].header
prihdr['COMMENT'] = "= Cropped fits file")
# Write cropped frame to new fits file.
fits.writeto('crop.fits', hdu_crop, prihdr)

The original (left) and cropped (right) images look like this:

The (ra, dec) equatorial coordinates for the star in the center of the frame are:

Original frame: 12:10:32  +39:24:17
Cropped frame:  12:12:07  +39:06:50

Why are the coordinates different in the cropped frame?

These are the two ways to resolve this, using two different methods.

from astropy.io import fits
from photutils.utils import cutout_footprint
from astropy.wcs import WCS
from astropy.nddata.utils import Cutout2D
import datetime

# Read fits file.
hdulist = fits.open('test.fits')
hdu = hdulist[0].data
# Header
hdr = hdulist[0].header
hdulist.close()

# Some center and box to crop
xc, yc, xbox, ybox = 267., 280., 50., 100.

# First method using cutout_footprint
# Crop image.
hdu_crop = cutout_footprint(hdu, (xc, yc), (ybox, xbox))[0]
# Read original WCS
wcs = WCS(hdr)
# Cropped WCS
wcs_cropped = wcs[yc - ybox // 2:yc + ybox // 2, xc - xbox // 2:xc + xbox // 2]
# Update WCS in header
hdr.update(wcs_cropped.to_header())
# Add comment to header
hdr['COMMENT'] = "= Cropped fits file ({}).".format(datetime.date.today())
# Write cropped frame to new fits file.
fits.writeto('crop.fits', hdu_crop, hdr)

# Second method using Cutout2D
# Crop image
hdu_crop = Cutout2D(hdu, (xc, yc), (xbox, ybox), wcs=WCS(hdr))
# Cropped WCS
wcs_cropped = hdu_crop.wcs
# Update WCS in header
hdr.update(wcs_cropped.to_header())
# Add comment to header
hdr['COMMENT'] = "= Cropped fits file ({}).".format(datetime.date.today())
# Write cropped frame to new fits file.
fits.writeto('crop.fits', hdu_crop.data, hdr)

Gabriel - Ah, OK, thanks for letting me know. When I click the photutils tag now, I get "The photutils tag has no usage guidance, can you help us create it?" and when I click "details" the info is confusing, I don't see what you proposed. I'm not a SO power-user, if you are and know where to report this, please do. It looks like something they could improve in their process or new proposed tags. — Christoph, Nov 26 '16 at 23:06
The tag description I gave is awaiting peer review. I assumed you could see this. Can you post a question regarding this issue in Meta? http://meta.stackoverflow.com/ — Gabriel, Nov 26 '16 at 23:19
@Gabriel Instead of editing the solution into your question please post an answer containing your solution (for example Cutout2D) and remove the solution part from the question again. You can even accept it as answer! — MSeifert, Nov 27 '16 at 16:35
I add them for completeness sake. Christoph could add the code to his own answer, which is the one that best resolved the issue. — Gabriel, Nov 27 '16 at 18:51

score 3 · Accepted Answer · edited Nov 27 '16 at 05:04

3

photutils.utils.cutout_footprint only cuts out the pixels, it doesn't update the WCS which is used to convert between pixel and world coordinates.

Use astropy.nddata.utils.Cutout2D instead.

edited Nov 27 '16 at 05:04

MSeifert

145,886
38
333
352

answered Nov 26 '16 at 23:13

Christoph

2,790
2
18
23

I tried to use `Cutout2D` (`cutout = Cutout2D(hdu_data, (xc, yc), (xbox, ybox))`) but it returns an object that raises a `KeyError: 'object'` when I try to save it: `fits.writeto('crop.fits', cutout, prihdr)`. I tried converting it using `hdu_crop = fits.PrimaryHDU(cutout.data)` but it still results in the same error. – Gabriel Nov 26 '16 at 23:22
Read the `Cutout2D` docstring. It says you have to pass wcs separately, i.e. `from astropy.wcs import WCS; cutout = Cutout2D(data=hdu.data, wcs=WCS(hdu.header)`. Then look at the output `cutout` object. According to the docs, it should have a new `data` and `wcs` attribute for your cutout. You can again turn wcs to a header and make an HDU out of data and header if you need that. Does that help? I'm offline now, but can post working code tomorrow if not. Actually there should be an example with WCS in the `Cutout2D` docs, someone should send a pull request. – Christoph Nov 26 '16 at 23:32
There's two examples using Cutout2D with wcs in the tests: https://github.com/astropy/astropy/blob/master/astropy/nddata/tests/test_utils.py#L441 Should definitely be in the docs, but for now might be helpful for you. – Christoph Nov 26 '16 at 23:36
I've added to my question the two ways to accomplish this, using `cutout_footprint` and `Cutout2D`. Your solution indeed feels simpler, so I'm changing the accepted answer to this one. Thank you both! – Gabriel Nov 27 '16 at 15:01

MSeifert · Answer 2 · 2016-11-26T19:17:31.023

2

The coordinates changed because you sliced the image but you did not alter the WCS informations (especially the reference pixel values).

One way would be to use astropy.WCS:

from astropy.wcs import WCS
wcs = WCS(hdulist[0].header)
wcs_cropped = wcs[280-50 : 280+50 , 267-25 : 267+25]

and then copy this updated wcs to your header:

prihdr.update(wcs_cropped.to_header())

before saving the file.

I'm not sure about what cutout_footprint does so maybe you need to change the slice indices when creating wcs_cropped.

There is an convenience functionality in astropy.nddata named Cutout2D which updates the WCS by default.

edited Nov 26 '16 at 19:17

answered Nov 26 '16 at 18:37

MSeifert

145,886
38
333
352

This brings the coords closer to the original values, but still off: `12:10:37 +39:25:13`. Have you tried this solution? Does it work as expected for you? Also a warning is triggered: `WARNING: FITSFixedWarning: RADECSYS= 'FK5 ' / WORLD COORDINATE FRAME the RADECSYS keyword is deprecated, use RADESYSa. [astropy.wcs.wcs] WARNING: FITSFixedWarning: 'datfix' made the change 'Changed '13/03/95' to '1995-03-13''. [astropy.wcs.wcs]` – Gabriel Nov 26 '16 at 18:53
I think the reason is that photutils uses the (x, y) convention when using `cutout_footprint`. Can you try using `[280-50 : 280+50 , 267-25 : 267+25]` when slicing the WCS? – MSeifert Nov 26 '16 at 19:08
Yes, that was it. The `Cutout2D` class also looks interesting, I'll check it out. Thank you MSeifert! – Gabriel Nov 26 '16 at 19:15
1

@Gabriel I just noticed something weird, your center is `(x, y)` but your box is `(y, x)`. Not sure if that's correct. – MSeifert Nov 26 '16 at 19:19
MSeifert that's correct. That's apparently how `cutout_footprint` works. If you notice the cropped image to the right in my original question, you'll see the crop was applied using a `(y, x)` box. This looked weird to me too at first. Perhaps I'll ask over at the project's Github. – Gabriel Nov 26 '16 at 22:35
1

@MSeifert - I've posted an answer recommending directly to use `Cutout2D`. Doing the WCS changes by hand is more complicated and easy to get wrong, especially switching x/y like discussed here. `Cutout2D` helps, although you still have to read the docs to get x/y order right. – Christoph Nov 26 '16 at 23:16

score 1 · Answer 3 · answered Nov 28 '16 at 19:46

Another answer because it requires an additional package: ccdproc especially the class CCDData which based on astropy.nddata.NDData:

It simplifies reading the file:

from ccdproc import CCDData
ccd = CCDData.read(filename, unit='erg/cm**2/s/AA')

The unit has to be specified because the header unit isn't compatible with the astropy.units module.

The important thing about the CCDData is that you (mostly) don't need to take care of units, wcs, header and masks by yourself, they are saved as attributes and the most basic operations preserve (and update) them accordingly. For example slicing:

xc, yc, xbox, ybox = 267., 280., 50., 100.
ccd_cropped = ccd[int(yc - ybox // 2) : int(yc + ybox // 2), 
                  int(xc - xbox // 2) : int(xc + xbox // 2)]

This sliced ccd_cropped already updated the WCS information, so you can simply continue processing it (like saving it again):

ccd_cropped.write('cropped.fits')

And it should have correctly updated coordinates. The slicing part is actually also possible using astropy.nddata.NDDataRef, only the read and write part are implemented explicitly in ccdproc.CCDData

Coordinates are not conserved after cropping fits file

3 Answers3