<p>The code in question outputs:</p>
<pre><code>Name,Distance (AU),Velocity (m/s),Direction (radians)
Earth,0.90583,30207.28054,1.76660
</code></pre>
<p>Direction is reasonable for January 1 but Distance is too short.
<a href="https://docs.astropy.org/en/stable/api/astropy.coordinates.get_body_barycentric_posvel.html" rel="nofollow noreferrer">get_body_barycentric_posvel</a> returns <a href="https://en.wikipedia.org/wiki/International_Celestial_Reference_System_and_its_realizations" rel="nofollow noreferrer">ICRS</a> coordinates, which are <a href="https://en.wikipedia.org/wiki/Equatorial_coordinate_system" rel="nofollow noreferrer">equatorial</a>.
If I compute <code>pos.z / pos.y</code> for each planet, all are between 0.30 and 0.45, and Earth's is 0.4336 = tan 23.44°.</p>
<p>It sounds like you want <a href="https://en.wikipedia.org/wiki/Ecliptic_coordinate_system" rel="nofollow noreferrer">ecliptic</a> coordinates instead.
<a href="https://docs.astropy.org/en/stable/api/astropy.coordinates.BarycentricMeanEcliptic.html" rel="nofollow noreferrer">BarycentricMeanEcliptic</a> 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:</p>
<pre><code>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
angle = pos_ecl.lon.to(u.radian).value

vel = vel_eq.norm().to(u.m / u.s).value
</code></pre>
<p>If you assume the planets' initial radial velocities are 0, that's true only at perihelion and aphelion, so their orbits won't be exactly right, but at least they won't look too strange.</p>

CC BY-SA 4.0
2023-12-30T00:14:01.670Z
2023-12-30T00:14:01.670Z
Mike G
<p><em>I first asked this at <a href="https://stackoverflow.com/questions/77732740/how-to-flatten-3d-astronomical-coordinates-into-2d-polar-form-in-python">Stack Overflow</a> but got redirected here</em></p>
<p>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 <code>color</code>, <code>radius</code>, <code>mass</code> and <code>Max_Trail_Length</code> columns as they are cosmetics for the model):</p>
<div class="s-table-container">
<table class="s-table">
<thead>
<tr>
<th style="text-align: left;">Name</th>
<th style="text-align: left;">Color (RGB)</th>
<th style="text-align: right;">Radius (km)</th>
<th style="text-align: right;">Mass (kg)</th>
<th style="text-align: right;">Distance (AU)</th>
<th style="text-align: right;">Velocity (m/s)</th>
<th style="text-align: left;">Direction (radians)</th>
<th style="text-align: right;">Max_Trail_Length</th>
</tr>
</thead>
<tbody>
<tr>
<td style="text-align: left;">Sun</td>
<td style="text-align: left;">(255, 255, 0)</td>
<td style="text-align: right;">6963400</td>
<td style="text-align: right;">1.98847e+30</td>
<td style="text-align: right;">0</td>
<td style="text-align: right;">0</td>
<td style="text-align: left;">0</td>
<td style="text-align: right;">100</td>
</tr>
<tr>
<td style="text-align: left;">Mercury</td>
<td style="text-align: left;">(255, 0, 0)</td>
<td style="text-align: right;">2440</td>
<td style="text-align: right;">3.3011e+23</td>
<td style="text-align: right;">0.393</td>
<td style="text-align: right;">47362</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">200</td>
</tr>
<tr>
<td style="text-align: left;">Venus</td>
<td style="text-align: left;">(255, 165, 0)</td>
<td style="text-align: right;">6052</td>
<td style="text-align: right;">4.8675e+24</td>
<td style="text-align: right;">0.723</td>
<td style="text-align: right;">35021.4</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">600</td>
</tr>
<tr>
<td style="text-align: left;">Earth</td>
<td style="text-align: left;">(0, 0, 255)</td>
<td style="text-align: right;">6371</td>
<td style="text-align: right;">5.9722e+24</td>
<td style="text-align: right;">1</td>
<td style="text-align: right;">29784.8</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">800</td>
</tr>
<tr>
<td style="text-align: left;">Mars</td>
<td style="text-align: left;">(255, 0, 0)</td>
<td style="text-align: right;">3390</td>
<td style="text-align: right;">6.4171e+23</td>
<td style="text-align: right;">1.52</td>
<td style="text-align: right;">24130.8</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">2000</td>
</tr>
<tr>
<td style="text-align: left;">Jupiter</td>
<td style="text-align: left;">(255, 222, 173)</td>
<td style="text-align: right;">69911</td>
<td style="text-align: right;">1.8982e+27</td>
<td style="text-align: right;">5.203</td>
<td style="text-align: right;">13070</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">10000</td>
</tr>
<tr>
<td style="text-align: left;">Saturn</td>
<td style="text-align: left;">(210, 180, 140)</td>
<td style="text-align: right;">58232</td>
<td style="text-align: right;">5.6834e+26</td>
<td style="text-align: right;">9.737</td>
<td style="text-align: right;">9690</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">20000</td>
</tr>
<tr>
<td style="text-align: left;">Uranus</td>
<td style="text-align: left;">(0, 0, 128)</td>
<td style="text-align: right;">25362</td>
<td style="text-align: right;">8.681e+25</td>
<td style="text-align: right;">19.61</td>
<td style="text-align: right;">6810</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">50000</td>
</tr>
<tr>
<td style="text-align: left;">Neptune</td>
<td style="text-align: left;">(0, 0, 255)</td>
<td style="text-align: right;">24622</td>
<td style="text-align: right;">1.02413e+26</td>
<td style="text-align: right;">29.897</td>
<td style="text-align: right;">5430</td>
<td style="text-align: left;">pi / 2</td>
<td style="text-align: right;">200000</td>
</tr>
</tbody>
</table>
</div>
<p>As you can see, the direction is the same, so all planets are aligned, and the distance and velocity data is just the average.</p>
<p>For the record, I do not plan to make this into a 3d model, even though it would mitigate this problem entirely.</p>
<p>I want to collect real data from a date (January 1. 1750) using <code>astropy</code>, however when i collect the position and velocity using <code>pos, vel = coordinates.get_body_barycentric_posvel(name, time)</code> (<code>coordinates</code> from <code>astropy</code>), I get 3d data, as you can expect. But I want to convert it into data similar as I at the moment have.</p>
<p>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).</p>
<p>Here is the code I have right now for collecting this data.</p>
<pre><code>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()
</code></pre>
<p>I want to know if I can either modify my <code>transform_data</code> function, or collect the data in some sort of other way.</p>
<p>I couldn't find any 2d data from <code>astropy</code>, and I haven't found a way to convert the coordinates in a good way yet. I researched some alternatives to <code>astropy</code>, like just collecting the data from <a href="https://ssd.jpl.nasa.gov/horizons/app.html#/" rel="nofollow noreferrer">NASA's Horizons System</a>, however I there run into the same problem of flattening the data.</p>

CC BY-SA 4.0
2023-12-29T17:07:03.513Z
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.
Greg Miller
CC BY-SA 4.0
2023-12-29T17:29:21.223Z
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.
Tmpecho
CC BY-SA 4.0
2023-12-29T23:31:57.540Z
Have you considered just using two of the position/velocity vector elements, and ignoring the "up" dimension entirely?
Michael Bonnet
CC BY-SA 4.0
2023-12-30T00:46:16.137Z
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.
Tmpecho
<solar-system><coordinate><modeling><raw-data>
How to flatten 3d astronomical coordinates into 2d polar form in Python
Tmpecho
# How to flatten 3d astronomical coordinates into 2d polar form in Python

(Asked by: username_0 on 2023-12-29)

*I first asked this at [Stack Overflow](https://stackoverflow.com/questions/77732740/how-to-flatten-3d-astronomical-coordinates-into-2d-polar-form-in-python) 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](https://ssd.jpl.nasa.gov/horizons/app.html#/), however I there run into the same problem of flattening the data.

**Tags:** solar-system, coordinate, modeling, raw-data

> **Comments:**
> **username_1**: 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.  
> **username_0**: 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.  
> **username_2**: Have you considered just using two of the position/velocity vector elements, and ignoring the "up" dimension entirely?  
> **username_0**: 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.  

---

## 1 answers

### Answer by username_3 on 2023-12-30. Score: 2 | Accepted answer

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](https://docs.astropy.org/en/stable/api/astropy.coordinates.get_body_barycentric_posvel.html) returns [ICRS](https://en.wikipedia.org/wiki/International_Celestial_Reference_System_and_its_realizations) coordinates, which are [equatorial](https://en.wikipedia.org/wiki/Equatorial_coordinate_system).
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](https://en.wikipedia.org/wiki/Ecliptic_coordinate_system) coordinates instead.
[BarycentricMeanEcliptic](https://docs.astropy.org/en/stable/api/astropy.coordinates.BarycentricMeanEcliptic.html) 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
angle = pos_ecl.lon.to(u.radian).value

vel = vel_eq.norm().to(u.m / u.s).value

```

If you assume the planets' initial radial velocities are 0, that's true only at perihelion and aphelion, so their orbits won't be exactly right, but at least they won't look too strange.


---

