Creating a new measurement class Python wrapper

Dear all,

I am trying to create a new measurement class in Python since the measurement is not included with Orekit measurement classes. I am using Python wrapper. I would appreciate if you could direct me to some relevant documentation to create own measurement class in python.

As far as I understand, the Python wrapper has its own set of classes. Therefore the new class has to inherit from PythonAbstractMeasurement.

class MyClass(PythonAbstractMeasurement):
         def __init__(self, date, observed, sigma, base_weight, satellite):
                  super().__init__(date, observed, sigma, base_weight, satellite)
    
    
        def theoreticalEvaluation(self, iteration, evaluation, states):
               # override the method

But unlike Java, I am unable to find the python implementations of measurement classes.

Does python support signalTimeOfFlight function?

Thank you.

Hi @niluj,

I’m not sure I understand… The Java implementation of the measurements is wrapped in Python.
There is technically no implementation of the measurements in Python.
In the __init__.py file of the Python equivalent of the estimation package (look here in @petrus.hyvonen Github repo) you’ll see that the methods of the measurements are declared but not implemented.

Yes, as long as you inherit from PythonAbstractMeasurement you shoudl have access to it.
Try self.signalTimeOfFlight(...) inside your theoreticalEvaluation method.

Hope this helps,
Maxime

Hi @MaximeJ ,

Thank for your replay. I am fairly new to Orekit and this was immensely helpful. I do have another question though.
I am trying to write my own TDOA class as to understand the process of creating a new measurement class. Following is what I have been doing.

class myTDOA(PythonAbstractMeasurement):
    def __init__(self, primaryGS,secondaryGS,date, observedValue, sigma, baseWeight, satellite):
        super().__init__(date,observedValue, sigma, baseWeight, satellite)
        self._station = gs
        self._twoWay = twoWay

data3=JArray_double(3)
data3[0] = -0.05683626922643456
data3[1] = 0.01
data3[2] = 0.1
# Orekit tdoa class
tdoa_measurement = TDOA(stations[0],stations[1],orekit_measurement_times[0],data3[0],data3[1],data3[2],sat)
# New tdoa class
tdoa_measurement2 = myTDOA(stations[0],stations[1],orekit_measurement_times[0],data3[0],data3[1],data3[2],sat)

I am getting an Argument error with the line

super().__init__(date,observedValue, sigma, baseWeight, satellite)

I think above is syntax wise correct. But it seems the last parameter should be a list of satellites based on the Java implementation. I am uncertain whether the error is coming from that.

I appreciate some insights in to this.

Thank you very much.

Hi @niluj,

satellite must be a list of Observable satellite
Try:

satellite = ObservableSatellite(0)
from java.util import ArrayList
sat_list=ArrayList()
sat_list.add(satellite)

# And then
tdoa_measurement2 = myTDOA(stations[0],stations[1],orekit_measurement_times[0],data3[0],data3[1],data3[2],sat_list)

Another option, using Java Collections class:

from java.util import Collections
sat_list = Collections.singletonList(satellite)

(this one is actually the one used in the Java constructor of the TDOA class).

Note that you can directly create the list of one ObservableSatellite when calling the super().__init__ method in the constructor of the class myTDOA.

Thank you very much!! I did not know that would work in Python.

@MaximeJ

Thank you for your suggestion. I was able to create my own TDOA class without syntax error. But I am still debugging to check whether the class calculates values correctly.

There is an offset between the TDOA estimated values from my measurement class and inbuild class.

  1. I have not selected any clock offset parameters to estimate
  2. I omitted parameterDriver
  3. I have followed the TDOA class (java) exactly, except for the getOffsetToInertial method call. I am calling the method with two argument and as a result getting a Transform object( Not fieldTransform). Made some modification to account for that and get a FieldAboluteDate
  4. Omitted the lines 193 and 194 from TDOA class
        final Gradient offset1 = primeStation.getClockOffsetDriver().getValue(nbParams, indices);
        final Gradient offset2 = secondStation.getClockOffsetDriver().getValue(nbParams, indices);

The first part of code of my class is

class MyTDOA(PythonAbstractMeasurement):
    
    def __init__(self,primeStation,secondStation,date, tdoa,sigma,base_weight, satellite):
        sat_list = Collections.singletonList(satellite)
        super().__init__(date, tdoa, sigma, base_weight, sat_list)
        self._primeStation  = primeStation
        self._secondStation = secondStation
        self._observedTdoa = tdoa

    def theoreticalEvaluation(self,iteration,evaluation,states):
        state = states[0] # space craft state
        nbParams = 6;
        zero = FieldVector3D.getZero(GradientField.getField(nbParams))
        filedType =  Gradient.variable(6, 0, 0.0).getField()
      
        pvaG = self.getCoordinates(state, 0, nbParams);
      
        primeToInert = self._primeStation.getOffsetToInertial(state.getFrame(), self.getDate())
        measurementDateG = FieldAbsoluteDate(filedType, primeToInert.getDate())
        # get PV coordinates from field transform
        primePV = primeToInert.transformPVCoordinates(TimeStampedFieldPVCoordinates(measurementDateG,zero, zero, zero))
        primePVPosition =  primePV.getPosition()
        tau1 = self.signalTimeOfFlight(pvaG, primePVPosition, measurementDateG)
        dtMtau1 = Gradient.cast_(measurementDateG.durationFrom(state.getDate()).subtract(tau1))
        emitterState = state.shiftedBy(dtMtau1.getValue())
        emitterPV = pvaG.shiftedBy(dtMtau1)
         # Repeat the similar code for second station
        #  excluded as it is lengthy
        # Following is the derivative. as per another post.
        derivatives = tdoaG.getGradient()
        jArrayArray = JArray('object')(1)
        jArrayArray[0] =  Arrays.copyOfRange(derivatives, 0, 6)
        testt = JArray('double')(Arrays.copyOfRange(derivatives, 0, 6))
        estimated.setStateDerivatives(0, jArrayArray)
        return estimated

I tried BatchLSwith both TDOA classes. As mentioned before there is an offset.
measured value : [0.0031269635433629085] date = 2022-08-15T00:00:00.000Z
estimated value from original TDOA class = [0.0031266958108103066] date = 2022-08-15T00:00:00.000Z
estimated value from my TDOA class = [0.0031269378741726195] date = 2022-08-15T00:00:00.000Z

I would like some input to what are the reasons for this. Is this due to ignoring the clock offset.?
If non of the parameters are selected for estimation what is the implication of lines 193 and 194?

Thank you.

1 Like

Hi @niluj,

That’s very interesting, and good job for writing the class.

I can’t say, did you actually define a non-zero value for the clock offset of the station(s) ?
If yes, it must come from there. If no, it shouldn’t change the value.

Selection will play on the derivatives with respect to the parameters, not on the value itself.
If the driver is selected, the derivatives with respect to it must be computed in the Batch LS so that a new estimated value can be computed after each iteration.
In both cases the values of the parameters are used to correct the measurement value at line 195/196.

That’s interesting, you’re actually closer to the measurement with your implementation than Orekit.
Is that real measurements or simulated ones ?
If the difference is not due to the clock offsets, we will be very interested in some data, like a working test case to eventually fix Orekit’s implementation.

Maxime

@MaximeJ
Thank for the reply.
I have not set clock offset values for BatchLSEstimator.
I am using simulated noisy measurements.
I am only using TDOA measurements

I got a prime meridian error similar to this post..
Interestingly, when I run a BatchLSEstimator with with inbuild TDOA class and then run another BatchLSEstimator ( I maintained two separate variables) with my TDOA class the prime meridian error disappears. Do you know why this happens ?
I change the code as suggested in the post.

Also, with my TDOA class, the BatchLSEstimator terminates after first evaluation and giving a Java run time error.The first iteration runs without any error, I can see the values with the Observer .

EXCEPTION_ACCESS_VIOLATION (0xc0000005) at pc=0x0000000056ead040, pid=16488, tid=0x0000000000002d74
Failed to write core dump. Minidumps are not enabled by default on client versions of Windows

I am running the python file within Anaconda using command prompt as Admin. The java version is

openjdk version “1.8.0_332”
OpenJDK Runtime Environment (Zulu 8.62.0.19-CA-win64) (build 1.8.0_332-b09)
OpenJDK 64-Bit Server VM (Zulu 8.62.0.19-CA-win64) (build 25.332-b09, mixed mode)
Python version is 3.10.5

I am not sure what is the cause of this error. Appreciate any ideas.

Is the error related to this post and reply by @petrus.hyvonen ?

by

Thank you.

Hi @niluj,

Do you simulate measurements with Orekit’s TDOA or your own custom TDOA class ?

You probably didn’t add the prime meridian offset as a ParameterDriver of your measurement ?
In Orekit’s TDOA, they are added in the constructor.

I have no idea, it’s the first time I’ve seen this issue. I’ve googled it and it doesn’t look like an easy one to solve, maybe changing the JDK version would help…

Hi @niluj,

Hm, don’t recognize that error no. The linked post is something triggering (discovering?) an wrapping error by PyCharm’s inspection routines in the debugger. I have not been able to isolate that one but I would hesitate to think that that is related to your error if you are running it from the normal python without debug. Do you have a longer dump, do you get a feeling if it is the java side that is crashing or the wrapper part?

Did you try if you have any changes in the error by changing JDK version?

Regards
/Petrus

Hi @petrus.hyvonen ,

Thank you for the clarification.
I have attached the error log file.
I tried the code on Linux pc as well and get the same error.
I tried to use a different openjdk version. But when I uninstall the exiting version in anaconda, it somehow removed orekit as well. I tried to manually install the orekit (Orkeit default installation commands in conda seems to install jdk1.8 also) , but was not successful with installation of a different openjdk version.
The code works with the default TDOA class. But with Python class, the code crashes. I am assuming since multiple methods are accessed via Python, the code requires more memory allocations ?

hs_err_pid17556.log (38.6 KB)

Also the code works with out error, when I create an class instance with PythonTDOA class and call .theoreticalEvaluation method seperately. (as below)

state = SpacecraftState(tle_orbit_ECI)
tdoa_test = MyTDOA(stations[0],stations[1],orekit_measurement_times[0],dd[0],50e-9,0.1,sat)
states=[state]
output_test = tdoa_test.theoreticalEvaluation(1,1,states)

The error comes with when I use Python TDOA class with BatchLSEstimator.

Update :slight_smile:
It seems I am getting the error from the self.date. I had override method and implement the getDate method in my measurement class. It seems the above error disappears. I am still doing testing and fixing the clockoffsets.

Thank you.

@petrus.hyvonen @MaximeJ

The Java run time error of my code was due to the line,

estimated = EstimatedMeasurement(self, iteration, evaluation, [emitterState],timeStampedCoordList)

I had to create a TDOA instance and pass the instance to the EstimatedMeasurement.

# Create a new instance
tdoaInstance = MyTDOA(self.getPrimeStation(),self.getSecondStation(),self.getDate(), float(self.getObservedValue()[0]),float(self.getTheoreticalStandardDeviation()[0]),float(self.getBaseWeight()[0]), self.getSatellites().get(0))
estimated = EstimatedMeasurement(tdoaInstance , iteration, evaluation, [emitterState],timeStampedCoordList)

The above error seems to resolve for now.

But the values I am getting are vastly different for the second iteration. Is is necessary to add clock offset drivers if the values are all zero for the ground stations ?
I have attached the data I get from BatchLSObserver . (observations_TDOA : data from Orekit TDOA class)
What could be the reason for the difference ?

The theoraticalEvaluation method is -

def theoreticalEvaluation(self,iteration,evaluation,states):
        state = states[0] 
        nbParams = 6
        pvaG = self.getCoordinates(state, 0, nbParams)    
        primePV = self._primeStation.getBaseFrame().getPVCoordinates(self.getDate(), state.getFrame())
        primePVFieldPosition, measurementDateG = self.getFiedStationCoordinates(primePV.getPosition(),primePV.getDate(),nbParams,0)
        tau1 = self.signalTimeOfFlight(pvaG, primePVFieldPosition, measurementDateG)
        dtMtau1 = Gradient.cast_(measurementDateG.durationFrom(state.getDate()).subtract(tau1))
        emitterState = state.shiftedBy(dtMtau1.getValue())
        emitterPV = pvaG.shiftedBy(dtMtau1)
        tau2 = tau1
        previous = tau2.getValue()                                    
        dateAt2 = emitterState.getDate().shiftedBy(previous)
        secondPV = self._secondStation.getBaseFrame().getPVCoordinates(dateAt2, state.getFrame())
        secondPVFieldPosition, measurementDateG = self.getFiedStationCoordinates(secondPV.getPosition(),secondPV.getDate(),nbParams,0)
        tau2 = Gradient.cast_(self.linkDelay(emitterPV.getPosition(), secondPVFieldPosition))
        delta = FastMath.abs(tau2.getValue() - previous)
        count = 1
        while (count < 10 and delta >= 2 * FastMath.ulp(tau2.getValue())):
        # While loop code
        tdoaG = Gradient.cast_(tau1.subtract(tau2))
        tdoa      = tdoaG.getValue()
        pv1 = primePV
        pv2 = secondPV
        timeStampedCoordList = [emitterPV.toTimeStampedPVCoordinates(), pv2 if tdoa>0 else pv1, pv1 if tdoa>0 else pv2]
        tdoaInstance = MyTDOA(self.getPrimeStation(),self.getSecondStation(),self.getDate(), float(self.getObservedValue()[0]),float(self.getTheoreticalStandardDeviation()[0]),float(self.getBaseWeight()[0]), self.getSatellites().get(0))
        estimated = EstimatedMeasurement(tdoaInstance, iteration, evaluation, [emitterState],timeStampedCoordList)
        estimated.setEstimatedValue(tdoa)
        derivatives = tdoaG.getGradient()
        jArrayArray = JArray('object')(1)
        jArrayArray[0] =  Arrays.copyOfRange(derivatives, 0, 6)
        estimated.setStateDerivatives(0, jArrayArray)
        
        return estimated





observations_TDOA.txt (2.2 MB)
observations_MyTdoa.txt (552.3 KB)

I am getting "Unable to solve , Singular problem " error due this.

Thank you.