add code and testing data

.gitignore

@@ -1,3 +1,5 @@
+.DS_Store
+
 # Byte-compiled / optimized / DLL files
 __pycache__/
 *.py[cod]

LICENSE

@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2015 junzis
+Copyright (c) 2015 Junzi Sun (TuDelft)
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -1,2 +1,7 @@
-# py-adsb-decoder
-Another ADS-B decoder written in python
+# A Python ADS-B Decoder
+
+This is a python ADS-B message decoder written by Junzi Sun (TuDelft). It is related with the ADS-B decode guide we have published:
+
+http://adsb-decode-guide.readthedocs.org
+
+https://github.com/junzis/adsb-decode-guide

decoder.py

@@ -0,0 +1,233 @@
+""" 
+This is a decoder of ABS-D date from Mode-S receiver. The inputs 
+    of most functions are the Hexdecial strings.
+
+Created by  : Junzi Sun (TU Delft)
+Date        : March 2015
+"""
+
+import math
+
+def hex2bin(hexstr):
+    """Convert a hexdecimal string to binary string, with zero fillings. """
+    length = len(hexstr) * 4
+    msgbin = bin(int(hexstr, 16))[2:]
+    while ((len(msgbin)) < length):
+        msgbin = '0' + msgbin
+    return msgbin
+
+def bin2int(msgbin):
+    return int(msgbin, 2)
+
+def get_df(msg):
+    """Decode Downlink Format vaule, bits 1 to 5."""
+    msgbin = hex2bin(msg)
+    return bin2int( msgbin[0:5] )
+
+def get_ca(msg):        
+    """Decode CA vaule, bits: 6 to 8."""
+    msgbin = hex2bin(msg)
+    return bin2int( msgbin[5:8] )
+
+def get_icao_addr(msg):
+    """Get the ICAO 24 bits address, bytes 3 to 8. """
+    return msg[2:8]
+
+def get_tc(msg):
+    """Get Type Code, bits 33 to 37 """
+    msgbin = hex2bin(msg)
+    return bin2int(msgbin[32:37])
+
+def get_oe_flag(msg):
+    """Check the odd/even flag. Bit 54, 0 for even, 1 for odd."""
+    msgbin = hex2bin(msg)
+    return msgbin[53]
+
+def get_alt(msg):
+    """Calculate the altitude from the message. Bit 41 to 52, Q-bit at 48"""
+    msgbin = hex2bin(msg)
+    q = msgbin[47]
+    if q:
+        n = bin2int(msgbin[40:47]+msgbin[48:52])
+        alt = n * 25 - 1000
+        return alt
+    else:
+        return None
+
+def get_cprlat(msg):
+    msgbin = hex2bin(msg)
+    return bin2int(msgbin[54:71])
+
+def get_cprlon(msg):
+    msgbin = hex2bin(msg)
+    return bin2int(msgbin[71:88])
+
+def get_position(msg0, msg1, t0, t1):
+    cprlat0 = get_cprlat(msg0)
+    cprlat1 = get_cprlat(msg1)
+    cprlon0 = get_cprlon(msg0)
+    cprlon1 = get_cprlon(msg1)
+    return cpr2position(cprlat0, cprlat1, cprlon0, cprlon1, t0, t1)
+
+def cpr2position(cprlat0, cprlat1, cprlon0, cprlon1, t0, t1):
+    '''
+     This algorithm comes from:
+     http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html.
+     
+     131072 is 2^17 since CPR latitude and longitude are encoded in 17 bits.
+    '''
+    cprlat_even = cprlat0 / 131072.0
+    cprlat_odd  = cprlat1 / 131072.0
+    cprlon_even = cprlon0 / 131072.0
+    cprlon_odd  = cprlon0 / 131072.0
+
+    air_d_lat_even = 360.0 / 60 
+    air_d_lat_odd = 360.0 / 59 
+
+    # compute latitude index 'j'
+    j = int(59 * cprlat_even - 60 * cprlat_odd + 0.5)
+
+    lat_even = float(air_d_lat_even * (j % 60 + cprlat_even))
+    lat_odd  = float(air_d_lat_odd  * (j % 59 + cprlat_odd))
+
+    if lat_even >= 270:
+        lat_even = lat_even - 360
+
+    if lat_odd >= 270:
+        lat_odd = lat_odd - 360
+
+    # check if both are in the same latidude zone, exit if not
+    if cprNL(lat_even) != cprNL(lat_odd):
+      return None
+
+    # compute ni, longitude index m, and longitude
+    if (t0 > t1):
+      ni = cprN(lat_even, 0)
+      m = math.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
+    else:
+      ni = cprN(lat_odd, 1)
+      m = math.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
+
+    if lon > 180:
+        lon = lon - 360 
+
+    return [lat, lon]
+
+
+def get_speed_heading(msg):
+    """Calculate the speed and heading."""
+    msgbin = hex2bin(msg)
+
+    v_ew_dir = bin2int(msgbin[45])
+    v_ew     = bin2int(msgbin[46:56])       # east-west velocity
+
+    v_ns_dir = bin2int(msgbin[56])
+    v_ns     = bin2int(msgbin[57:67])       # north-south velocity
+
+    v_ew = -1*v_ew if v_ew_dir else v_ew
+    v_ns = -1*v_ns if v_ns_dir else v_ns
+
+    # vr       = bin2int(msgbin[68:77])       # vertical rate
+    # vr_dir   = bin2int(msgbin[77])
+
+    speed = math.sqrt(v_ns*v_ns + v_ew*v_ew)    # unit in kts
+
+    heading = math.atan2(v_ew, v_ns)
+    heading = heading * 360.0 / (2 * math.pi)   #convert to degrees
+    heading = heading if heading >= 0 else heading + 360     # no negative val
+    return [speed, heading]
+
+def get_callsign(msg):
+    """Decode aircraft identification, aka. Callsign"""
+    
+    charset = '#ABCDEFGHIJKLMNOPQRSTUVWXYZ#####_###############0123456789######'
+    msgbin = hex2bin(msg)
+    csbin = msgbin[40:96]
+
+    cs = ''
+    cs += charset[ bin2int(csbin[0:6]) ]
+    cs += charset[ bin2int(csbin[6:12]) ]
+    cs += charset[ bin2int(csbin[12:18]) ]
+    cs += charset[ bin2int(csbin[18:24]) ]
+    cs += charset[ bin2int(csbin[24:30]) ]
+    cs += charset[ bin2int(csbin[30:36]) ]
+    cs += charset[ bin2int(csbin[36:42]) ]
+    cs += charset[ bin2int(csbin[42:48]) ]
+
+    # clean string, remove spaces and marks, if any.
+    cs = cs.replace('_', '')
+    cs = cs.replace('#', '')
+    return cs
+
+
+def cprN (lat, isodd):
+    nl = cprNL(lat) - isodd
+    return nl if nl > 1 else 1
+
+
+def cprNL(lat):
+    """Lookup table to convert the latitude to index. """
+    if lat < 0 : lat = -lat             # Table is simmetric about the equator.
+    if lat < 10.47047130 : return 59
+    if lat < 14.82817437 : return 58
+    if lat < 18.18626357 : return 57
+    if lat < 21.02939493 : return 56
+    if lat < 23.54504487 : return 55
+    if lat < 25.82924707 : return 54
+    if lat < 27.93898710 : return 53
+    if lat < 29.91135686 : return 52
+    if lat < 31.77209708 : return 51
+    if lat < 33.53993436 : return 50
+    if lat < 35.22899598 : return 49
+    if lat < 36.85025108 : return 48
+    if lat < 38.41241892 : return 47
+    if lat < 39.92256684 : return 46
+    if lat < 41.38651832 : return 45
+    if lat < 42.80914012 : return 44
+    if lat < 44.19454951 : return 43
+    if lat < 45.54626723 : return 42
+    if lat < 46.86733252 : return 41
+    if lat < 48.16039128 : return 40
+    if lat < 49.42776439 : return 39
+    if lat < 50.67150166 : return 38
+    if lat < 51.89342469 : return 37
+    if lat < 53.09516153 : return 36
+    if lat < 54.27817472 : return 35
+    if lat < 55.44378444 : return 34
+    if lat < 56.59318756 : return 33
+    if lat < 57.72747354 : return 32
+    if lat < 58.84763776 : return 31
+    if lat < 59.95459277 : return 30
+    if lat < 61.04917774 : return 29
+    if lat < 62.13216659 : return 28
+    if lat < 63.20427479 : return 27
+    if lat < 64.26616523 : return 26
+    if lat < 65.31845310 : return 25
+    if lat < 66.36171008 : return 24
+    if lat < 67.39646774 : return 23
+    if lat < 68.42322022 : return 22
+    if lat < 69.44242631 : return 21
+    if lat < 70.45451075 : return 20
+    if lat < 71.45986473 : return 19
+    if lat < 72.45884545 : return 18
+    if lat < 73.45177442 : return 17
+    if lat < 74.43893416 : return 16
+    if lat < 75.42056257 : return 15
+    if lat < 76.39684391 : return 14
+    if lat < 77.36789461 : return 13
+    if lat < 78.33374083 : return 12
+    if lat < 79.29428225 : return 11
+    if lat < 80.24923213 : return 10
+    if lat < 81.19801349 : return 9
+    if lat < 82.13956981 : return 8
+    if lat < 83.07199445 : return 7
+    if lat < 83.99173563 : return 6
+    if lat < 84.89166191 : return 5
+    if lat < 85.75541621 : return 4
+    if lat < 86.53536998 : return 3
+    if lat < 87.00000000 : return 2
+    else : return 1

run.py

@@ -0,0 +1,39 @@
+import decoder
+import re
+
+with open('data.txt', 'r') as f:
+    lines = f.readlines()
+
+for line in lines:
+    m = re.search('^(\d+\-\d+\-\d+T\d+\:\d+\:\d+\.\d+)\s*(\S*).*$', line)
+
+    time = m.group(1)
+    msg = m.group(2)
+
+    df = decoder.get_df(msg)
+    tc = decoder.get_tc(msg)
+    ca = decoder.get_ca(msg)
+
+    if df==17:
+        addr = decoder.get_icao_addr(msg)
+
+        if tc>=1 and tc<=4:
+            # aircraft identification
+            callsign = decoder.get_callsign(msg)
+            print msg, '->', 'ID:', callsign
+
+        elif tc>=9 and tc<=18:
+            # airbone postion frame
+            alt = decoder.get_alt(msg)
+            oe = decoder.get_oe_flag(msg)  # odd or even frame
+            # cprlat = decoder.get_cprlat(msg)
+            # cprlon = decoder.get_cprlon(msg)
+            print msg, '->',  'CPR odd frame' if oe else 'CPR even frame'
+
+        elif tc==19:
+            # airbone velocity frame
+            sh = decoder.get_speed_heading(msg)
+            print msg, '->', 'Speed:', sh[0], 'Heading:', sh[1]
+    else:
+        print df
+        pass