junzis / pyModeS

Python decoder for Mode S and ADS-B signals
GNU General Public License v3.0
527 stars 151 forks source link

quality of airbone positions #60

Closed wrobell closed 4 years ago

wrobell commented 4 years ago

Please consider the script below. Its output is as follows (columns are timestamp, typecode, oe flag and message)

2019-03-23 20:08:47.011677+00:00 11 1 8d4ca5415809670448d60d2fc835
2019-03-23 20:08:47.124150+00:00 19 1 8d4ca541990c9f02183c1e845693
2019-03-23 20:08:47.730598+00:00 29 0 8d4ca541ea05f9271b1c0891d457
2019-03-23 20:08:47.978328+00:00 11 0 8d4ca541580953cf23c531144f9b
message pair: (71.71367, 177.70279)
even message with ref: (53.71367, -1.18142)
odd message with ref: (53.41534, -6.16163)

I have received sequence of above messages. If I calculate the position using message pairs, then the position is far, far away from the airport near by I am getting the messages. When calculating the positions using reference location (Dublin airport) the result looks more plausible.

Is that normal?

The script

from dateutil.parser import isoparse
import pyModeS as pms

data = [
    (isoparse('2019-03-23 20:08:47.011677+00'), '8d4ca5415809670448d60d2fc835'),
    (isoparse('2019-03-23 20:08:47.124150+00'), '8d4ca541990c9f02183c1e845693'),
    (isoparse('2019-03-23 20:08:47.730598+00'), '8d4ca541ea05f9271b1c0891d457'),
    (isoparse('2019-03-23 20:08:47.978328+00'), '8d4ca541580953cf23c531144f9b'),
]

for ts, msg in data:
    print(ts, pms.adsb.typecode(msg), pms.adsb.oe_flag(msg), msg)

ts_even, msg_even = data[-1]
ts_odd, msg_odd = data[0]
pos = pms.adsb.position(msg_even, msg_odd, ts_even, ts_odd)
print('message pair:', pos)

pos = pms.adsb.position_with_ref(msg_even, 53.421389, -6.27)
print('even message with ref:', pos)

pos = pms.adsb.position_with_ref(msg_odd, 53.421389, -6.27)
print('odd message with ref:', pos)
wrobell commented 4 years ago

I also recalculated all positions (over 40mln positions) in my database with the following algorithm

pairs = defaultdict(lambda: [(None, None), (None, None)])
for ts, icao, odd, msg in data:
    pairs[icao][odd] = (ts, msg)
    (ts1, m1), (ts2, m2) = pairs[icao]
    if ts1 is not None and ts2 is not None and abs(ts1 - ts2).total_seconds() < 180:
        yield pms.adsb.position(m1, m2, ts1, ts2)

and plotted it on a map

pairs

As you can see I can provide a lot of examples where there is some issue with calculation or my data or the algorithm above.

For comparison, below all the positions recalculated using reference position

ref

junzis commented 4 years ago

I think the CPR algorithm is correctly implemented in pyModeS. However, it is possible that the decoded globally unambiguous airborne or surface position is wrong using the CPR encoding algorithm.

According to ICAO Doc 9871, we have to perform a "reasonableness test" for these decoding results.

This is something worth to be added to pyModeS. However, we might have to implement it in the streaming module, because the check needs additional information regarding the receive position, or aircraft positions decoded previously or later in time.

wrobell commented 4 years ago

Thanks for pointing me to the document. First of all, the airborne section says that the even/odd messages should be no more than 10s apart[1]. After correction of my algorithm the map looks bit better (and it yields few million less positions, which is understandable)

pairs

[1] I can see you use the same value for odd/even message in streaming module as well.

wrobell commented 4 years ago

I looked more into the document. It says that first, a pair of messages should be used to calculate position (global). Then the position should be used as reference position for next message (local) and new position becomes new reference position.

If I have not missed anything your streaming module implements above.

Using such algorithm, the map looks much better (and I get few million positions back)

pairs

Probably what's left is the reasonableness test, which is different matter. Thanks for the pointers. Closing this issue.

junzis commented 4 years ago

I looked more into the document. It says that first, a pair of messages should be used to calculate position (global). Then the position should be used as reference position for next message (local) and new position becomes new reference position.

If I have not missed anything your streaming module implements above.

Using such algorithm, the map looks much better (and I get few million positions back)

Probably what's left is the reasonableness test, which is different matter. Thanks for the pointers. Closing this issue.

Looks good! It is also possible the input data from FMS to the transponder is wrong, for example, the straight line through Africa. Normally, if you know the position of your receiver, these kinds of errors can be easily filtered out.

wrobell commented 4 years ago

I have hacked some version of reasonableness test in my scripts and got very good results. I will try to convert it into a proper function in next few days. One can follow commits in my project

https://gitlab.com/wrobell/sq1090/-/commits/master

My intention is to implement a reusable function. If I am successful, I can propose it be included in pyModeS via a pull request. Time will show. :)

wrobell commented 4 years ago

The algorithm is implemented, but for airborne positions only at the moment (see #63)

https://gitlab.com/wrobell/sq1090/-/blob/master/sq1090/pos_parser.py

The map of positions looks much better (bit different type of map this time, the darker color shows more positions) map

wrobell commented 4 years ago

I have finished implementing the algorithm for both surface and airborne positions (see also my previous comment). The reasonableness tests, as described in the ICAO document, are very important in removal of anomalies. Still, I am getting some positions, which I suspect are not within range of my receiver. Therefore, when plotting maps like below, I am removing 1% of outliers (outliers in terms of distance). If you have questions about the algorithm, please contact me via my project's GitLab: https://gitlab.com/wrobell/sq1090/.

map