Conversion from ECI to geodetic coordinates

[Amaury Dehecq]

Hello,

I'm not really an astropy user at the moment, but I could be interested in some of its functionalities.
I need to convert artificial satellite coordinates in Earth-Centered Inertial (Earth rotation model J2000, referenced to the WGS84 ellipsoid) coordinates to geodetic, i.e lat/lon/alt coordinates. I found a solution using Matlab, with the function eci2lla, which I'm calling from within Python, but I would really appreciate if I could do it directly in Python. I just saw your project and I thought you could help. Is there a relatively easy way to do it using astropy? If yes, how would you do it? Otherwise, could you point out a Python or C/C++ package to help me do it?
Thank you very much!

Amaury

[Marten van Kerkwijk]

Hi Amaury,

Have a look at EarthCoordinate -- while this may sound strange, it
does seem to do exactly what you want, convert X,Y,Z relative to
Earth's centre to lon/lat/altitude relative to some ellipsoid (WGS84
being the default).

See http://docs.astropy.org/en/latest/api/astropy.coordinates.EarthLocation.html#astropy.coordinates.EarthLocation

All the best,

Marten

p.s. Internally, this uses the IAU standards of astronomy routines.

[Amaury Dehecq]

Hi Marten,

thank you for your answer. This is a first step to what I'm looking for. I successfully converted my coordinates to lat/lon but it looks like the conversion goes from ECEF (Earth fixed) to lat/lon and not ECI to lat/lon as I asked.
For the ECI conversion, the time is also required as an input, and I can't find how I could input this information to the function. Any idea?

Below is what I've done so far:
###
from astropy.coordinates import EarthLocation
from astropy import units as u

# x,y,z are my ECI coordinates

# convert array to astropy quantities
xq=u.quantity.Quantity(x,u.km)
yq=u.quantity.Quantity(y,u.km)
zq=u.quantity.Quantity(z,u.km)

# create an Earthlocation instance
test=EarthLocation(xq,yq,zq)

# run the conversion
lon2, lat2, alt2 = test.to_geodetic()

Best,

Amaury

[Marten van Kerkwijk]

I fear I simply don't quite know what the astronomical equivalent of
your ECI is. Could it be the same as what astronomers would call
geocentric celestial reference system or GCRS? (centred on earth, but
aligned with celestial sphere) We do have an
EarthLocation.get_gcrs_posvel()` method, which gives coordinates in
that system and requires a time.

-- Marten

> --
> You received this message because you are subscribed to the Google Groups
> "astropy-dev" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to astropy-dev...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Amaury Dehecq]

Yes, you're exactly right, Earth-Centered Inertial frame is centered on Earth and aligned with Celestial objects.
I need to do some checks but the method you're talking about seem to do the job. I'll update the solution once it's working.
Thanks!

Amaury

[Amaury Dehecq]

Ok, I wrote a function to do the conversion from ECI to geodetic, shown below, but it doesn't seem to work. I was able to check for a few points only, but the altitude and latitude seem correct, but the longitude is not, so it probably comes from the earth inertial to earth fixed conversion. Moreover, the error in longitude is not constant, which would mean that the rotation is not correct...

def eci2lla(x,y,z,dt):
    """
    Convert Earth-Centered Inertial (ECI) cartesian coordinates to latitude, longitude, and altitude, using astropy
    Inputs :
    x = ECI X-coordinate (m)
    y = ECI Y-coordinate (m)
    z = ECI Z-coordinate (m)
    dt = UTC time (datetime object)
    Outputs :
    lat = geodetic latitude (radians)
    lon = longitude (radians)
    alt = height above WGS84 ellipsoid (m)
    """
    # create astropy quantifies objects
    xq=u.quantity.Quantity(x,u.m)
    yq=u.quantity.Quantity(y,u.m)
    zq=u.quantity.Quantity(z,u.m)

    # initiate the earthlocation object
    obj=EarthLocation.from_geocentric(xq,yq,zq)

    # convert datetime objects to astropy time objects
    tt=[Time(t,format='datetime') for t in dt]

    # return the position and velocity vector in earth fixed coordinates
    p, v = obj.get_gcrs_posvel(tt)
   
    # now initiate a new object with the earth-fixed coordinates
    obj2 = EarthLocation(p[0],p[1],p[2])

    # conversion to geodetic
    lon, lat, alt = obj2.geodetic

    return lon, lat, alt

Something that does not seem correct to you?
Thanks,

Amaury

[Marten van Kerkwijk]

Hi Amaury, [Other, better coordinate experts: please check this makes sense!]

I think the conversion may simply be the wrong way around: it may be
that you should start with a GCRS, transform that to ICRS, input that
in EarthLocation, and finally convert to geodetic.

I.e., something like:

```
In [1]: from astropy.coordinates import GCRS, ITRS, EarthLocation,
CartesianRepresentation
import
In [2]: from astropy.time import Time

In [3]: import astropy.units as u

In [4]: t = Time('2015-01-01')

In [5]: gcrs = GCRS(CartesianRepresentation(x=1e8*u.m, y=0.*u.m,
z=0.*u.m), obstime=t)

In [6]: itrs = gcrs.transform_to(ITRS(obstime=t))

In [7]: el = EarthLocation.from_geocentric(itrs.x, itrs.y, itrs.z)

In [8]: el.to_geodetic()
Out[8]:
(<Longitude -100.1356126819113 deg>,
<Latitude 0.0839866338805836 deg>,
<Quantity 93621863.04585256 m>)
```

I should stress, though, that I'm really not quite sure whether I'm
identifying frames correctly here; please do check whether ECI is
indeed the same as GCRS, or whether it has another definition
(oriented like FK5??).

All the best,

Marten

[Amaury Dehecq]

Hi Marten,

This seem to do the job, thanks! Below is the function I wrote if this can be useful to anyone:

from astropy.coordinates import GCRS, ITRS, EarthLocation, CartesianRepresentation

from astropy import units as u

from astropy.time import Time

def eci2lla(x,y,z,dt):
    """
    Convert Earth-Centered Inertial (ECI) cartesian coordinates to latitude, longitude, and altitude, using astropy.

    Inputs :
    x = ECI X-coordinate (m)
    y = ECI Y-coordinate (m)
    z = ECI Z-coordinate (m)
    dt = UTC time (datetime object)
    Outputs :

    lon = longitude (radians)

    lat = geodetic latitude (radians)

    alt = height above WGS84 ellipsoid (m)
    """

    # convert datetime object to astropy time object
    tt=Time(dt,format='datetime')

    # Read the coordinates in the Geocentric Celestial Reference System
    gcrs = GCRS(CartesianRepresentation(x=x*u.m, y=y*u.m,z=z*u.m), obstime=tt)

    # Convert it to an Earth-fixed frame
    itrs = gcrs.transform_to(ITRS(obstime=tt))

    el = EarthLocation.from_geocentric(itrs.x, itrs.y, itrs.z)

    # conversion to geodetic
    lon, lat, alt = el.to_geodetic()

    return lon, lat, alt

Best,

Amaury

[Marten van Kerkwijk]

Happy to hear that worked! Might there be a good place in the
documentation that could have put you on the right track (even just
knowing that GCRS is the right frame, I guess...)? If so, would you be
willing to submit a PR with it?

-- Marten