Document how to convert coordinate from ECI to LLA

[Cory-Kramer]

I've been reading over the documentation and I can't quite figure out how to leverage skyfield to convert a coordinate from earth-centered inertial (ECI, given as x, y, z) to LLA (lat, long, altitude). Consider the following example

from datetime import datetime, timedelta
obs_time = datetime(2000, 1, 1) + timedelta(seconds=240)   # observe at 2000-01-01 00:04:00

# Given x, y, z coordinates in ECI [meters]
eci = (2183457.758387, 1924098.065987, 5898755.881731)

# What I would like to do
lat, lon, alt = eci_to_lla(eci, obs_time)

# expected output should be somewhere over eastern europe
# lat = 63.7398
# lon = 40.3853
# alt = 199473.0 

Is there a class in skyfield.positionlib that I can use to accomplish this calculation?

[brandon-rhodes]

I don't have an editor open on this desktop to try it out, but I'm imagining something like:

from skyfield.api import Distance, load, wgs84
from skyfield.positionlib import Geocentric

ts = load.timescale()
t = ts.utc(2020, 11, 6)
d = Distance(m=[2183457.758387, 1924098.065987, 5898755.881731])
p = Geocentric(d.au, t=t)
g = wgs84.subpoint(p)
print(g.latitude.degrees, 'degrees latitude')
print(g.longitude.degrees, 'degrees longitude')
print(g.elevation.m, 'meters above WGS84 mean sea level')

See whether something like that does what you intend. I'm guessing here that GCRS is a good enough ECI frame for you? You don't mention whether you have a particular one in mind or whether any old ECI frame will do. :)

[Cory-Kramer]

Thanks for the code snippet! In that example where does the GCRS class come from? I'm getting a NameError: name 'GCRS' is not defined and am not getting a hit for that class name in the docs.

[brandon-rhodes]

@Cory-Kramer — Now that I'm back at the keyboard, I have made a few edits, and corrected a couple of names; try out the example again and see if it works better this time.

I suppose the documentation needs a new section on coordinate conversion to handle these kinds of cases.

[Cory-Kramer]

Thank you @brandon-rhodes that worked perfectly! Just had to tweak the Distance line to the following

d = [Distance(m=i).au for i in (2183457.758387, 1924098.065987, 5898755.881731)]
p = Geocentric(d, t=t)

Thanks again for the help.

[brandon-rhodes]

Thank you @brandon-rhodes that worked perfectly!

Great! Now I just need to figure out where this will go in the documentation. Might take a week or two.

Just had to tweak the Distance line to the following…

But, if I'm reading correctly, that would leave d as a plain Python list. Wouldn't that raise an AttributeError in the next line where d is asked for its .au attribute?

[Cory-Kramer]

You're correct, I had to update the line below it as well, edited my comment above to reflect that.