Styjun

I have a FITS image with a header containing information about the coordinates along each axis

NAXIS = 3
NAXIS1 = 259
NAXIS2 = 272
NAXIS3 = 100
CDELT1 = -0.08333333330000001
CDELT2 = 0.08333333330000001
CDELT3 = 0.1
CRPIX1 = 130.5
CRPIX2 = 132.1906015909634
CRPIX3 = 0
CRVAL1 = 255.0
CRVAL2 = 60.0
CRVAL3 = 15.45

Is there an easy way (e.g. and Astropy function) to get the actual (physical) coordinates from the pixel coordinates?

Conversely, is there a built-in function in Astropy to get the coordinates of the closest pixel from actual coordinates?

Edit: I found the WCS.all_pix2world function, but I don't understand how to use it and I didn't find any examples of use.

From the documentation:

There are two accepted forms for the positional arguments:

 - 2 arguments: An N x naxis array of coordinates, and an origin.
 - more than 2 arguments: An array for each axis, followed by an origin. These arrays must be broadcastable to one another.

In the first case, what is N? Wouldn't origin always be [0,0,0]?

• N is the number of points for which you want to do the conversion, i.e. you can compute the pixel or world coordinates for an array of N points.




• origin is just of matter of convention, about the coordinate of the upper left corner. With Python arrays use 0 based indexing, so it is 0. See example below:

In [1]: from astropy.wcs import WCS 

In [2]: from astropy.io import fits 

In [3]: hdr = fits.Header.fromstring(""" 
 ...: NAXIS = 3 
 ...: NAXIS1 = 259 
 ...: NAXIS2 = 272 
 ...: NAXIS3 = 100 
 ...: CDELT1 = -0.08333333330000001 
 ...: CDELT2 = 0.08333333330000001 
 ...: CDELT3 = 0.1 
 ...: CRPIX1 = 130.5 
 ...: CRPIX2 = 132.1906015909634 
 ...: CRPIX3 = 0 
 ...: CRVAL1 = 255.0 
 ...: CRVAL2 = 60.0 
 ...: CRVAL3 = 15.45 
 ...: """, sep='n') 

In [4]: wcs = WCS(hdr) 

In [5]: wcs 
Out[5]: 
WCS Keywords

Number of WCS axes: 3
CTYPE : '' '' '' 
CRVAL : 255.0 60.0 15.45 
CRPIX : 130.5 132.1906015909634 0.0 
PC1_1 PC1_2 PC1_3 : 1.0 0.0 0.0 
PC2_1 PC2_2 PC2_3 : 0.0 1.0 0.0 
PC3_1 PC3_2 PC3_3 : 0.0 0.0 1.0 
CDELT : -0.0833333333 0.0833333333 0.1 
NAXIS : 259 272 100

Now you can compute the world coordinate of the upper left corner:

In [6]: wcs.all_pix2world([[0, 0, 0]], 0) 
Out[6]: array([[265.79166666, 49.06744987, 15.55 ]])

Here you could use an array of Nx3 indices.

So if you have pixel indices and want to convert them to sky coordinates, you need to use origin=0, and the same for the opposite, to convert sky coordinates to pixel indices with wcs.all_world2pix.

Using origin=1 may sometimes be useful if you have pixel indices stored in a catalog using the FITS/Fortran convention.