# How to flatten 3d astronomical coordinates into 2d polar form in Python

I first asked this at Stack Overflow but got redirected here

I am making a 2d simulation of the solar system. For that you need a starting point for all planets. Right now I am storing that data in this form (ignore the color, radius, mass and Max_Trail_Length columns as they are cosmetics for the model):

Name Color (RGB) Radius (km) Mass (kg) Distance (AU) Velocity (m/s) Direction (radians) Max_Trail_Length
Sun (255, 255, 0) 6963400 1.98847e+30 0 0 0 100
Mercury (255, 0, 0) 2440 3.3011e+23 0.393 47362 pi / 2 200
Venus (255, 165, 0) 6052 4.8675e+24 0.723 35021.4 pi / 2 600
Earth (0, 0, 255) 6371 5.9722e+24 1 29784.8 pi / 2 800
Mars (255, 0, 0) 3390 6.4171e+23 1.52 24130.8 pi / 2 2000
Jupiter (255, 222, 173) 69911 1.8982e+27 5.203 13070 pi / 2 10000
Saturn (210, 180, 140) 58232 5.6834e+26 9.737 9690 pi / 2 20000
Uranus (0, 0, 128) 25362 8.681e+25 19.61 6810 pi / 2 50000
Neptune (0, 0, 255) 24622 1.02413e+26 29.897 5430 pi / 2 200000

As you can see, the direction is the same, so all planets are aligned, and the distance and velocity data is just the average.

For the record, I do not plan to make this into a 3d model, even though it would mitigate this problem entirely.

I want to collect real data from a date (January 1. 1750) using astropy, however when i collect the position and velocity using pos, vel = coordinates.get_body_barycentric_posvel(name, time) (coordinates from astropy), I get 3d data, as you can expect. But I want to convert it into data similar as I at the moment have.

I can't seem to do it right, as flattening the z axis doesn't work (I though it might be ignored at first as the solar system is relatively flat).

Here is the code I have right now for collecting this data.

import csv
import os

import numpy as np
from astropy import coordinates as coord
from astropy import units as u
from astropy.time import Time

# FIXME: Temporary wrong solution
def transform_data(pos, vel):
pos = coord.CartesianRepresentation(pos.x, pos.y, 0 * u.au)
vel = coord.CartesianDifferential(vel.x, vel.y, 0 * u.au / u.s)
vel = vel.norm().to(u.m / u.s).value
dist = pos.norm().to(u.AU).value

return dist, vel

def get_path():
parent_dir = os.path.dirname(os.path.dirname(os.path.abspath(__file__)))
data_dir = os.path.join(parent_dir, 'data')
os.makedirs(data_dir, exist_ok=True)
file_path = os.path.join(data_dir, 'solsystem_data_1750.csv')
return file_path

def process_data(time, writer):
for name in ['Sun', 'Mercury', 'Venus', 'Earth', 'Mars', 'Jupiter', 'Saturn', 'Uranus', 'Neptune']:
pos, vel = coord.get_body_barycentric_posvel(name, time)

dist, vel = transform_data(pos, vel)

angle = np.arctan2(pos.y, pos.x).value

writer.writerow([name, f'{dist:.5f}', f'{vel:.5f}', f'{angle:.5f}'])

def main():
coord.solar_system_ephemeris.set('jpl')

time = Time('1970-01-01')  # TODO: Change to 1750-01-01, but 1750 is a dubious year

file_path = get_path()

with open(file_path, 'w', newline='') as file:
writer = csv.writer(file)
writer.writerow(['Name', 'Distance (AU)', 'Velocity (m/s)', 'Direction (radians)'])

process_data(time, writer)

if __name__ == '__main__':
main()


I want to know if I can either modify my transform_data function, or collect the data in some sort of other way.

I couldn't find any 2d data from astropy, and I haven't found a way to convert the coordinates in a good way yet. I researched some alternatives to astropy, like just collecting the data from NASA's Horizons System, however I there run into the same problem of flattening the data.

• It's not clear what the end goal is. Since all the planets are lined up in your original data, it is not 2d, it's 1d. Some search terms that might help is "orthographic projection" to get from 3d to 2d. And "convert cartesian coordinates to polar coordinates" to get from x,y to radius and angle. Dec 29, 2023 at 17:07
• Well i don't want the planets to be lined up like it is now. I get that it could be interpreted as 2d, but it represents 2d positions. I'll look into orthographic projection. Dec 29, 2023 at 17:29
• Have you considered just using two of the position/velocity vector elements, and ignoring the "up" dimension entirely? Dec 29, 2023 at 23:31
• Yes. I thought that you could ignore it, but I found that it had too much of an impact. Therefore I couldn't simply just set the z-axis to 0 as I did in the code I provided. Dec 30, 2023 at 0:46

The code in question outputs:

Name,Distance (AU),Velocity (m/s),Direction (radians)
Earth,0.90583,30207.28054,1.76660


Direction is reasonable for January 1 but Distance is too short. get_body_barycentric_posvel returns ICRS coordinates, which are equatorial. If I compute pos.z / pos.y for each planet, all are between 0.30 and 0.45, and Earth's is 0.4336 = tan 23.44°.

It sounds like you want ecliptic coordinates instead. BarycentricMeanEcliptic may help. Its default representation is spherical. If you just want the planets' speeds, you can leave their velocities in equatorial form, but include the z component. I suggest something like:

pos_eq, vel_eq = coord.get_body_barycentric_posvel(name, time)

pos_ecl = coord.BarycentricMeanEcliptic(pos_eq)
dist = pos_ecl.distance.to(u.AU).value