Cropping very large fits files using specified boundaries

Question

I have a large fits file (over 30,000 x 30,000) pixels. IRAF cannot handle this size of image. How can one crop a file of this size while retaining correct header information, as IRAF does when using its standard cropping mode?

Answer 1

You can do this sort of cropping with astropy.io.fits, though it's not trivial yet. Since astropy.io.fits uses memory mapping by default, it should be able to handle arbitrarily large files (within some practical limits). If you want non-python solutions, look here for details about postage stamp creation.

from astropy.io import fits
from astropy import wcs
f = fits.open('file.fits')
w = wcs.WCS(f[0].header)
newf = fits.PrimaryHDU()
newf.data = f[0].data[100:-100,100:-100]
newf.header = f[0].header
newf.header.update(w[100:-100,100:-100].to_header())

See also this pull request, which implements a convenience Cutout2D function, though this is not yet available in a released version of astropy. Its usage can be seen in the documentation, modified to include WCS:

from astropy.nddata import Cutout2D
position = (49.7, 100.1)
shape = (40, 50)
cutout = Cutout2D(f[0].data, position, shape, wcs=w)

There are more examples here