complete surface position decoding

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:https://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:https://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:
 

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)

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
 

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
+ assert ehs.ins_vr("A000029CFFBAA11E2004727281F1") == -3648