CODEWITHSUNDEEP
pythonprojectionastropy
Gabriel
Gabriel

Reputation: 42329

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)

Upvotes: 4

Views: 1163

Answers (2)

Emilio Miralles
Emilio Miralles

Reputation: 129

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.

Upvotes: 3

MSeifert
MSeifert

Reputation: 152587

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)]>

Upvotes: 1

Related Questions