Preserving the WCS information of a FITS file when rebinned

Question

Aim : Rebin an existing image (FITS file) and write the new entries into a new rebinned image (also a FITS file).

Issue : Rebinned FITS file and the original FITS file seem to have mismatched co-ordinates (figure shown later in the question).

Process : I will briefly describe my process to shed more light. The first step is to read the existing fits file and define numpy arrays

from math import *
import numpy as np
import matplotlib.pyplot as plt 
from astropy.visualization import astropy_mpl_style
from astropy.io import fits 
import matplotlib.pyplot as plt
%matplotlib notebook
import aplpy
from aplpy import FITSFigure



file = 'F0621_HA_POL_0700471_HAWDHWPD_PMP_070-199Augy20.fits'
hawc = fits.open(file)

stokes_i = np.array(hawc[0].data)
stokes_i_rebinned = congrid(stokes_i,newdim,method="neighbour", centre=False, minusone=False)

Here "congrid" is a function I have used for near-neigbhour rebinning that rebins the original array to a new dimension given by "newdim". Now the goal is to write this rebinned array back into the FITS file format as a new file. I have several more such arrays but for brevity, I just include one array as an example. To keep the header information same, I read the header information of that array from the existing FITS file and use that to write the new array into a new FITS file. After writing, the rebinned file can be read just like the original :-

header_0= hawc[0].header
fits.writeto("CasA_HAWC+_rebinned_congrid.fits", rebinned_stokes_i, header_0, overwrite=True)

rebinned_file = 'CasA_HAWC+_rebinned_congrid.fits'
hawc_rebinned= fits.open(rebinned_file)

To check how the rebinned image looks now I plot them

cmap = 'rainbow'

stokes_i = hawc[0]
stokes_i_rebinned = hawc_rebinned[0]

axi = FITSFigure(stokes_i, subplot=(1,2,1))  # generate FITSFigure as subplot to have two axes together
axi.show_colorscale(cmap=cmap)              # show I


axi_rebinned = FITSFigure(stokes_i_rebinned, subplot=(1,2,2),figure=plt.gcf())
axi_rebinned.show_colorscale(cmap=cmap)              # show I rebinned

# FORMATTING
axi.set_title('Stokes I (146 x 146)')
axi_rebinned.set_title('Rebinned Stokes I (50 x 50)')
axi_rebinned.axis_labels.set_yposition('right')
axi_rebinned.tick_labels.set_yposition('right')
axi.tick_labels.set_font(size='small')
axi.axis_labels.set_font(size='small')
axi_rebinned.tick_labels.set_font(size='small')
axi_rebinned.axis_labels.set_font(size='small')

As you see for the original and rebinned image, the X,Y co-ordinates seem mismatched and my best guess was that WCS (world co-ordinate system) for the original FITS file wasn't properly copied for the new FITS file, thus causing any mismatch. So how do I align these co-ordinates ?

Any help will be deeply appreciated ! Thanks

Answer 1

As @astrochun already said, your re-binning function does not adjust the WCS of the re-binned image. In addition to reproject and Montage, astropy.wcs.WCSobject has slice() method. You could try using it to "re-bin" the WCS like this:

from astropy.wcs import WCS
import numpy as np
wcs = WCS(hawc[0].header, hawc)
wcs_rebinned = wcs.slice((np.s_[::2], np.s_[::2]))
wcs_hdr = wcs_rebinned.to_header()
header_0.update(wcs_hdr)  # but watch out for CD->PC conversion

You should also make a "real" copy of hawc[0].header in header_0= hawc[0].header, for example as header_0= hawc[0].header.copy() or else header_0.update(wcs_hdr) will modify hawc[0].header as well.

Answer 2

I'm posting my answer in an astropy slack channel here should this be useful for others.

congrid will not work because it doesn't include information about the WCS. For example, your CD matrix is tied to the original image, not the re-binned set. There are a number of way to re-bin data with proper WCS. You might consider reproject although this often requires another WCS header to re-bin to.

Montage (though not a Python tool but has Python wrappers) is potentially another way.