Search code examples
pythonprojectionastropy

Perform coordinates projection with astropy

I have a set of points in equatorial coordinates and I need to project them onto a plane, ie: a zenithal or azimuthal projection.

astropy is apparently able to perform this type of projection, among many others.

The problem is that I don't know how to apply those modules to my data, and the documentation is rather scarce.

For this example, assume the equatorial ra, dec coordinates in degrees can be generated via:

ra = np.random.uniform(130., 135., 1000)
dec = np.random.uniform(-55., -57., 1000)
Solution

  • These projections are used with a wcs (world coordinate system) object and the FITS format. The main steps are:

    import numpy as np
from astropy import wcs

#create data
ra = np.random.uniform(130., 135., 1000)
dec = np.random.uniform(-55., -57., 1000)

#create wcs object
w = wcs.WCS(naxis=2)
# set coordinate and projection types for each axis
w.wcs.ctype = ["RA---ARC", "DEC--ARC"]

#do the projection
pixel_coords = w.wcs_world2pix(ra, dec, 1)

    Note that there are several additional configuration parameters for the wcs item that may be important, such as crpixj, the image reference point for axis j, and pvi_m, intermediate world coordinate axis i parameter value.

    For some usage examples see the astropy.wcs documentation and the projection unit tests.