Coordinate conversion problem of a FITS file

Question

I have loaded and plotted a FITS file in python. With the help of a previous post, I have managed to get the conversion of the axis from pixels to celestial coordinates. But I can't manage to get them in milliarcseconds (mas) correctly. The code is the following

import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.wcs import WCS
from astropy.io import fits
from astropy.utils.data import get_pkg_data_filename

filename = get_pkg_data_filename('hallo.fits')


hdu = fits.open(filename)[0]
wcs = WCS(hdu.header).celestial
wcs.wcs.crval = [0,0]

plt.subplot(projection=wcs) 
plt.imshow(hdu.data[0][0], origin='lower') 
plt.xlim(200,800)
plt.ylim(200,800)
plt.xlabel('Relative R.A ()')
plt.ylabel('Relative Dec ()')
plt.colorbar()

The output looks like

The y-label is cut for some reason, I do not know.

As it was shown in another post, one could use

wcs.wcs.ctype = [ 'XOFFSET' , 'YOFFSET' ]

to switch it to milliarcsecond, and I get

but the scale is incorrect!. For instance, 0deg00min00.02sec should be 20 mas and not 0.000002! Did I miss something here?

Answer 1

Looks like a spectral index map. Nice! I think the issue might be that FITS implicitly uses degrees for values like CDELT. And they should be converted to mas explicitly for the plot. The most straightforward way is to multiply CDELT values by 3.6e6 to convert from degrees to mas. However, there is a more general approach which could be useful if you want to convert to different units at some point:

import astropy.units as u
w.wcs.cdelt = (w.wcs.cdelt * u.deg).to(u.mas)

So it basically says first that the units of CDELT are degrees and then converts them to mas.

The whole workflow is like this:

def make_transform(f):
    '''use already read-in FITS file object f to build pixel-to-mas transformation'''
    print("Making a transformation out of a FITS header")

    w = WCS(f[0].header)
    w = w.celestial
    
    w.wcs.crval = [0, 0]
    w.wcs.ctype = [ 'XOFFSET' , 'YOFFSET' ]
    w.wcs.cunit = ['mas' , 'mas']
    w.wcs.cdelt = (w.wcs.cdelt * u.deg).to(u.mas)
    print(w.world_axis_units)
    return w

def read_fits(file):
    '''read fits file into object'''
    try:
        res =  fits.open(file)
        return res
    except:
        return None

def start_plot(i,df=None, w=None, xlim = [None, None], ylim=[None, None]):
    '''starts a plot and returns fig,ax .
    xlim, ylim - axes limits in mas
    '''
       
    # make a transformation
    # Using a dataframe
    if df is not None:
        w = make_transform_df(df)         
    # using a header    
    if w is not None:
        pass
    # not making but using one from the arg list
    else:
        w = make_transform(i)

#    print('In start_plot using the following transformation:\n {}'.format(w))


    fig = plt.figure()
    
    if w.naxis == 4:
        ax = plt.subplot(projection = w, slices = ('x', 'y', 0 ,0 ))
    elif w.naxis == 2:
        ax = plt.subplot(projection = w)
        
    
    # convert xlim, ylim to coordinates of BLC and TRC, perform transformation, then return back to xlim, ylim in pixels
    if any(xlim) and any(ylim):
        xlim_pix, ylim_pix = limits_mas2pix(xlim, ylim, w)
        ax.set_xlim(xlim_pix)
        ax.set_ylim(ylim_pix)
        
    
    fig.add_axes(ax)  # note that the axes have to be explicitly added to the figure
    return fig, ax 


rm = read_fits(file)
wr = make_transform(rm)
fig, ax = start_plot(RM, w=wr, xlim = xlim, ylim = ylim)

Then just plot to the axes ax with imshow or contours or whatever. Of course, this piece of code could be reduced to meet your particular needs.