Perform coordinates projection with astropy

Question

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)

Emilio Miralles · Accepted Answer · 2016-03-07T18:36:46.330

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.

MSeifert · Answer 2 · 2016-02-22T14:54:55.993

I'm not really sure if that is exactly what you want but I give it a try:

You probably need to use the SkyCoord class:

import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord

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

You must use astropy units to set it up:

a = SkyCoord(ra * u.degree, dec * u.degree)
a
<SkyCoord (ICRS): (ra, dec) in deg
    [(134.26176984, -55.707478), (134.58510684, -55.87649941),
     (134.87524059, -56.48554659), (133.13341559, -56.81440026),
     (132.87143675, -56.44936089), (131.14191795, -56.68300147),
     (133.50910855, -56.07594072), (132.05561955, -56.90231565),
     (134.62961065, -55.72207914), (133.64757046, -55.59277667)]>

Now you have a coordinate class (not sure what coordinate system is most appropriate. SkyCoord defaults to ICRS). If you want to calculate it's azimutal projection you can use other representations I think PhysicsSphericalRepresentation is what you want?!

a.represent_as('physicsspherical')
<PhysicsSphericalRepresentation (phi, theta, r) in (deg, deg, )
    [(134.26176984, 145.707478, 1.0), (134.58510684, 145.87649941, 1.0),
     (134.87524059, 146.48554659, 1.0), (133.13341559, 146.81440026, 1.0),
     (132.87143675, 146.44936089, 1.0), (131.14191795, 146.68300147, 1.0),
     (133.50910855, 146.07594072, 1.0), (132.05561955, 146.90231565, 1.0),
     (134.62961065, 145.72207914, 1.0), (133.64757046, 145.59277667, 1.0)]>

Thanks MSeifert but this is not what I want. This gives coordinates on a sphere of unit radius, I need to project equatorial coordinates onto a plane. — Gabriel, Feb 22 '16 at 14:58
Just as question: Do you want something like [Altitude-Azimuth system](http://docs.astropy.org/en/stable/api/astropy.coordinates.AltAz.html) or do you just want to apply the formulas given in [wikipedia](https://en.wikipedia.org/wiki/Azimuthal_equidistant_projection#Mathematical_definition) applied to your coordinates? — MSeifert, Feb 22 '16 at 15:46
I want to apply a projection on my data from the many defined in `astropy`, for example the [Zenithal equidistant projection](http://astropy.readthedocs.org/en/latest/api/astropy.modeling.projections.Sky2Pix_ZenithalEquidistant.html#astropy.modeling.projections.Sky2Pix_ZenithalEquidistant) — Gabriel, Feb 22 '16 at 15:54
See [here](http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html) for the definition of the "Azimuthal Equidistant Projection" (equivalent to the "Zenithal Equidistant Projection") — Gabriel, Feb 22 '16 at 15:56

Perform coordinates projection with astropy

2 Answers2