ICRS to galactocentric conversion using Astropy

Astropy
, , ,
1

Hi all ,
I am using astropy package for converting RA, Dec and distance ( Line of Sight) to Galactocentric coordinates. The problem is the following:

I generated a sphere of randomly distributed 1 million points in RA , Dec and distance plot. It is like a synthetic cluster located at 4000 pc and having radius of 1 pc. Now when converted to Galactocentric frame, I am getting all points flattened to an ellipse. I intuitively expected it to still remain a sphere, as we are just looking from a different location , thus a sphere should, remain a sphere. Is this output expected , or am I doing something wrong ? Following is the code and corresponding output:

from astropy.io.votable import parse, writeto
from astropy.table import Table, Column
import numpy as np
from astropy.io import fits
from astropy.coordinates import SkyCoord, Galactic, CartesianRepresentation
from astropy import units as u

Generating Synthetic Cluster in ICRS

Center coordinates

center_ra = 343.81252 # degrees
center_dec = 62.53608 # degrees
center_distance = 4000 # parsecs

Number of stars

num_stars = 10**6

Radius of the sphere

radius = 1 # this is in pc

distances = (np.random.uniform(0, 1, num_stars))**(1/3) * radius

Generating random spherical coordinates

theta = np.random.uniform(0, 2*np.pi, num_stars)
phi = np.arccos(2 * np.random.uniform(0, 1, num_stars) - 1)

Converting spherical coordinates to Cartesian coordinates

ra = center_ra + np.arctan(distances * np.sin(phi) * np.cos(theta)/center_distance)
dec = center_dec + np.arctan(distances * np.sin(phi) * np.sin(theta)/center_distance)
distance = center_distance + distances * np.cos(phi)

Creating a FITS table with the generated data (RA, Dec and Distance)

hdu = fits.BinTableHDU.from_columns([
fits.Column(name=‘RA’, format=‘D’, array=ra),
fits.Column(name=‘Dec’, format=‘D’, array=dec), This text will be hidden
fits.Column(name=‘Distance’, format=‘D’, array=distance)
])

Save the FITS table to a file

hdu.writeto(‘synthetic_star_cluster.fits’, overwrite=True)

Output in RA, Dec and Distance plot is as expected ( a perfect spherical synthetic cluster, refer the figure at the end.)

Now I am reading this fits file to convert coordinates to galactocentric frame :

import astropy.units as u
import astropy.coordinates as coord

fits_data = fits.getdata(‘synthetic_star_cluster.fits’)
table = Table(fits_data)

Extracting data from the loaded data array

ra = table[‘RA’]
dec = table[‘Dec’]
distance = table[‘Distance’]

c = coord.SkyCoord(ra=ra * u.degree,
dec=dec * u.degree,
distance=(distance) * u.pc,
frame=‘icrs’)
galactocentric_cartesian = c.transform_to(coord.Galactocentric(galcen_distance = 8.3 * u.kpc , z_sun = 27 * u.pc))

Define new columns for galactocentric coordinates

X_column = Column(data=galactocentric_cartesian.x.value, name=‘X’, dtype=‘float64’)
Y_column = Column(data=galactocentric_cartesian.y.value, name=‘Y’, dtype=‘float64’)
Z_column = Column(data=galactocentric_cartesian.z.value, name=‘Z’, dtype=‘float64’)

Append the new column to the table

table.add_column(X_column)
table.add_column(Y_column)
table.add_column(Z_column)

Write the modified Astropy Table to a FITS file

table.write(‘synthetic_XYZ_out’, format=‘fits’, overwrite=True)

The figure below clearly shows that a spherical synthetic cluster in ICRS frame is a flattened ellipse in galactocentric frame , shouldn’t it remain a sphere? . Why is this happening?

image
image989×739 42.2 KB

2

Your equations for the modified RA and declination are not correct for adding two spherical vectors. RA and declination are coupled, and you’re missing the cos(dec) factor that converts between an angular offset in the RA direction and the corresponding offset in RA units.

I suggest you leverage Astropy’s spherical representations and being able to add them:

from astropy.coordinates import SphericalRepresentation

center = SphericalRepresentation(center_ra, center_dec, center_distance)
offset = SphericalRepresentation(theta, phi - np.pi/2, distances)

c = coord.SkyCoord(center + offset, ...)
3

Thanks @ayshih , that was the right solution :slight_smile:
Sphere does transform to a sphere !