initial commit for v2.0

This commit is contained in:
Junzi Sun
2017-12-20 17:00:02 +01:00
parent 2b1f2a5878
commit d13e1bd12e
15 changed files with 756 additions and 9 deletions

View File

@@ -1,5 +1,5 @@
The Python Mode-S Decoder
=========================
The Python Mode-S Decoder (2.0-Dev branch)
==========================================
Python library for Mode-S message decoding. Support Downlink Formats (DF) are:
@@ -28,6 +28,13 @@ Detailed manual on Mode-S decoding is published by the author, at:
http://adsb-decode-guide.readthedocs.io
New features in v2.0
---------------------
- ADS-B and EHS data streaming
- Active aircraft viewing (in terminal)
- More advanced BDS identification in Enhanced Mode-S
Source code
-----------
Checkout and contribute to this open source project at:

View File

@@ -1,6 +1,8 @@
from __future__ import absolute_import, print_function, division
from .util import *
from . import adsb
from . import ehs
from . import els
from .decoder.util import *
from .decoder import adsb
from .decoder import ehs
from .decoder import els
from .decoder import util
from .decoder import modes_common

View File

@@ -0,0 +1 @@
from __future__ import absolute_import, print_function, division

178
pyModeS/streamer/aero.py Normal file
View File

@@ -0,0 +1,178 @@
"""
Functions for aeronautics in this module
- physical quantities always in SI units
- lat,lon,course and heading in degrees
International Standard Atmosphere
p,rho,T = atmos(H) # atmos as function of geopotential altitude H [m]
a = vsound(H) # speed of sound [m/s] as function of H[m]
p = pressure(H) # calls atmos but retruns only pressure [Pa]
T = temperature(H) # calculates temperature [K]
rho = density(H) # calls atmos but retruns only pressure [Pa]
Speed conversion at altitude H[m] in ISA:
Mach = tas2mach(Vtas,H) # true airspeed (Vtas) to mach number conversion
Vtas = mach2tas(Mach,H) # true airspeed (Vtas) to mach number conversion
Vtas = eas2tas(Veas,H) # equivalent airspeed to true airspeed, H in [m]
Veas = tas2eas(Vtas,H) # true airspeed to equivent airspeed, H in [m]
Vtas = cas2tas(Vcas,H) # Vcas to Vtas conversion both m/s, H in [m]
Vcas = tas2cas(Vtas,H) # Vtas to Vcas conversion both m/s, H in [m]
Vcas = mach2cas(Mach,H) # Mach to Vcas conversion Vcas in m/s, H in [m]
Mach = cas2mach(Vcas,H) # Vcas to mach copnversion Vcas in m/s, H in [m]
"""
import numpy as np
"""Aero and geo Constants """
kts = 0.514444 # knot -> m/s
ft = 0.3048 # ft -> m
fpm = 0.00508 # ft/min -> m/s
inch = 0.0254 # inch -> m
sqft = 0.09290304 # 1 square foot
nm = 1852. # nautical mile -> m
lbs = 0.453592 # pound -> kg
g0 = 9.80665 # m/s2, Sea level gravity constant
R = 287.05287 # m2/(s2 x K), gas constant, sea level ISA
p0 = 101325. # Pa, air pressure, sea level ISA
rho0 = 1.225 # kg/m3, air density, sea level ISA
T0 = 288.15 # K, temperature, sea level ISA
gamma = 1.40 # cp/cv for air
gamma1 = 0.2 # (gamma-1)/2 for air
gamma2 = 3.5 # gamma/(gamma-1) for air
beta = -0.0065 # [K/m] ISA temp gradient below tropopause
r_earth = 6371000. # m, average earth radius
a0 = 340.293988 # m/s, sea level speed of sound ISA, sqrt(gamma*R*T0)
def atmos(H):
# H in metres
T = np.maximum(288.15 - 0.0065 * H, 216.65)
rhotrop = 1.225 * (T / 288.15)**4.256848030018761
dhstrat = np.maximum(0., H - 11000.0)
rho = rhotrop * np.exp(-dhstrat / 6341.552161)
p = rho * R * T
return p, rho, T
def temperature(H):
p, r, T = atmos(H)
return T
def pressure(H):
p, r, T = atmos(H)
return p
def density(H):
p, r, T = atmos(H)
return r
def vsound(H):
"""Speed of sound"""
T = temperature(H)
a = np.sqrt(gamma * R * T)
return a
def distance(lat1, lon1, lat2, lon2, H=0):
"""
Compute spherical distance from spherical coordinates.
For two locations in spherical coordinates
(1, theta, phi) and (1, theta', phi')
cosine( arc length ) =
sin phi sin phi' cos(theta-theta') + cos phi cos phi'
distance = rho * arc length
"""
# phi = 90 - latitude
phi1 = np.radians(90.0 - lat1)
phi2 = np.radians(90.0 - lat2)
# theta = longitude
theta1 = np.radians(lon1)
theta2 = np.radians(lon2)
cos = np.sin(phi1) * np.sin(phi2) * np.cos(theta1 - theta2) \
+ np.cos(phi1) * np.cos(phi2)
arc = np.arccos(cos)
dist = arc * (r_earth + H) # meters, radius of earth
return dist
def bearing(lat1, lon1, lat2, lon2):
lat1 = np.radians(lat1)
lon1 = np.radians(lon1)
lat2 = np.radians(lat2)
lon2 = np.radians(lon2)
x = np.sin(lon2-lon1) * np.cos(lat2)
y = np.cos(lat1) * np.sin(lat2) \
- np.sin(lat1) * np.cos(lat2) * np.cos(lon2-lon1)
initial_bearing = np.arctan2(x, y)
initial_bearing = np.degrees(initial_bearing)
bearing = (initial_bearing + 360) % 360
return bearing
# -----------------------------------------------------
# Speed conversions, altitude H all in meters
# -----------------------------------------------------
def tas2mach(Vtas, H):
"""True Airspeed to Mach number"""
a = vsound(H)
Mach = Vtas/a
return Mach
def mach2tas(Mach, H):
"""Mach number to True Airspeed"""
a = vsound(H)
Vtas = Mach*a
return Vtas
def eas2tas(Veas, H):
"""Equivalent Airspeed to True Airspeed"""
rho = density(H)
Vtas = Veas * np.sqrt(rho0/rho)
return Vtas
def tas2eas(Vtas, H):
"""True Airspeed to Equivalent Airspeed"""
rho = density(H)
Veas = Vtas * np.sqrt(rho/rho0)
return Veas
def cas2tas(Vcas, H):
"""Calibrated Airspeed to True Airspeed"""
p, rho, T = atmos(H)
qdyn = p0*((1.+rho0*Vcas*Vcas/(7.*p0))**3.5-1.)
Vtas = np.sqrt(7.*p/rho*((1.+qdyn/p)**(2./7.)-1.))
return Vtas
def tas2cas(Vtas, H):
"""True Airspeed to Calibrated Airspeed"""
p, rho, T = atmos(H)
qdyn = p*((1.+rho*Vtas*Vtas/(7.*p))**3.5-1.)
Vcas = np.sqrt(7.*p0/rho0*((qdyn/p0+1.)**(2./7.)-1.))
return Vcas
def mach2cas(Mach, H):
"""Mach number to Calibrated Airspeed"""
Vtas = mach2tas(Mach, H)
Vcas = tas2cas(Vtas, H)
return Vcas
def cas2mach(Vcas, H):
"""Calibrated Airspeed to Mach number"""
Vtas = cas2tas(Vcas, H)
Mach = tas2mach(Vtas, H)
return Mach

163
pyModeS/streamer/client.py Normal file
View File

@@ -0,0 +1,163 @@
'''
Stream beast raw data from a TCP server, convert to mode-s messages
'''
from __future__ import print_function, division
import sys
import socket
import time
from threading import Thread
if (sys.version_info > (3, 0)):
PY_VERSION = 3
else:
PY_VERSION = 2
class BaseClient(Thread):
def __init__(self, host, port):
Thread.__init__(self)
self.host = host
self.port = port
self.buffer = []
def connect(self):
while True:
try:
s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
s.settimeout(10) # 10 second timeout
s.connect((self.host, self.port))
print("Server connected - %s:%s" % (self.host, self.port))
print("collecting ADS-B messages...")
return s
except socket.error as err:
print("Socket connection error: %s. reconnecting..." % err)
time.sleep(3)
def read_beast_buffer(self):
'''
<esc> "1" : 6 byte MLAT timestamp, 1 byte signal level,
2 byte Mode-AC
<esc> "2" : 6 byte MLAT timestamp, 1 byte signal level,
7 byte Mode-S short frame
<esc> "3" : 6 byte MLAT timestamp, 1 byte signal level,
14 byte Mode-S long frame
<esc> "4" : 6 byte MLAT timestamp, status data, DIP switch
configuration settings (not on Mode-S Beast classic)
<esc><esc>: true 0x1a
<esc> is 0x1a, and "1", "2" and "3" are 0x31, 0x32 and 0x33
timestamp:
wiki.modesbeast.com/Radarcape:Firmware_Versions#The_GPS_timestamp
'''
messages_mlat = []
msg = []
i = 0
# process the buffer until the last divider <esc> 0x1a
# then, reset the self.buffer with the remainder
while i < len(self.buffer):
if (self.buffer[i:i+2] == [0x1a, 0x1a]):
msg.append(0x1a)
i += 1
elif (i == len(self.buffer) - 1) and (self.buffer[i] == 0x1a):
# special case where the last bit is 0x1a
msg.append(0x1a)
elif self.buffer[i] == 0x1a:
if i == len(self.buffer) - 1:
# special case where the last bit is 0x1a
msg.append(0x1a)
elif len(msg) > 0:
messages_mlat.append(msg)
msg = []
else:
msg.append(self.buffer[i])
i += 1
# save the reminder for next reading cycle, if not empty
if len(msg) > 0:
reminder = []
for i, m in enumerate(msg):
if (m == 0x1a) and (i < len(msg)-1):
# rewind 0x1a, except when it is at the last bit
reminder.extend([m, m])
else:
reminder.append(m)
self.buffer = [0x1a] + msg
else:
self.buffer = []
# extract messages
messages = []
for mm in messages_mlat:
msgtype = mm[0]
# print(''.join('%02X' % i for i in mm))
if msgtype == 0x32:
# Mode-S Short Message, 7 byte, 14-len hexstr
msg = ''.join('%02X' % i for i in mm[8:15])
elif msgtype == 0x33:
# Mode-S Long Message, 14 byte, 28-len hexstr
msg = ''.join('%02X' % i for i in mm[8:22])
else:
# Other message tupe
continue
if len(msg) not in [14, 28]:
# incomplete message
continue
ts = time.time()
messages.append([msg, ts])
return messages
def handle_messages(self, messages):
"""re-implement this method to handle the messages"""
for msg, t in messages:
print("%f %s" % (t, msg))
def run(self):
sock = self.connect()
while True:
try:
received = sock.recv(1024)
if PY_VERSION == 2:
received = [ord(i) for i in received]
self.buffer.extend(received)
# print(''.join(x.encode('hex') for x in self.buffer))
# process self.buffer when it is longer enough
# if len(self.buffer) < 2048:
# continue
# -- Removed!! Cause delay in low data rate scenario --
messages = self.read_beast_buffer()
if not messages:
continue
else:
self.handle_messages(messages)
time.sleep(0.001)
except Exception as e:
print("Unexpected Error:", e)
try:
sock = self.connect()
except Exception as e:
print("Unexpected Error:", e)
if __name__ == '__main__':
# for testing purpose only
host = sys.argv[1]
port = int(sys.argv[2])
client = BaseClient(host=host, port=port)
client.daemon = True
client.run()

View File

@@ -0,0 +1,95 @@
from __future__ import print_function, division
import os
import sys
import curses
import numpy as np
import time
import pyModeS as pms
from threading import Lock
from client import BaseClient
from stream import Stream
from screen import Screen
LOCK = Lock()
ADSB_MSG = []
ADSB_TS = []
EHS_MSG = []
EHS_TS = []
HOST = 'sil.lr.tudelft.nl'
PORT = 30334
LAT0 = 51.9899
LON0 = 4.3754
class ModesClient(BaseClient):
def __init__(self, host, port):
super(ModesClient, self).__init__(host, port)
def handle_messages(self, messages):
local_buffer_adsb_msg = []
local_buffer_adsb_ts = []
local_buffer_ehs_msg = []
local_buffer_ehs_ts = []
for msg, t in messages:
if len(msg) < 28: # only process long messages
continue
df = pms.df(msg)
if df == 17 or df == 18:
local_buffer_adsb_msg.append(msg)
local_buffer_adsb_ts.append(t)
elif df == 20 or df == 21:
local_buffer_ehs_msg.append(msg)
local_buffer_ehs_ts.append(t)
else:
continue
LOCK.acquire()
ADSB_MSG.extend(local_buffer_adsb_msg)
ADSB_TS.extend(local_buffer_adsb_ts)
EHS_MSG.extend(local_buffer_ehs_msg)
EHS_TS.extend(local_buffer_ehs_ts)
LOCK.release()
sys.stdout = open(os.devnull, 'w')
client = ModesClient(host=HOST, port=PORT)
client.daemon = True
client.start()
stream = Stream(lat0=LAT0, lon0=LON0)
try:
screen = Screen()
screen.daemon = True
screen.start()
while True:
if len(ADSB_MSG) > 200:
LOCK.acquire()
stream.process_raw(ADSB_TS, ADSB_MSG, EHS_TS, EHS_MSG)
ADSB_MSG = []
ADSB_TS = []
EHS_MSG = []
EHS_TS = []
LOCK.release()
acs = stream.get_aircraft()
# try:
screen.update_data(acs)
screen.update()
# except KeyboardInterrupt:
# raise
# except:
# continue
except KeyboardInterrupt:
sys.exit(0)
finally:
curses.endwin()

133
pyModeS/streamer/screen.py Normal file
View File

@@ -0,0 +1,133 @@
from __future__ import print_function, division
import curses
import numpy as np
import time
from threading import Thread
COLUMNS = ['lat', 'lon', 'alt', 'gs', 'tas', 'ias', 'mach', 'roc', 'trk', 'hdg', 't']
class Screen(Thread):
def __init__(self):
Thread.__init__(self)
self.screen = curses.initscr()
curses.noecho()
self.screen.keypad(True)
self.y = 3
self.x = 1
self.offset = 0
self.acs = {}
self.lock_icao = None
def reset_cursor_pos(self):
self.screen.move(self.y, self.x)
def update_data(self, acs):
self.acs = acs
def draw_frame(self):
self.screen.border(0)
self.screen.addstr(0, 2, "Online aircraft ('crtl+c' to exit, 'enter' to select)")
def update(self):
if len(self.acs) == 0:
return
resized = curses.is_term_resized(self.scr_h, self.scr_w)
if resized is True:
self.scr_h, self.scr_w = self.screen.getmaxyx()
self.screen.clear()
curses.resizeterm(self.scr_h, self.scr_w)
self.screen.refresh()
self.draw_frame()
row = 1
header = 'icao'
for c in COLUMNS:
c = 'updated' if c=='t' else c
header += '%10s' % c
if len(header) > self.scr_w - 2:
header = header[:self.scr_w-3] + '>'
self.screen.addstr(row, 1, header)
row +=1
self.screen.addstr(row, 1, '-'*(self.scr_w-2))
icaos = np.array(list(self.acs.keys()))
icaos = np.sort(icaos)
for row in range(3, self.scr_h - 3):
idx = row + self.offset
if idx > len(icaos) - 1:
line = ' '*(self.scr_w-2)
else:
line = ''
icao = icaos[idx]
ac = self.acs[icao]
line += icao
for c in COLUMNS:
if c == 't':
val = str(int(ac[c]))
line += '%12s' % val
else:
val = '' if ac[c] is None else ac[c]
line += '%10s' % val
if len(line) > self.scr_w - 2:
line = line[:self.scr_w-3] + '>'
if self.lock_icao == icao:
self.screen.addstr(row, 1, line, curses.A_STANDOUT)
elif row == self.y:
self.screen.addstr(row, 1, line, curses.A_BOLD)
else:
self.screen.addstr(row, 1, line)
self.screen.addstr(self.scr_h-3, 1, '-'*(self.scr_w-2))
total_page = len(icaos) // (self.scr_h - 4) + 1
current_page = self.offset // (self.scr_h - 4) + 1
self.screen.addstr(self.scr_h-2, 1, '(%d / %d)' % (current_page, total_page))
self.reset_cursor_pos()
def run(self):
self.draw_frame()
self.scr_h, self.scr_w = self.screen.getmaxyx()
while True:
c = self.screen.getch()
if c == curses.KEY_HOME:
self.x = 1
self.y = 1
elif c == curses.KEY_NPAGE:
offset_intent = self.offset + (self.scr_h - 4)
if offset_intent < len(self.acs) - 5:
self.offset = offset_intent
elif c == curses.KEY_PPAGE:
offset_intent = self.offset - (self.scr_h - 4)
if offset_intent > 0:
self.offset = offset_intent
else:
self.offset = 0
elif c == curses.KEY_DOWN :
y_intent = self.y + 1
if y_intent < self.scr_h - 3:
self.y = y_intent
elif c == curses.KEY_UP:
y_intent = self.y - 1
if y_intent > 2:
self.y = y_intent
elif c == curses.KEY_ENTER or c == 10 or c == 13:
self.lock_icao = (self.screen.instr(self.y, 1, 6)).decode()

169
pyModeS/streamer/stream.py Normal file
View File

@@ -0,0 +1,169 @@
from __future__ import print_function, division
import numpy as np
import time
import pyModeS as pms
class Stream():
def __init__(self, lat0, lon0):
self.acs = dict()
self.cache_new_acs = False
self.__new_acs = set()
self.lat0 = lat0
self.lon0 = lon0
self.t = 0
self.cache_timeout = 60 # seconds
def process_raw(self, adsb_ts, adsb_msgs, ehs_ts, ehs_msgs, tnow=None):
"""process a chunk of adsb and ehs messages recieved in the same
time period.
"""
if tnow is None:
tnow = time.time()
self.t = tnow
local_updated_acs_buffer = []
# process adsb message
for t, msg in zip(adsb_ts, adsb_msgs):
icao = pms.adsb.icao(msg)
tc = pms.adsb.typecode(msg)
if icao not in self.acs:
self.acs[icao] = {
'lat': None,
'lon': None,
'alt': None,
'gs': None,
'trk': None,
'roc': None,
'tas': None,
'ias': None,
'mach': None,
'hdg': None,
}
self.acs[icao]['t'] = t
if 1 <= tc <= 4:
self.acs[icao]['callsign'] = pms.adsb.callsign(msg)
if (5 <= tc <= 8) or (tc == 19):
vdata = pms.adsb.velocity(msg)
if vdata is None:
continue
spd, trk, roc, tag = vdata
if tag != 'GS':
continue
self.acs[icao]['gs'] = spd
self.acs[icao]['trk'] = trk
self.acs[icao]['roc'] = roc
self.acs[icao]['tv'] = t
if (5 <= tc <= 18):
oe = pms.adsb.oe_flag(msg)
self.acs[icao][oe] = msg
self.acs[icao]['t'+str(oe)] = t
if ('tpos' in self.acs[icao]) and (t - self.acs[icao]['tpos'] < 180):
# use single message decoding
rlat = self.acs[icao]['lat']
rlon = self.acs[icao]['lon']
latlon = pms.adsb.position_with_ref(msg, rlat, rlon)
elif ('t0' in self.acs[icao]) and ('t1' in self.acs[icao]) and \
(abs(self.acs[icao]['t0'] - self.acs[icao]['t1']) < 10):
# use multi message decoding
try:
latlon = pms.adsb.position(
self.acs[icao][0],
self.acs[icao][1],
self.acs[icao]['t0'],
self.acs[icao]['t1'],
self.lat0, self.lon0
)
except:
# mix of surface and airborne position message
continue
else:
latlon = None
if latlon is not None:
self.acs[icao]['tpos'] = t
self.acs[icao]['lat'] = latlon[0]
self.acs[icao]['lon'] = latlon[1]
self.acs[icao]['alt'] = pms.adsb.altitude(msg)
local_updated_acs_buffer.append(icao)
# process ehs message
for t, msg in zip(ehs_ts, ehs_msgs):
icao = pms.ehs.icao(msg)
if icao not in self.acs:
continue
bds = pms.ehs.BDS(msg)
if bds == 'BDS50':
tas = pms.ehs.tas50(msg)
if tas:
self.acs[icao]['t50'] = t
self.acs[icao]['tas'] = tas
elif bds == 'BDS60':
ias = pms.ehs.ias60(msg)
hdg = pms.ehs.hdg60(msg)
mach = pms.ehs.mach60(msg)
if ias or hdg or mach:
self.acs[icao]['t60'] = t
if ias:
self.acs[icao]['ias'] = ias
if hdg:
self.acs[icao]['hdg'] = hdg
if mach:
self.acs[icao]['mach'] = mach
# clear up old data
for icao in list(self.acs.keys()):
if self.t - self.acs[icao]['t'] > self.cache_timeout:
del self.acs[icao]
continue
if self.cache_new_acs:
self.add_new_aircraft(local_updated_acs_buffer)
return
def get_aircraft(self):
"""all aircraft that are stored in memeory"""
acs = self.acs
icaos = list(acs.keys())
for icao in icaos:
if acs[icao]['lat'] is None:
acs.pop(icao)
return acs
def add_new_aircraft(self, acs):
"""add new aircraft to the list"""
self.__new_acs.update(acs)
return
def get_new_aircraft(self):
"""update aircraft from last iteration"""
newacs = dict()
for ac in self.__new_acs:
newacs[ac] = self.acs[ac]
return newacs
def reset_new_aircraft(self):
"""reset the updated icao buffer once been read"""
self.__new_acs = set()

View File

@@ -30,7 +30,7 @@ setup(
# Versions should comply with PEP440. For a discussion on single-sourcing
# the version across setup.py and the project code, see
# https://packaging.python.org/en/latest/single_source_version.html
version='1.2.0',
version='2.0-dev',
description='Python Mode-S Decoder',
long_description=long_description,
@@ -63,7 +63,6 @@ setup(
# Specify the Python versions you support here. In particular, ensure
# that you indicate whether you support Python 2, Python 3 or both.
'Programming Language :: Python :: 2',
'Programming Language :: Python :: 2.6',
'Programming Language :: Python :: 2.7',
'Programming Language :: Python :: 3',
'Programming Language :: Python :: 3.3',

View File

@@ -1,6 +1,6 @@
[tox]
toxworkdir=/tmp/tox
envlist = py26,py27,py35
envlist = py27,py35
[testenv]
deps=pytest
commands=py.test