Gabriel Gabriel - 12 days ago 13
Python Question

Coordinates are not conserved after cropping fits file

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:

enter image description here

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?

Answer

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.