Discrepency during converting ITRF to Topocentric AZEL between OREKIT and STK

[jlipinski]

Hi,

I have an ephemeris SP3 file with coords of one satellite in ITRF with time in TAI. I’ve constructed form it an sp3 propagator like this (code in python but I accept responses in both Python and Java):

data = DataSource("ssaja320.b23130.e23140.DG_.sp3.001")
sp3_file = sp3.SP3Parser().parse(data)
sp3_eph = sp3_file.getSatellites().get("L39")
sp3_prop = sp3_eph.getPropagator()

next I defined a little function that will help me convert python datetimes to AbsoluteDate:

def datetime_to_AbsoluteDate(row, colname):
    tmp_date = row[colname]
    y = tmp_date.year
    m = tmp_date.month
    d = tmp_date.day
    h = tmp_date.hour
    mi = tmp_date.minute
    s = tmp_date.second
    ms = float(tmp_date.microsecond/1e6)
    s = s + ms
    return AbsoluteDate(y,m,d,h,mi,s, TimeScalesFactory.getUTC())

I have a dataframe in which one of the column is named “UTC” and it contains datetimes of observations. With it I can construct new column with AbsoluteDates:

df["AbsDate"] = df.apply(lambda row: datetime_to_AbsoluteDate(row, "UTC"), axis = 1)

Then I define frame and observing station:

itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                         Constants.WGS84_EARTH_FLATTENING, 
                         itrf)

longitude = math.radians(lon)
latitude  = math.radians(lat)
altitude  = el
station = GeodeticPoint(latitude, longitude, altitude)
topo = TopocentricFrame(earth, station, "Observer")

Next I’m getting PV coordinates for my AbsDates and converting them to topocentric AZEL:

df["timestamped"] = df.apply(lambda row: sp3_prop.getPVCoordinates(row["AbsDate"], itrf), axis = 1)
df["SP3_alt_noref"] = df.apply(lambda row: math.degrees(topo.getElevation(row["timestamped"].getPosition(), itrf, row["timestamped"].getDate())), axis = 1)
df["SP3_az"] = df.apply(lambda row: math.degrees(topo.getAzimuth(row["timestamped"].getPosition(), itrf, row["timestamped"].getDate())), axis = 1)

Later I apply refraction to elevation.

I have a set of data generated with STK for comparison. In ITRF data from SP3 propagated with Orekit are consistent with data from STK (there are around 2-3 cm differences which is negligible for satellite of altitude ~1330km). Unfortunately data converted to AZEL using Orekit differs from data from STK. Residuals Orekit-STK look like this:

Does anybody has any idea if there is a mistake in my code? Parameters like latitude, longitude, elevation are the same in both STK and my code. Refraction model should be roughly the same as well, but even if there are small differences in parameters used for refraction it should only affect elevation, but the biggest differences are in azimuth.

[MaximeJ]

Hi @jlipinski,

@jlipinski

Does anybody has any idea if there is a mistake in my code? Parameters like latitude, longitude, elevation are the same in both STK and my code. Refraction model should be roughly the same as well, but even if there are small differences in parameters used for refraction it should only affect elevation, but the biggest differences are in azimuth.

I don’t know why you get these differences. It’s hard because they are rather small (30µrad at most if I’m not wrong).
Some suggestions:

• Deactivate refraction to remove this potential source of error and see if the error in elevation is better;
• Make sure that the Earth ellipsoids are exactly the same between STK and Orekit;
• Check the precision of the numbers in input for the (lat, long, alt) of the stations and make sure they are exactly the same;
• Light time delay? You’re computing a purely geometrical AZEL with Orekit. Is that the same with STK? Or is it supposed to be a “real” measurement from a station taking light time delay into account?
  If this is the latter, then you should probably use an “AngularAzel” object in Orekit instead of using the geometric function “TopocentricFrame.getAzimuth/Elevation”.

Hope this helps,
Maxime

[Serrof]

Hi there,

Maxime may very well be right, could be the light time delay. Otherwise, I’ll try different ways of building the ITRF frame.

On another note, do you know that the wrapper has built in functions to go between datetime and AbsoluteDate? It’s in the pyhelper package. Check the tutorials.

Cheers,
Romain

[Serrof]

Also, are your Earth orientation parameters up-to-date? There’s is an update.sh file in the “orekit-data” folder, check if you get different results after its execution.

Cheers,
Romain.

[jlipinski]

Thank for responses!

@MaximeJ

Deactivate refraction to remove this potential source of error and see if the error in elevation is better;

Unfortunately I cannot do that - I don’t have access to STK, I only have generated data. Without adding refraction on my side, errors are bigger.

@MaximeJ

• Make sure that the Earth ellipsoids are exactly the same between STK and Orekit;
• Check the precision of the numbers in input for the (lat, long, alt) of the stations and make sure they are exactly the same;

Everything is exactly the same.

@MaximeJ

Light time delay? You’re computing a purely geometrical AZEL with Orekit. Is that the same with STK? Or is it supposed to be a “real” measurement from a station taking light time delay into account?
If this is the latter, then you should probably use an “AngularAzel” object in Orekit instead of using the geometric function “TopocentricFrame.getAzimuth/Elevation”.

Can I get some help on how to do that? I’ve tried this approach:

generator = Generator()
sat = generator.addPropagator(sp3_prop)
azel = AngularAzElBuilder(None, groundStation, [0.], [1.], sat)
dateSelector = FixedStepSelector(0.2, TimeScalesFactory.getUTC())
scheduler = ContinuousScheduler(azel, dateSelector)
generator.addScheduler(scheduler)
t0 = datetime_to_absolutedate(df["UTC"].iloc[0])
t1 = datetime_to_absolutedate(df["UTC"].iloc[-1])
res = generator.generate(t0, t1)
az_arr = []
el_arr = []
dt_arr = []
for i in res:
    i = AngularAzEl.cast_(i)
    az = math.degrees(i.getObservedValue()[0])
    el = math.degrees(i.getObservedValue()[1])
    dt = absolutedate_to_datetime(i.date)
    az_arr.append(az)
    el_arr.append(el)
    dt_arr.append(dt)
new_df = pd.DataFrame(data = {"UTC":dt_arr, "new_AZ": az_arr, "new_ALT_noref":el_arr})

But the results are much worse:

@Serrof

Also, are your Earth orientation parameters up-to-date? There’s is an update.sh file in the “orekit-data” folder, check if you get different results after its execution.

I’ve updated orekit-data but the results stayed tha same.

[luc]

Concerning the light time delay, you could try as a rough estimate to use the position of both ground station and satellite at the same time to have a first guess of the distance, then get the delay by dividing this distance by sped of light. Then, you could fix either location of spacecraft or ground station by recomputing the position with a corrected date. Beware that currently, Orekit uses the signal arrival at ground station as the measurement date, so it is the spacecraft position that should be moved in the past with respect to the date, but I don’t know how STK works. If they chose to use the signal departure from satellite, then it is the ground station position that should be moved in the future with respect to the date. You can even do this correction a few times to improve accuracy, i.e. use a fixed point method to find the light time delay.

[jlipinski]

Thank you so much! Only the first iteration was enough to bring residuals down to ~0.1’‘-0.5’'.

[Serrof]

I’m glad you reduced the residuals “by hand”, but I don’t understand why the AngularAzElBuilder didn’t give the same results…

Romain.