Topocentric Position to Geocentric Frame Conversion

Orekit Python Wrapper usage
1

Hi all,

I am new to this topic and Orekit. I am trying perform IOD using the IodGibbs method. From my understanding this method can take the following call: estimate(Frame frame, Position p1, Position p2, Position p3)

I have right ascension and declination measurements of a satellite and the latitude, longitude and altitude of the ground station. The radec measurements are from the ground station and therefore taken from a topocentric frame. As seen in the code below I use the radec measurement to create a vector of the observation and am trying to create a Vector3D position in the ECI frame so that I can input the position for 3 observations into the iodGibbs function. Could I please get some validation on whether this method of doing a frame conversion is correct? I seem to be getting coordinates that are in the same order of magnitude of the actual J2000 coordinates of the satellite at the specific epoch that I am testing… but the numbers are still fairly off as shown below.

utc_time_scale = TimeScalesFactory.getUTC()

ECI = FramesFactory.getEME2000()  # Define ECI reference frame.
ECEF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)  # Define ECEF reference frame.

earth = OneAxisEllipsoid(
    Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
    Constants.WGS84_EARTH_FLATTENING,
    ECEF,
)

def create_sat_position(ra, dec, lat, lon, alt, rng, earth: OneAxisEllipsoid, date: AbsoluteDate):
    range_in_m = rng * 1000
    alt_in_m = alt * 1000

    ra = radians(ra)
    dec = radians(dec)
    lat = radians(lat)
    lon = radians(lon)

    # Create the topocentric frame of the observer
    observer_point = GeodeticPoint(lat, lon, alt_in_m)
    observer_topo_frame = TopocentricFrame(earth, observer_point, "Observer")

    # Convert Ra and Dec to a direction vector in the topocentric frame
    sat_direction_topo = Vector3D(
        cos(dec) * cos(ra),
        cos(dec) * sin(ra),
        sin(dec)
    )
    # Transform observer's position to ECEF frame
    observer_position_ecef = earth.transform(observer_point)

    # Transform observer's position from ECEF frame to ECI frame.
    ecef_to_eci_transform = FramesFactory.getITRF(IERSConventions.IERS_2010, True).getStaticTransformTo(ECI, date)
    observer_position_eci = ecef_to_eci_transform.transformPosition(observer_position_ecef)

    # Transform direction from topocentric to ECI frame
    topocentric_to_eci_transform = observer_topo_frame.getStaticTransformTo(ECI, date)
    satellite_direction_eci = topocentric_to_eci_transform.transformVector(sat_direction_topo)

    # Scale the direction vector by the range
    satellite_relative_position_eci = satellite_direction_eci.scalarMultiply(range_in_m)

    # Combine observer's ECEF position with the satellite's relative position
    satellite_position_eci = observer_position_eci.add(satellite_relative_position_eci)

    # The satellite's position in ECEF coordinates
    
    print("Satellite Position in ECI (X, Y, Z):", satellite_position_eci)
    g = input()
    
Position = create_sat_position(
    327.850223,
    36.185382,
    38.994194,
    -104.6365,
    2.237,
    1120.4,
    earth,
    AbsoluteDate(2023, 6, 23, 4, 46, 38.674272, utc_time_scale),
)

I am expecting a position output of X=-1873.775, Y = -4683.677, Z = 4658.770 (km) but I am getting:

Satellite Position in ECI (X, Y, Z): {-2,421,665.6076448183; -5,302,058.939778156; 4,041,127.0943075754}

Thank you for your time.

2

Hi and welcome to the forum,

I believe you’ll find a solution there.

Cheers,
Romain.

3

Hi @Serrof.

Thank you for your reply. I tried to implement what was tagged as the solution in the link that you sent. I tried to implement it as shown below. I am simply trying to get a position vector right now in the ECI (J2000) frame. I am still getting values that are off. Here is my code:

# create a timescale object for UTC time
utc_time_scale = TimeScalesFactory.getUTC()

# initialise our different frames
ECI = FramesFactory.getEME2000()  # Define ECI reference frame.
ECEF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)  # Define ECEF reference frame.

# we need a BodyShape object to reference when we make our Topocentric Frame
earth = OneAxisEllipsoid(
    Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
    Constants.WGS84_EARTH_FLATTENING,
    ECEF,
)

def create_sat_position(ra, dec, lat, lon, alt, rng, earth: OneAxisEllipsoid, date: AbsoluteDate):
    """
    Returns 'Position' of observation in ECI frame.
    """
    range_in_m = rng * 1000
    alt_in_m = alt * 1000

    ra = radians(ra)
    dec = radians(dec)
    lat = radians(lat)
    lon = radians(lon)

    # Create the topocentric frame of the observer
    groundStationPoint = GeodeticPoint(lat, lon, alt_in_m)
    groundStationFrame = TopocentricFrame(earth, groundStationPoint, "Observer")

    # Convert Ra and Dec to a direction vector in the topocentric frame
    sat_direction_topo = Vector3D(
        range_in_m * cos(dec) * cos(ra),
        range_in_m * cos(dec) * sin(ra),
        range_in_m * sin(dec)
    )

    transform = groundStationFrame.getTransformTo(earth.getBodyFrame(), date)

    # Transform RADEC to ECI coordinates
    pvECI = groundStationFrame.getTransformTo(earth.getBodyFrame(), date).transformPVCoordinates(PVCoordinates(sat_direction_topo,
                         Vector3D.ZERO)); #Velocity is assumed to be zero in this example

    print(pvECI.getPosition().scalarMultiply(1/1000))


position = create_sat_position(
    327.850223,
    36.185382,
    38.994194,
    -104.6365,
    2.237,
    1120.4,
    earth,
    AbsoluteDate(2023, 6, 23, 4, 46, 38.674272, utc_time_scale),
)

I am concerned that the frame conversion is still not correct which might be what is causing this disparity. I am still confident that the date, time and measurement values are all accurately inputted into the function call. For your reference:

Expected Result: X=-1873.775, Y = -4683.677, Z = 4658.770 (km)

After implementing the solution outlined in your suggestion, I am getting:

X = -720.3461726758, Y = -5,788.190605978, Z = 4,035.4618581839

I would appreciate any advice. Thank you.

4

Hi there,

You’re transforming to ECEF here, not ECI. Just before your call to transformPVCoordinates, the call to getTransformTo must use the inertial frame as input, not earth.getBodyFrame.

NB: to build your Vector3D, you could do Vector3D(ra, dec).scalarMultiply(range_in_m)

Cheers,
Romain.

5

Hi @Serrof

I really appreciate your advice. I am still however facing issues where the output is not accurate after performing the frame conversion. I have taken your last advice into account and used the inertial frame (ECI) as the input.

def create_sat_position(ra, dec, lat, lon, alt, rng, earth: OneAxisEllipsoid, date: AbsoluteDate):
    """
    Returns 'Position' of observation in ECI frame.
    """
    range_in_m = rng * 1000
    alt_in_m = alt * 1000

    ra = radians(ra)
    dec = radians(dec)
    lat = radians(lat)
    lon = radians(lon)

    # Create the topocentric frame of the observer
    groundStationPoint = GeodeticPoint(lat, lon, alt_in_m)
    groundStationFrame = TopocentricFrame(earth, groundStationPoint, "Observer")

    # Convert Ra and Dec to a direction vector in the topocentric frame
    sat_direction_topo = Vector3D(ra, dec).scalarMultiply(range_in_m)

    # Transform RADEC to ECI coordinates
    pvECI = groundStationFrame.getTransformTo(ECI, date).transformPVCoordinates(
        PVCoordinates(sat_direction_topo, Vector3D.ZERO)
    )  # Velocity is assumed to be zero in this example

I am using the same test and am getting the old result I was getting in my first query:

Expected: X=-1873.775, Y = -4683.677, Z = 4658.770 (km)
Result: X = -2,421,665.6076448183 Y = -5,302,058.939778156, Z = 4,041,127.0943075754 (km)

Could this be an issue with the Python wrapper? I am not sure where I am going wrong with the frame conversion or how to further troubleshoot the problem. I am confident that the expected result is correct and the inputs are correctly applied. Thank you once again for your help.

6

I guess your ra and dec stand for right ascension and declination. If this is the case, the sat_direction_topo vector computed is not expressed in groundStationFrame but in some intermediate frame centered at ground station but oriented along the axes of inertial frame. A TopoCentricFrame orientation is:

  • X axis in the local horizontal plane (normal to zenith direction) and following the local parallel towards East
  • Y axis in the horizontal plane (normal to zenith direction) and following the local meridian towards North
  • Z axis towards Zenith direction

station_satellite_inert = Vector3D(ra, dec).scalarMultiply(range_in_m) will give you a station-to-satellite vector in inertial frame. You should add the earth-center-to-station vector in inertial frame to it to get the vector you need in inertial frame. You don’t need the full groundStationFrame, only the position of its origin (something like earth_station_inert = groundStationFrame.getTransformTo(ECI, date).transformPosition(Vector3D.ZERO) and then p1 = earth_station_inert.add(station_satellite_inert).

7

My bad, the link I gave you @panik was for converting AzEl, which are in a topocentric frame with “geodetic” axes. As Luc pointed out, topocentric RaDec usually use inertially oriented axes.

Cheers,
Romain

8

Thank you @luc for this detailed clarification and explanation. This seems to fix my problem. Very much appreciated.