ECEF to LLA

[karatemir]

I was trying to see whether I can jump between ECEF and LLA and although I was able to; I feel like my solution is not very clean:
I use the following script to do the conversion LLA --> ECEF --> LLA:

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

ts = load.timescale()
t0 = ts.now()

p = wgs84.latlon(50, 10)
d = p.itrs_xyz.m
print(f'ECEF:{d}')

ecef_to_wgs = wgs84.geographic_position_of( Geocentric(Distance(m=ITRF_to_GCRS(t0, d)).au, t=t0) )
print(f'LLA: {ecef_to_wgs}')

which prints:

ECEF:[4045456.40534089  713323.11354148 4862789.03770643]
LLA: WGS84 latitude +50.0000 N longitude 10.0000 E elevation 3.4e-08 m

Is there a more straightforward way to do this? I feel like wgs84.geographic_position_of( Geocentric(Distance(m=ITRF_to_GCRS(t0, d)).au, t=t0) ) is not the proper way to do this.

Thanks.

[brandon-rhodes]

Alas, as Skyfield is basically a library for building and combining vectors in the fixed ICRS reference frame, it does not provide quick conversions between the two different forms of non-ICRS coordinate you are dealing with: an $(x,y,z)$ ITRS vector and the equivalent latitude-longitude-altitude designation for that same ITRS position.

The logic for computing latitude-longitude-altitude, alas, is currently hidden in a method Geoid._compute_latitude() which is currently only for the use of the public Geoid methods—all of which want an ICRS Position as their input, not a vector that’s already rotated into the ITRS.

This limitation is reflected in an asymmetry in Skyfield's public interface: if you create an ITRS position by supplying a lat-lon and getting back a GeographicPosition, then it will remember its latitude and longitude and you can ask it for those angles again later. But if you create a simple ITRSPosition by supplying a plain $(x,y,z)$ vector, then Skyfield can't be asked what its latitude and latitude is; not out of deliberate spite, but because no one ever asked for that functionality.

So my guess is that we should supplement ITRSPosition so that it can compute its corresponding GeographicPosition.

It almost half feels like they could have been the same class, with two different constructors, one taking $(x,y,z)$ and one taking lat-lon, and then which provided some way to auto-construct the lat-lon on-the-fly if the user tried asking for it. But I've tried to keep automatic conversions in Skyfield to a minimum. The cost here will be the creation of two objects, both an ITRSPosition and a GeographicPosition, for what is really a single idea, a single vector, that's being represented.

Or maybe ITRSPosition should remain very simple, and GeographicPosition should gain a new constructor that takes an $(x,y,z)$?

I'll weigh the options and maybe try writing some code later this month or next, and see which feels the most natural.

[karatemir]

Thanks for detailed response! I also think they should be one class with multiple constructors.

[hunterdsp]

I was looking for the same feature and came up with this for now:

from skyfield.api import wgs84
from skyfield.units import Distance
from skyfield.toposlib import ITRSPosition

class MyITRSPosition(ITRSPosition):
    def init(self, xyz):
        super().__init__(xyz)
    def frame_xyz(self, frame):
        return self.itrs_xyz

position = MyITRSPosition(Distance(m=[-2.66e5, 1.18e6, 7.09e6]))
lat, lon = wgs84.latlon_of(position)
height = wgs84.height_of(position)

print(' {:.4f}° latitude'.format(lat.degrees))
print(' {:.4f}° longitude'.format(lon.degrees))
print(' {:.0f} km above sea level'.format(height.km))