Best performance for decomposing SkyCoord?

Astropy
1

I have an N x N pixel image and I’d like to compute the numerical values of RA and Dec for each of the N**2 pixels.

import astropy.wcs

n = 200
wcs_dict = {
     'CTYPE1': 'RA---SIN', 'CUNIT1': 'arcmin', 'CDELT1':-1.0,
     'CRPIX1': n/2 + 1, 'CRVAL1': 0, 'NAXIS1': n,
     'CTYPE2': 'DEC--SIN', 'CUNIT2': 'arcmin', 'CDELT2': 1.0,
     'CRPIX2': n/2 + 1, 'CRVAL2': 0, 'NAXIS2': n
}
w = astropy.wcs.WCS(wcs_dict)
d = list(range(n))
x = []
y = []
for i in range(n):
    x.extend(n * [i])
    y.extend(d)
coords = w.pixel_to_world(x,y)
ra = []
dec = []
for c in coords:
     ra.append(c.ra.arcmin)
     dec.append(c.dec.arcmin)

Computing coords takes a couple of tens of milliseconds, which is great. However, assembling the ra and dec lists takes several tens of seconds, and of course gets quickly worse as n increases.

Is there a way to do this that gives significantly higher performance?
Thanks

1 Like
2

You don’t need to iterate through the SkyCoord instance (coords) returned by w.pixel_to_world(x, y). You can access the component arrays directly:

ra = coords.ra.arcmin
dec = coords.dec.arcmin

Note that these are NumPy arrays rather than lists, but if you need lists for whatever reason, you can simply wrap each of the above calls in a list() constructor.

Incidentally, you can also speed up the construction of x and y by using NumPy rather than creating lists:

import numpy as np
x, y = np.indices(w.array_shape)

Note that these are 2D rather than 1D, but if you need 1D, you can simply:

x = x.ravel()
y = y.ravel()
3 Likes
3

Additionally, if you aren’t interested in the actual SkyCoord object and just want arrays / lists of the coordinates you can use the wcs.pixel_to_world_values method.

Making this the fastest way I can think of to do what you want:

x, y = np.indicies(w.array_shape)
ra, dec = w.world_to_pixel_values(x, y)
3 Likes