Add Earth Center Earth Fixed (ECEF) reference frames

[jtegedor]

It would be useful to be able to define an orbit using a state vector in Earth Center Earth Fixed (ECEF) reference frames, such as ITRF2014.

🐞
This is a new feature to facilitate the usage of the package.

🎯
Precise GPS orbits are normally given in ITRF reference frame. Therefore it would be of great help to be able to use ECEF state vectors to define an orbit, to be able to perform mission analysis for GPS-based positioning.

💡
This is just sample code using astropy to do the reference frame transformation

     # These are state vectors in Earth Center Earth Fixed (ECEF reference frame)
    pos = [SV.X[0], SV.Y[0], SV.Z[0]]
    vel = [SV.VX[0], SV.VY[0], SV.VZ[0]]

    from astropy import coordinates as coord
    cartdiff = coord.CartesianDifferential(*vel, unit='m/s')
    cartrep  = coord.CartesianRepresentation(*pos, unit = u.m, differentials=cartdiff)

    # Transform state vector from ECEF to J2000
    itrs = coord.ITRS(cartrep, obstime = epoch)
    gcrs = itrs.transform_to(coord.GCRS(obstime=epoch))

    #This should be coordinates in J2000
    pos = gcrs.cartesian.xyz
    vel = [gcrs.velocity.d_x.to_value(u.km/u.s), gcrs.velocity.d_y.to_value(u.km/u.s), gcrs.velocity.d_z.to_value(u.km/u.s)]

    # Generate orbit using J2000 state vector
    orbit = Orbit.from_vectors(Earth, pos, vel * u.km / u.s )

📋
I am not familiar with the package yet to be able to propose a solution.

[astrojuanlu]

Hi @jtegedor, thanks for your interest in poliastro! (I edited the issue for formatting)

Adding easy frame conversion to poliastro is the focus of our next release, which will be out in mid-September, before OSCW:

https://2018.oscw.space/event/1/contributions/13/

The way to do it right now in poliastro is exactly as you did. We had to use a similar strategy to convert to Heliocentric coordinates here:

http://docs.poliastro.space/en/latest/examples/Catch%20that%20asteroid!.html

I'll update you in this issue when we do the work. We still have to work on an API, but probably we will link the attractors (Sun, Earth, Moon) to the reference frames (ICRS, GCRS, ...) and create a Orbit.transform method that takes an astropy Reference Frame. We will have to take care of the propagation code too, to make sure we don't use the two body equations in a non-inertial frame.

[jtegedor]

Thanks for feedback, I have limited experience with GitHub & poliastro but hope to be able to contribute when possible. This is a very nice package, thanks for the initiative!

[astrojuanlu]

With the improvements in reference frames treatment (see #439 and #446), this particular issue will be the focus of the net release. We will add support for all the IAU body-fixed frames.

[jorgepiloto]

After a deep reading on how Orbit class inherits from BaseState and the different methods and attributes it contains, I decided to start working on this. I thought it would be useful to have a method called 'to_icrs()' that converts orbit in GCRS and returns it into ITRS.

Following the procedure above I ended up with the following code:

def to_itrs(orbit):
    if isinstance(orbit.frame, coord.GCRS):
        
        # These vectors are in GCRS frame
        r = orbit.r
        v = orbit.v
        
        # We build up a reference frame with them
        cart_diff = coord.CartesianDifferential(*v, unit=u.m/u.s)
        cart_rep = coord.CartesianRepresentation(*r, unit=u.m, differentials=cart_diff)
        
        # Conversion is done
        gcrs = coord.GCRS(cart_rep, obstime=orbit.epoch)
        itrs = gcrs.transform_to(coord.ITRS(obstime=orbit.epoch))
        
        # Orbit is defined with proper frame
        r_itrs = itrs.cartesian.xyz
        v_itrs = [itrs.velocity.d_x.value, itrs.velocity.d_y.value, itrs.velocity.d_z.value] * u.km / u.s
        
        return Orbit.from_vectors(orbit.attractor, r_itrs, v_itrs, frame=coord.ITRS())        
        
    elif isinstance(orbit.frame, coord.ITRS):
        print("Already in ITRS frame")
        
    else:
        print("Unknown frame. Cannot make conversion")

This function seems to work, but the problem comes when checking the frame of returned Orbit:

from poliastro.examples import iss
iss_gcrs = iss
iss_itrs = iss.to_itrs()

print(iss_gcrs)
>>> 6772 x 6790 km x 51.6 deg (GCRS) orbit around Earth (♁) at epoch 2013-03-18 12:00:00.000 (UTC)
print(iss_itrs)
>>> 5800 x 6775 km x 51.6 deg (ITRS) orbit around Earth (♁) at epoch J2000.000 (TT)

Not sure if this is the proper way to transform between GCRS and ITRS.

[astrojuanlu]

Hi @IngenieroDeAviones, thanks for looking into this! The problem is here:

        itrs = gcrs.transform_to(coord.ITRS(obstime=orbit.epoch))
        
        # Orbit is defined with proper frame
        r_itrs = itrs.cartesian.xyz
        v_itrs = [itrs.velocity.d_x.value, itrs.velocity.d_y.value, itrs.velocity.d_z.value] * u.km / u.s
        
        return Orbit.from_vectors(orbit.attractor, r_itrs, v_itrs, frame=coord.ITRS())   

You are transforming to ITRS(obstime=orbit.epoch), but then doing frame=ITRS() at the end. You should be consistent and use the same frame.

On the other hand, what would be the expected result of this?

Orbit.from_body_ephem(Earth).to_itrs()

?

[jorgepiloto]

Oh, it is true! Sorry for the error, this should be the code for the return:

return Orbit.from_vectors(orbit.attractor, r_itrs, v_itrs, epoch=orbit.epoch, frame=itrs)

I added the conversion from ICRS to ITRS in the same function for the proposed case, which can be used as a test, since position for Earth in ITRS should led {0, 0, 0}. However got this error:

# Earth from ICRS
earth_icrs = Orbit.from_body_ephem(Earth)
print(earth_icrs)
>>> 1 x 1 AU x 23.4 deg (ICRS) orbit around Sun (☉) at epoch 2019-01-23 11:26:34.929918 (TDB)

earth_itrs = earth_icrs.to_itrs()
print(earth_itrs)
>>> TypeError: Coordinate frame got unexpected keywords: ['obstime']

The code for this conversion is:

icrs = coord.ICRS(cart_rep, obstime=orbit.epoch) #ERROR LINE
itrs = icrs.transform_to(coord.ITRS(obstime=orbit.epoch))

[astrojuanlu]

That's because ICRS has a fixed definition that does not depend on the time of observation :) Further reading in https://github.com/poliastro/poliastro/wiki/Reference-frames#links

[astrojuanlu]

Hi @jtegedor! With the change we just landed in master, thanks to @hrishikeshgoyal in #554, now you can do this (copying and pasting from your original example):

     # These are state vectors in Earth Center Earth Fixed (ECEF reference frame)
    pos = [SV.X[0], SV.Y[0], SV.Z[0]]
    vel = [SV.VX[0], SV.VY[0], SV.VZ[0]]

    from astropy import coordinates as coord
    cartdiff = coord.CartesianDifferential(*vel, unit='m/s')
    cartrep  = coord.CartesianRepresentation(*pos, unit = u.m, differentials=cartdiff)

    # Transform state vector from ECEF to J2000
    itrs = coord.ITRS(cartrep, obstime = epoch)

    # Generate orbit using J2000 state vector
    orbit = Orbit.from_coords(Earth, itrs)

Hopefully, a bit easier. For your information, the orbit will be internally converted to ECI (GCRS). We will release poliastro 0.12 with this functionality in two weeks, but you can test it using pip install git+https://github.com/poliastro/poliastro.git if you want to give us early feedback. With your permission, I am closing the issue, but feel free to comment. Thanks again!