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.

[jlipinski]

Hi, I don’t know if is it a good place to bring this up, as my question isn’t exactly related to OREKIT but to general astronomy, but it is somewhat of a follow up to problems I had discussed above. Right now I have some ground-based observations of LEO satellites, mostly satellites like Sentinels which have their precise ephemeris published regularly. I compare my observations with those ephemeris and I’m seeing some discrepancies - my observations are shifted in the orbital plane. Depending on the day and position from which observation takes place, the magnitude of those shifts changes. I’ve found that introducing a small time bias (0.001-0.01s) will make my observations and data from ephemeris align perfectly. Unfortunately, I can’t understand the nature of those time biases. I don’t believe that this is a hardware problem, I am triggering exposure using GPS. For data reduction, I’m correcting for atmospheric refraction and light time delay. Is there another physical effect I am forgetting about, something like aberration?

[MaximeJ]

Hi @jlipinski,

@jlipinski

For data reduction, I’m correcting for atmospheric refraction and light time delay. Is there another physical effect I am forgetting about, something like aberration?

Are you measurements AZEL or RADEC?
I don’t think we have an atmospheric refraction correction for RADEC, that’s why I’m asking.

I think modeling aberrations would be a good lead to explore if you’re using RADEC.

[jlipinski]

@MaximeJ

Are you measurements AZEL or RADEC?

I have ephemeris file in ITRF and my observations are in AZEL.
First I am calculating light travel time from satellite to observation site, then I’m shifting data from ephemeris by that time and then I am transforming ITRF to AZEL. Lastly I am applying refracion correction using ERFA (something like SOFA but on BSD license) method with parameters (temperature, pressure, humidity) measured on site. Then I am comparing ephemeris (transformed to AZEL, corrected for light travel time and refraction) with my observations.

[MaximeJ]

@jlipinski

I have ephemeris file in ITRF and my observations are in AZEL.

Ok thanks.
I thought you were using Orekit measurements package that’s why I asked.

@jlipinski

First I am calculating light travel time from satellite to observation site

Are you taking into account the small movement of the station due to Earth’s rotation during that light time delay? That cannot explain the (0.001-0.01s) time bias though.
Aberrations apply to telescope-based measurements (i.e. when you’re determining the satellite angles) by comparing to neighboring stars, given a star catalog).
Is that your case? (Although generally these measurements are given in RADEC).

[jlipinski]

@MaximeJ

Are you taking into account the small movement of the station due to Earth’s rotation during that light time delay? That cannot explain the (0.001-0.01s) time bias though.

Yes and no. I was testing for this but corrections were neglegible in comparison to my target precision (which is 5’') so I ditched it.

@MaximeJ

Aberrations apply to telescope-based measurements (i.e. when you’re determining the satellite angles) by comparing to neighboring stars, given a star catalog).
Is that your case? (Although generally these measurements are given in RADEC).

Yes, I am correcting for stellar, annual and diurnal aberration during astrometric solving of my observations. Moreover - I think that this would cause offsets in different directions, depending on satellite direction and location of observer. In my case I only have offsets in the orbital plane, colinear with satellite movement.

I can’t (with 100% certainty) say that this is not a hardware problem, but designing and constructing a testing rig to check for mili-, and micro-second precision is not an easy task.

[luc]

If time shifts of a few milliseconds order of magnitude helps reducing differences, I guess Earth Orientation Parameters may be the culprit.
Perhaps you could look what happens if you change your EOP sources for example to rapid data instead of final data?

[MaximeJ]

Thanks @luc for the hint!
@jlipinski I was thinking of two things:

1. If your measurement system is a telescope and you’re doing an astrometric reduction, then you don’t need to correct for atmospheric refraction. Removing it won’t correct the date bias though.
2. Could the date bias come from your telescope?
   I know that for the European SSA telescopes, a calibration campaign is done every 4 months using GALILEO satellites and precise ephemeris as references.
   Biases are estimated, in particular date biases, and they are in the order of magnitude of what you observe (tens to hundreds of milliseconds).

[jlipinski]

Hi,
Thanks @luc and and @MaximeJ for your replies. Sorry for late response, I had a busy month and then holiday break happened.

@luc, unfortunately changing EOP did not help.

@MaximeJ 1. Yes, we are using a kind of telescope. We are measuring in AZEL and transforming it to ITRF, that’s why (I think) we need to correct for refraction and light time delay. I was hoping that there is another effect, like aberration, that we forgot of.
2. We are suspecting that this is the case, but it is potentially easier and cheaper to firstly inspect our math and software. Unfortunately our system does not allow observations of GALILEO or other GNSS satellites. This system is calibrated only for LEO satellites, luckily there are satellites like Sentinels, Jasons, Haiyangs with quite precise ephemeris. For now we will be thoroughly testing our hardware, but detecting delays that short is neither easy nor cheap

[MaximeJ]

Hi @jlipinski,

@jlipinski

Yes, we are using a kind of telescope. We are measuring in AZEL and transforming it to ITRF, that’s why (I think) we need to correct for refraction and light time delay

If you are making the measurements by comparing to a star catalog then you don’t need to correct refraction (at least for elevations above a few degrees) but you will probably need to correct aberrations, .

@jlipinski

2. We are suspecting that this is the case, but it is potentially easier and cheaper to firstly inspect our math and software. Unfortunately our system does not allow observations of GALILEO or other GNSS satellites. This system is calibrated only for LEO satellites, luckily there are satellites like Sentinels, Jasons, Haiyangs with quite precise ephemeris. For now we will be thoroughly testing our hardware, but detecting delays that short is neither easy nor cheap

But you are already detecting these delays, right?
I think that using several Sentinels SP3 ephemeris on a few days and a batch-least square estimator using EphemerisPropagatorBuilders you could probably estimate the date bias.
We wanted to do a tutorial about this with @bcazabonne but we lack real measurements data.
The parameter to estimate would be GroundStation.clockOffsetDriver.

Once you’ve estimated the delay precisely you can use it as a correction to perform your everyday OD.

Maxime

[jlipinski]

@MaximeJ

If you are making the measurements by comparing to a star catalog then you don’t need to correct refraction (at least for elevations above a few degrees) but you will probably need to correct aberrations, .

Doesn’t Earth’s diurnal aberration total about 0.15-0.2’'? It’s at least an order of magnitude lower than what we are looking for. And is secular aberration not applicable as satellites rotate around the Sun at more or less the same speed as the Earth? Or my way of thinking is completely off? It’s first time I try to apply aberrations to fast-moving satellites, not “stationary” stars.

@MaximeJ

But you are already detecting these delays, right?
I think that using several Sentinels SP3 ephemeris on a few days and a batch-least square estimator using EphemerisPropagatorBuilders you could probably estimate the date bias.

Yes, but it is not that easy. Our telescopes do not track. They have a very wide field of view and are fixed in one position, so we have to hope that a particular satellite happens to fly through our field of view. That’s why it takes some time to calculate those offsets. Furthermore, we need to do this for multiple telescopes in a sample large enough to have statistical significance. While we are working on this, I wanted to explore another possible reasons for time bias.

[MaximeJ]

@jlipinski

Doesn’t Earth’s diurnal aberration total about 0.15-0.2’'? It’s at least an order of magnitude lower than what we are looking for.

Indeed it’s not what you’re looking for; correcting for aberrations was more a comment on the overall measurements’ configuration.
Plus I don’t think you would have such a clear time bias if it came from aberrations.

@jlipinski

And is secular aberration not applicable as satellites rotate around the Sun at more or less the same speed as the Earth?

It depends on the output frame of your astrometric reduction I think. Usually star catalog are given in ICRS, so you’ll have to correct for annual aberration to get an Earth-centered measurement.
See the astrometric transformation chain of the SOFA for reference (page 2 of this document).

@jlipinski

Yes, but it is not that easy. Our telescopes do not track. They have a very wide field of view and are fixed in one position, so we have to hope that a particular satellite happens to fly through our field of view. That’s why it takes some time to calculate those offsets. Furthermore, we need to do this for multiple telescopes in a sample large enough to have statistical significance.

Thanks for the clarifications, I understand better why it’s not as easy as I first thought