Search code examples
pythonastropyfits

Preserving the WCS information of a FITS file when rebinned

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

Solution

  • 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.