How to obtain 2D Cutout of an image from a SkyCoord position?

Question

I am following the example from Astropy docs for 2D Cutout.

The header of my FITS file:

SIMPLE  =                    T / file does conform to FITS standard             
BITPIX  =                  -32 / number of bits per data pixel                  
NAXIS   =                    3 / number of data axes                            
NAXIS1  =                  512 / length of data axis 1                          
NAXIS2  =                  512 / length of data axis 2                          
NAXIS3  =                    3 / length of data axis 3                          
EXTEND  =                    T / FITS dataset may contain extensions            
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H 
SURVEY  = 'DECaLS  '                                                            
VERSION = 'DR8-south'                                                           
IMAGETYP= 'image   '                                                            
BANDS   = 'grz     '                                                            
BAND0   = 'g       '                                                            
BAND1   = 'r       '                                                            
BAND2   = 'z       '                                                            
CTYPE1  = 'RA---TAN'           / TANgent plane                                  
CTYPE2  = 'DEC--TAN'           / TANgent plane                                  
CRVAL1  =            186.11382 / Reference RA                                   
CRVAL2  =           0.15285422 / Reference Dec                                  
CRPIX1  =                256.5 / Reference x                                    
CRPIX2  =                256.5 / Reference y                                    
CD1_1   = -7.27777777777778E-05 / CD matrix                                     
CD1_2   =                   0. / CD matrix                                      
CD2_1   =                   0. / CD matrix                                      
CD2_2   = 7.27777777777778E-05 / CD matrix                                      
IMAGEW  =                 512. / Image width                                    
IMAGEH  =                 512. / Image height                                   

So far what I have tried :

from astropy.coordinates import SkyCoord
from astropy.wcs import WCS

position = SkyCoord(hdu[0].header['CRVAL1']*u.deg,hdu[0].header['CRVAL2']*u.deg)
size = 200*u.pixel

wcs1 = WCS(hdu[0].header)

cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1 )

I run into error in the last line. Error :

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-142-7cc21a13e941> in <module>
----> 1 cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1 )

/Applications/anaconda3/lib/python3.7/site-packages/astropy/nddata/utils.py in __init__(self, data, position, size, wcs, mode, fill_value, copy)
    714         if wcs is not None:
    715             self.wcs = deepcopy(wcs)
--> 716             self.wcs.wcs.crpix -= self._origin_original_true
    717             self.wcs.array_shape = self.data.shape
    718             if wcs.sip is not None:

ValueError: operands could not be broadcast together with shapes (3,) (2,) (3,) 

My guess is it is because naxis =3 in my file and the documentation assumes naxis = 2. Though I am not sure if that is the actual issue here. Can anybody help fix this error?

Answer 1

This might help to retain the other axes from the original image - e.g. frequency and polarization axes - and the cutout image can be viewed using ds9 without error.

from astropy.io import fits
from astropy.coordinates import SkyCoord
from astropy.wcs import WCS
from astropy.nddata import Cutout2D

# Load the FITS file and create a SkyCoord object for 
# the cutout center
sample_coord = [30.1, 40.2]
im_path = "input.fits"
cut_im_path = "cutout.fits"
c = SkyCoord(sample_coord[0], sample_coord[1], unit="deg")
hdu = fits.open(im_path)[0]

cutout_header = hdu.header.copy()

# Create a 2D WCS specifically for the spatial cutout
wcs_full = WCS(cutout_header)
wcs_2d = wcs_full.celestial

# Extract only the 2D spatial slice of the data
data_2d = hdu.data[0, 0, :, :]
size = (256, 256)

# Create the cutout with the 2D WCS
cutout = Cutout2D(data_2d, position=c, size=size, wcs=wcs_2d)

cutout_wcs_to_header = cutout.wcs.to_header()
# Essential line to retain other axes (e.g frequency and Stokes).
cutout_wcs_to_header["WCSAXES"] = hdu.header['NAXIS']

cutout_header.update(cutout_wcs_to_header)

# Save the cutout data with the updated header
hdu_cutout = fits.PrimaryHDU(data=cutout.data, header=cutout_header)
try:
    hdu_cutout.writeto(cut_im_path, overwrite=False)
except OSError:
    print("Output file already exists.")

Answer 2

Since your WCS is 3d, but you're getting a 2d cutout, you need to drop the 3rd dimension. Try cutout = Cutout2D(hdu[0].data[0], position ,size, wcs = wcs1.celestial ), where .celestial is a convenience tool to drop the third dimension in the WCS.