The implementation of the World Coordinate System (WCS) in python through astropy to handle fits files is quite useful and neat. However, I often find the documentation on it a little lacking. So here I assemble some useful snippets that I often need.

 

Here is a typical WCS for a 2D image: (1D and 3D is also possible but looks the same)

Number of WCS axes: 2
CTYPE : ‘RA—TAN’ ‘DEC–TAN’  #Projection type
CRVAL : 16.65733 3.54336 #central value
CRPIX : 447.0 452.0 #central pixel
PC1_1 PC1_2 : 0.0 1.0   #rotation matrix describing the orientation
PC2_1 PC2_2 : -1.0 0.0
CDELT : -8e-05 8.00000000000003e-05 #scaling
NAXIS : 901 906 #picture size

It works by matching one position on the image with a sky position and then adding a rotation matrix describing the orientation and scaling. This allows to uniquely transform between any image pixel and the corresponding sky position. Three is different conventions for rotation and scaling. Often The two are combined into the CD matrix. Internally astropy always converts to the PC+CDELT convention and therefore you don’t have to worry about it if you use the built-in methods of astropy to access the transformation.

In python, we can handle the WCS with astropy. Typically I end up needing the following imports:

from astropy.wcs import WCS
from astropy.wcs import Wcsprm
from astropy.io import fits
from astropy.wcs import utils

Note that there is a high-level WCS object that is great when just using the transformation. When you want to change the transformation, for example, calibrate the astrometry then I often find Wcsprm more convenient (see further down).

Creating a WCS object

…open fits file in hdu
hdu.header
wcs = WCS(hdr)

Transforming between pixel and sky coordinates

coord = SkyCoord(ra=10*u.deg, dec=10*u.deg)
target = utils.skycoord_to_pixel(coord, wcs) #can handle list of objects
–
coord_again = utils.pixel_to_skycoord(xp, yp, wcs) #can handle list of objects ([x1,x2],[y1,y2])

 

Plotting the data

plt.subplot(projection=wcs)
plt.imshow(hdu.data,cmap=”Greys”, origin=’lower’,vmin=15, vmax=1000)

coord = SkyCoord(ra=10*u.deg, dec=10*u.deg) #overplotting with an object in the field of view
target = utils.skycoord_to_pixel(coord, wcs)
plt.scatter(target[0],target[1], s=50, facecolors=’none’, edgecolors=’b’)

Notes:

  • for quick exploration, I often see more on an image with logarithmic scaling (from matplotlib.colors import LogNorm; plt.imshow(… norm=LogNorm())
  • the projection seems to only apply to plt.imshow and maybe plt.contour. So for plt.plot or plt.scatter you have to give the pixel coordinates, not ra and dec!

Getting the pixel scale of the image

pixel_scale = utils.proj_plane_pixel_scales(wcs)
—> returns a array with pixel scale in x and y direction, both in deg/pixel