How do I align SDSS sky images in Python using Astropy?

Question

I have a problem with aligning SDSS sky images. I have downloaded numerous image bands from SDDS. The problem is I don't know how to align them.

They are stored in FITS format. Each image consists of 5 bands (r, i, g, u, z). Unfortunately, each photo points slightly in different directions, so in order to use them I need to align them first. Information where they are directed is stored in FITS header and coordinates are given in WCS - World Coordinate System:

CRPIX1  float   X of reference pixel
CRPIX2  float   Y of reference pixel
CRVAL1  float   RA of reference pixel (deg)
CRVAL2  float   Dec of reference pixel (deg)
CD1_1   float   RA deg per column pixel
CD1_2   float   RA deg per row pixel
CD2_1   float   Dec deg per column pixel
CD2_2   float   Dec deg per row pixel

The goal is to shift 4 bands to point exactly the same point as other one (r for example).

I am using Python's astropy package to solve the problem. Here is my code so far:

from astropy.io import fits
from astropy.nddata import Cutout2D
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from astropy.visualization import make_lupton_rgb
import matplotlib.pyplot as plt
from scipy.ndimage import shift, rotate

def align_spectral_bands(list_of_bands):
    g, i, r, u, z = tuple(sorted_list_of_bands)

    ref_file = fits.open(r)
    ref_header = ref_file[0].header
    ref_wcs = WCS(ref_header)
    ref_data = ref_file[0].data

    rotated_g = align_single(ref_header, g)
    rotated_r = align_single(ref_header, r)
    rotated_u = align_single(ref_header, u)
    rotated_i = align_single(ref_header, i)
    rotated_z = align_single(ref_header, z)

    <rest of function>
    

def align_single(ref_header, band_path):
    ref_wcs = WCS(ref_header)

    curr_band_file = fits.open(band_path)
    curr_header = curr_band_file[0].header
    curr_wcs = WCS(curr_header)
    curr_data = curr_band_file[0].data
    curr_shape = curr_data.shape

    target_ra = ref_header['CRVAL1']
    target_dec = ref_header['CRVAL2']

    target_coord = SkyCoord(target_ra, target_dec, unit='deg')
    cutout = Cutout2D(curr_data, target_coord, curr_data.shape, wcs=curr_wcs)
    cutout_header = cutout.wcs.to_header()
    res = cutout.data


    curr_band_file.close()

    return res

When I do run the code and then plot image using i, r, g bands, then I get result as below.

First image presents data after aligning, second one data without aligning. As you can see, result is better than original data, but still not perfect.

I tried using np.roll() with offsets calculated like cutout_header['CRPIX1'] - ref_header['CRPIX1'], but it only makes it worse.

Answer 1

I think you are looking for reprojection. If you have WCS coordinates the reproject package! should help.