From 4bafa1de19a24277e04e7127ca292fb22f164481 Mon Sep 17 00:00:00 2001 From: junzis Date: Wed, 19 Oct 2016 16:22:28 +0200 Subject: [PATCH] complete surface position decoding --- README.rst | 25 +++++++---- pyModeS/adsb.py | 110 ++++++++++++++++++++++++++++++++++++++------- tests/test_adsb.py | 7 +++ tests/test_ehs.py | 5 ++- 4 files changed, 120 insertions(+), 27 deletions(-) diff --git a/README.rst b/README.rst index 0da0716..0157577 100644 --- a/README.rst +++ b/README.rst @@ -15,16 +15,16 @@ develop to decode the following messages: true airspeed, indicated airspeed, mach number, track angle, heading, and roll angle, etc. -A detailed manuel on Mode-S decoding is published by the author, at: +A detailed manuel on Mode-S decoding is published by the author, at: http://adsb-decode-guide.readthedocs.io Source code ----------- -Checkourt and contribute to this open source project at: +Checkourt and contribute to this open source project at: https://github.com/junzis/pyModeS -API documentation at: +API documentation at: http://pymodes.readthedocs.io Install @@ -49,7 +49,7 @@ Common function for Mode-S message: .. code:: python pms.df(msg) # Downlink Format - pms.crc(msg, encode=False) # Perform CRC or generate parity bit + pms.crc(msg, encode=False) # Perform CRC or generate parity bit pms.hex2bin(str) # Convert hexadecimal string to binary string pms.bin2int(str) # Convert binary string to integer @@ -62,15 +62,24 @@ Core functions for ADS-B decoding: pms.adsb.icao(msg) pms.adsb.callsign(msg) - pms.adsb.position(msg_even, msg_odd, t_even, t_odd) + + pms.adsb.position(msg_even, msg_odd, t_even, t_odd, lat_ref=None, lon_ref=None) + pms.adsb.airborne_position(msg_even, msg_odd, t_even, t_odd) + pms.adsb.surface_position(msg_even, msg_odd, t_even, t_odd, lat_ref, lon_ref) + pms.adsb.position_with_ref(msg, lat_ref, lon_ref) + pms.adsb.airborne_position_with_ref(msg, lat_ref, lon_ref) + pms.adsb.surface_position_with_ref(msg, lat_ref, lon_ref) + pms.adsb.altitude(msg) pms.adsb.velocity(msg) pms.adsb.speed_heading(msg) -**Hint: When you have a fix position of the aircraft or you know the -location of your receiver, it is convinent to use `position_with_ref()` method -to decode with only one position message (either odd or even)** +**Hint: When you have a fix position of the aircraft, it is convenient to +use `position_with_ref()` method to decode with only one position message +(either odd or even). This works with both airborne and surface position +messages. But the reference position shall be with in 180NM (airborne) +or 45NM (surface) of the true position.** Core functions for EHS decoding: diff --git a/pyModeS/adsb.py b/pyModeS/adsb.py index a7e9d63..0bb6fad 100644 --- a/pyModeS/adsb.py +++ b/pyModeS/adsb.py @@ -167,8 +167,8 @@ def cprlon(msg): return util.bin2int(msgbin[71:88]) -def position(msg0, msg1, t0, t1): - """Decode position from a pair of even and odd position message +def position(msg0, msg1, t0, t1, lat_ref=None, lon_ref=None): + """Decode position from a pair of even and odd position message (works with both airborne and surface position messages) Args: @@ -181,7 +181,12 @@ def position(msg0, msg1, t0, t1): (float, float): (latitude, longitude) of the aircraft """ if (5 <= typecode(msg0) <= 8 and 5 <= typecode(msg1) <= 8): - return surface_position(msg0, msg1, t0, t1) + if (not lat_ref) or (not lon_ref): + raise RuntimeError("Surface position encountered, a reference \ + position lat/lon required. Location of \ + receiver can be used.") + else: + return surface_position(msg0, msg1, t0, t1, lat_ref, lon_ref) elif (9 <= typecode(msg0) <= 18 and 9 <= typecode(msg1) <= 18): return airborne_position(msg0, msg1, t0, t1) @@ -233,17 +238,17 @@ def airborne_position(msg0, msg1, t0, t1): # compute ni, longitude index m, and longitude if (t0 > t1): - ni = max(_cprNL(lat_even), 1) - m = util.floor(cprlon_even * (_cprNL(lat_even)-1) - - cprlon_odd * _cprNL(lat_even) + 0.5) - lon = (360.0 / ni) * (m % ni + cprlon_even) lat = lat_even + nl = _cprNL(lat) + ni = max(_cprNL(lat)- 0, 1) + m = util.floor(cprlon_even * (nl-1) - cprlon_odd * nl + 0.5) + lon = (360.0 / ni) * (m % ni + cprlon_even) else: - ni = max(_cprNL(lat_odd) - 1, 1) - m = util.floor(cprlon_even * (_cprNL(lat_odd)-1) - - cprlon_odd * _cprNL(lat_odd) + 0.5) - lon = (360.0 / ni) * (m % ni + cprlon_odd) lat = lat_odd + nl = _cprNL(lat) + ni = max(_cprNL(lat) - 1, 1) + m = util.floor(cprlon_even * (nl-1) - cprlon_odd * nl + 0.5) + lon = (360.0 / ni) * (m % ni + cprlon_odd) if lon > 180: lon = lon - 360 @@ -255,7 +260,9 @@ def position_with_ref(msg, lat_ref, lon_ref): """Decode position with only one message, knowing reference nearby location, such as previously calculated location, ground station, or airport location, etc. - (works with both airborne and surface position messages) + Works with both airborne and surface position messages. + The reference position shall be with in 180NM (airborne) or 45NM (surface) + of the true position. Args: msg0 (string): even message (28 bytes hexadecimal string) @@ -279,7 +286,8 @@ def position_with_ref(msg, lat_ref, lon_ref): def airborne_position_with_ref(msg, lat_ref, lon_ref): """Decode airborne position with only one message, knowing reference nearby location, such as previously calculated location, - ground station, or airport location, etc. + ground station, or airport location, etc. The reference position shall + be with in 180NM of the true position. Args: msg (string): even message (28 bytes hexadecimal string) @@ -317,15 +325,83 @@ def airborne_position_with_ref(msg, lat_ref, lon_ref): return round(lat, 5), round(lon, 5) -def surface_position(msg0, msg1, t0, t1): - # TODO: implement surface positon - raise RuntimeError('suface position decoding to be implemented soon...') +def surface_position(msg0, msg1, t0, t1, lat_ref, lon_ref): + """Decode surface position from a pair of even and odd position message, + the lat/lon of receiver must be provided to yield the correct solution. + + Args: + msg0 (string): even message (28 bytes hexadecimal string) + msg1 (string): odd message (28 bytes hexadecimal string) + t0 (int): timestamps for the even message + t1 (int): timestamps for the odd message + lat_ref (float): latitude of the receiver + lon_ref (float): longitude of the receiver + + Returns: + (float, float): (latitude, longitude) of the aircraft + """ + + msgbin0 = util.hex2bin(msg0) + msgbin1 = util.hex2bin(msg1) + + # 131072 is 2^17, since CPR lat and lon are 17 bits each. + cprlat_even = util.bin2int(msgbin0[54:71]) / 131072.0 + cprlon_even = util.bin2int(msgbin0[71:88]) / 131072.0 + cprlat_odd = util.bin2int(msgbin1[54:71]) / 131072.0 + cprlon_odd = util.bin2int(msgbin1[71:88]) / 131072.0 + + air_d_lat_even = 90.0 / 60 + air_d_lat_odd = 90.0 / 59 + + # compute latitude index 'j' + j = util.floor(59 * cprlat_even - 60 * cprlat_odd + 0.5) + + # solution for north hemisphere + lat_even_n = float(air_d_lat_even * (j % 60 + cprlat_even)) + lat_odd_n = float(air_d_lat_odd * (j % 59 + cprlat_odd)) + + # solution for north hemisphere + lat_even_s = lat_even_n - 90.0 + lat_odd_s = lat_odd_n - 90.0 + + # chose which solution corrispondes to receiver location + lat_even = lat_even_n if lat_ref > 0 else lat_even_s + lat_odd = lat_odd_n if lat_ref > 0 else lat_odd_s + + # check if both are in the same latidude zone, rare but possible + if _cprNL(lat_even) != _cprNL(lat_odd): + return None + + # compute ni, longitude index m, and longitude + if (t0 > t1): + lat = lat_even + nl = _cprNL(lat_even) + ni = max(_cprNL(lat_even) - 0, 1) + m = util.floor(cprlon_even * (nl-1) - cprlon_odd * nl + 0.5) + lon = (90.0 / ni) * (m % ni + cprlon_even) + else: + lat = lat_odd + nl = _cprNL(lat_odd) + ni = max(_cprNL(lat_odd) - 1, 1) + m = util.floor(cprlon_even * (nl-1) - cprlon_odd * nl + 0.5) + lon = (90.0 / ni) * (m % ni + cprlon_odd) + + # four possible longitude solutions + lons = [lon, lon + 90.0, lon + 180.0, lon + 270.0] + + # the closest solution to receiver is the correct one + dls = [abs(lon_ref - l) for l in lons] + imin = min(range(4), key=dls.__getitem__) + lon = lons[imin] + + return round(lat, 5), round(lon, 5) def surface_position_with_ref(msg, lat_ref, lon_ref): """Decode surface position with only one message, knowing reference nearby location, such as previously calculated location, - ground station, or airport location, etc. + ground station, or airport location, etc. The reference position shall + be with in 45NM of the true position. Args: msg (string): even message (28 bytes hexadecimal string) diff --git a/tests/test_adsb.py b/tests/test_adsb.py index 75bc6cc..d2c1f98 100644 --- a/tests/test_adsb.py +++ b/tests/test_adsb.py @@ -37,6 +37,13 @@ def test_adsb_surface_position_with_ref(): assert pos == (-43.48564, 175.87195) +def test_adsb_surface_position(): + pos = adsb.surface_position("8CC8200A3AC8F009BCDEF2000000", + "8FC8200A3AB8F5F893096B000000", + 0, 2, + -43.496, 172.558) + assert pos == (-43.48564, 172.53942) + def test_adsb_alt(): assert adsb.altitude("8D40058B58C901375147EFD09357") == 39000 diff --git a/tests/test_ehs.py b/tests/test_ehs.py index 9b6fbcc..e549295 100644 --- a/tests/test_ehs.py +++ b/tests/test_ehs.py @@ -12,7 +12,8 @@ def test_ehs_BDS(): assert ehs.BDS("A0001839CA3800315800007448D9") == 'BDS40' assert ehs.BDS("A000139381951536E024D4CCF6B5") == 'BDS50' assert ehs.BDS("A000029CFFBAA11E2004727281F1") == 'BDS60' - assert ehs.BDS("A0281838CAE9E12FA03FFF2DDDE5") is None + assert ehs.BDS("A0281838CAE9E12FA03FFF2DDDE5") == 'BDS44' + assert ehs.BDS("A00017B0C8480030A4000024512F") is None def test_ehs_BDS20_callsign(): @@ -39,4 +40,4 @@ def test_ehs_BDS60_functions(): assert ehs.ias("A000029CFFBAA11E2004727281F1") == 336 assert ehs.mach("A000029CFFBAA11E2004727281F1") == 0.48 assert ehs.baro_vr("A000029CFFBAA11E2004727281F1") == 0 - assert ehs.ins_vr("A000029CFFBAA11E2004727281F1") == -3648 \ No newline at end of file + assert ehs.ins_vr("A000029CFFBAA11E2004727281F1") == -3648