diff --git a/.idea/remote-mappings.xml b/.idea/remote-mappings.xml
new file mode 100644
index 0000000..31ee189
--- /dev/null
+++ b/.idea/remote-mappings.xml
@@ -0,0 +1,16 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
\ No newline at end of file
diff --git a/basic.py b/basic.py
index 969b2b5..1cd99f5 100644
--- a/basic.py
+++ b/basic.py
@@ -38,7 +38,7 @@ class ADSB:
This class defines fundamental constants and is not supposed to be instantiated
"""
- VERSION = "2.0"
+ VERSION = "2.1"
# basic constants
MODES_SIGMIN = 0
@@ -69,6 +69,8 @@ class ADSB:
METER_PER_FOOT = 0.3048
KPH_PER_KNOT = 1.852
+ MAX_17_BITS = float(2**17)
+
# Downlink formats
DF_SHORT_AIR2AIR_SURVEILLANCE_0 = "0"
DF_UNKNOWN_1 = "1"
diff --git a/cpr.c b/cpr.c
new file mode 100644
index 0000000..d470df2
--- /dev/null
+++ b/cpr.c
@@ -0,0 +1,318 @@
+#include
+#include
+
+//
+//=========================================================================
+//
+// Always positive MOD operation, used for CPR decoding.
+//
+static int cprModInt(int a, int b) {
+ int res = a % b;
+ if (res < 0) res += b;
+ return res;
+}
+
+static double cprModDouble(double a, double b) {
+ double res = fmod(a, b);
+ if (res < 0) res += b;
+ return res;
+}
+
+//
+//=========================================================================
+//
+// The NL function uses the precomputed table from 1090-WP-9-14
+//
+static int cprNLFunction(double lat) {
+ 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;
+}
+//
+//=========================================================================
+//
+static int cprNFunction(double lat, int fflag) {
+ int nl = cprNLFunction(lat) - (fflag ? 1 : 0);
+ if (nl < 1) nl = 1;
+ return nl;
+}
+//
+//=========================================================================
+//
+static double cprDlonFunction(double lat, int fflag, int surface) {
+ return (surface ? 90.0 : 360.0) / cprNFunction(lat, fflag);
+}
+//
+//=========================================================================
+//
+// This algorithm comes from:
+// http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html.
+//
+// A few remarks:
+// 1) 131072 is 2^17 since CPR latitude and longitude are encoded in 17 bits.
+//
+int decodeCPRairborne(int even_cprlat, int even_cprlon,
+ int odd_cprlat, int odd_cprlon,
+ int fflag,
+ double *out_lat, double *out_lon)
+{
+ double AirDlat0 = 360.0 / 60.0;
+ double AirDlat1 = 360.0 / 59.0;
+ double lat0 = even_cprlat;
+ double lat1 = odd_cprlat;
+ double lon0 = even_cprlon;
+ double lon1 = odd_cprlon;
+
+ double rlat, rlon;
+
+ // Compute the Latitude Index "j"
+ int j = (int) floor(((59*lat0 - 60*lat1) / 131072) + 0.5);
+ double rlat0 = AirDlat0 * (cprModInt(j,60) + lat0 / 131072);
+ double rlat1 = AirDlat1 * (cprModInt(j,59) + lat1 / 131072);
+
+ if (rlat0 >= 270) rlat0 -= 360;
+ if (rlat1 >= 270) rlat1 -= 360;
+
+ // Check to see that the latitude is in range: -90 .. +90
+ if (rlat0 < -90 || rlat0 > 90 || rlat1 < -90 || rlat1 > 90)
+ return (-2); // bad data
+
+ // Check that both are in the same latitude zone, or abort.
+ if (cprNLFunction(rlat0) != cprNLFunction(rlat1))
+ return (-1); // positions crossed a latitude zone, try again later
+
+ // Compute ni and the Longitude Index "m"
+ if (fflag) { // Use odd packet.
+ int ni = cprNFunction(rlat1,1);
+ int m = (int) floor((((lon0 * (cprNLFunction(rlat1)-1)) -
+ (lon1 * cprNLFunction(rlat1))) / 131072.0) + 0.5);
+ rlon = cprDlonFunction(rlat1, 1, 0) * (cprModInt(m, ni)+lon1/131072);
+ rlat = rlat1;
+ } else { // Use even packet.
+ int ni = cprNFunction(rlat0,0);
+ int m = (int) floor((((lon0 * (cprNLFunction(rlat0)-1)) -
+ (lon1 * cprNLFunction(rlat0))) / 131072) + 0.5);
+ rlon = cprDlonFunction(rlat0, 0, 0) * (cprModInt(m, ni)+lon0/131072);
+ rlat = rlat0;
+ }
+
+ // Renormalize to -180 .. +180
+ rlon -= floor( (rlon + 180) / 360 ) * 360;
+
+ *out_lat = rlat;
+ *out_lon = rlon;
+
+ return 0;
+}
+
+int decodeCPRsurface(double reflat, double reflon,
+ int even_cprlat, int even_cprlon,
+ int odd_cprlat, int odd_cprlon,
+ int fflag,
+ double *out_lat, double *out_lon)
+{
+ double AirDlat0 = 90.0 / 60.0;
+ double AirDlat1 = 90.0 / 59.0;
+ double lat0 = even_cprlat;
+ double lat1 = odd_cprlat;
+ double lon0 = even_cprlon;
+ double lon1 = odd_cprlon;
+ double rlon, rlat;
+
+ // Compute the Latitude Index "j"
+ int j = (int) floor(((59*lat0 - 60*lat1) / 131072) + 0.5);
+ double rlat0 = AirDlat0 * (cprModInt(j,60) + lat0 / 131072);
+ double rlat1 = AirDlat1 * (cprModInt(j,59) + lat1 / 131072);
+
+ // Pick the quadrant that's closest to the reference location -
+ // this is not necessarily the same quadrant that contains the
+ // reference location.
+ //
+ // There are also only two valid quadrants: -90..0 and 0..90;
+ // no correct message would try to encoding a latitude in the
+ // ranges -180..-90 and 90..180.
+ //
+ // If the computed latitude is more than 45 degrees north of
+ // the reference latitude (using the northern hemisphere
+ // solution), then the southern hemisphere solution will be
+ // closer to the refernce latitude.
+ //
+ // e.g. reflat=0, rlat=44, use rlat=44
+ // reflat=0, rlat=46, use rlat=46-90 = -44
+ // reflat=40, rlat=84, use rlat=84
+ // reflat=40, rlat=86, use rlat=86-90 = -4
+ // reflat=-40, rlat=4, use rlat=4
+ // reflat=-40, rlat=6, use rlat=6-90 = -84
+
+ // As a special case, -90, 0 and +90 all encode to zero, so
+ // there's a little extra work to do there.
+
+ if (rlat0 == 0) {
+ if (reflat < -45)
+ rlat0 = -90;
+ else if (reflat > 45)
+ rlat0 = 90;
+ } else if ((rlat0 - reflat) > 45) {
+ rlat0 -= 90;
+ }
+
+ if (rlat1 == 0) {
+ if (reflat < -45)
+ rlat1 = -90;
+ else if (reflat > 45)
+ rlat1 = 90;
+ } else if ((rlat1 - reflat) > 45) {
+ rlat1 -= 90;
+ }
+
+ // Check to see that the latitude is in range: -90 .. +90
+ if (rlat0 < -90 || rlat0 > 90 || rlat1 < -90 || rlat1 > 90)
+ return (-2); // bad data
+
+ // Check that both are in the same latitude zone, or abort.
+ if (cprNLFunction(rlat0) != cprNLFunction(rlat1))
+ return (-1); // positions crossed a latitude zone, try again later
+
+ // Compute ni and the Longitude Index "m"
+ if (fflag) { // Use odd packet.
+ int ni = cprNFunction(rlat1,1);
+ int m = (int) floor((((lon0 * (cprNLFunction(rlat1)-1)) -
+ (lon1 * cprNLFunction(rlat1))) / 131072.0) + 0.5);
+ rlon = cprDlonFunction(rlat1, 1, 1) * (cprModInt(m, ni)+lon1/131072);
+ rlat = rlat1;
+ } else { // Use even packet.
+ int ni = cprNFunction(rlat0,0);
+ int m = (int) floor((((lon0 * (cprNLFunction(rlat0)-1)) -
+ (lon1 * cprNLFunction(rlat0))) / 131072) + 0.5);
+ rlon = cprDlonFunction(rlat0, 0, 1) * (cprModInt(m, ni)+lon0/131072);
+ rlat = rlat0;
+ }
+
+ // Pick the quadrant that's closest to the reference location -
+ // this is not necessarily the same quadrant that contains the
+ // reference location. Unlike the latitude case, all four
+ // quadrants are valid.
+
+ // if reflon is more than 45 degrees away, move some multiple of 90 degrees towards it
+ rlon += floor( (reflon - rlon + 45) / 90 ) * 90; // this might move us outside (-180..+180), we fix this below
+
+ // Renormalize to -180 .. +180
+ rlon -= floor( (rlon + 180) / 360 ) * 360;
+
+ *out_lat = rlat;
+ *out_lon = rlon;
+ return 0;
+}
+
+//
+//=========================================================================
+//
+// This algorithm comes from:
+// 1090-WP29-07-Draft_CPR101 (which also defines decodeCPR() )
+//
+// Despite what the earlier comment here said, we should *not* be using trunc().
+// See Figure 5-5 / 5-6 and note that floor is applied to (0.5 + fRP - fEP), not
+// directly to (fRP - fEP). Eq 38 is correct.
+//
+int decodeCPRrelative(double reflat, double reflon,
+ int cprlat, int cprlon,
+ int fflag, int surface,
+ double *out_lat, double *out_lon)
+{
+ double AirDlat;
+ double AirDlon;
+ double fractional_lat = cprlat / 131072.0;
+ double fractional_lon = cprlon / 131072.0;
+ double rlon, rlat;
+ int j,m;
+
+ AirDlat = (surface ? 90.0 : 360.0) / (fflag ? 59.0 : 60.0);
+
+ // Compute the Latitude Index "j"
+ j = (int) (floor(reflat/AirDlat) +
+ floor(0.5 + cprModDouble(reflat, AirDlat)/AirDlat - fractional_lat));
+ rlat = AirDlat * (j + fractional_lat);
+ if (rlat >= 270) rlat -= 360;
+
+ // Check to see that the latitude is in range: -90 .. +90
+ if (rlat < -90 || rlat > 90) {
+ return (-1); // Time to give up - Latitude error
+ }
+
+ // Check to see that answer is reasonable - ie no more than 1/2 cell away
+ if (fabs(rlat - reflat) > (AirDlat/2)) {
+ return (-1); // Time to give up - Latitude error
+ }
+
+ // Compute the Longitude Index "m"
+ AirDlon = cprDlonFunction(rlat, fflag, surface);
+ m = (int) (floor(reflon/AirDlon) +
+ floor(0.5 + cprModDouble(reflon, AirDlon)/AirDlon - fractional_lon));
+ rlon = AirDlon * (m + fractional_lon);
+ if (rlon > 180) rlon -= 360;
+
+ // Check to see that answer is reasonable - ie no more than 1/2 cell away
+ if (fabs(rlon - reflon) > (AirDlon/2))
+ return (-1); // Time to give up - Longitude error
+
+ *out_lat = rlat;
+ *out_lon = rlon;
+ return (0);
+}
\ No newline at end of file
diff --git a/emitter.py b/emitter.py
index 9a54edc..238562a 100644
--- a/emitter.py
+++ b/emitter.py
@@ -12,8 +12,21 @@
def get_msg(message):
- sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
- sock.connect((cfg_server_address, cfg_server_port))
+ try:
+ sock = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
+ sock.setsockopt(socket.SOL_SOCKET, socket.SO_REUSEADDR, 1)
+ except socket.error as msg:
+ print "Sock err (1): {}".format(msg)
+ sock = None
+ return None
+
+ try:
+ sock.connect((cfg_server_address, cfg_server_port))
+ except socket.error as msg:
+ print "Sock err (2): {}".format(msg)
+ sock = None
+ return None
+
data = []
try:
sock.sendall(message)
@@ -24,6 +37,8 @@ def get_msg(message):
else:
data.append(stream)
data = ''.join(data)
+ except socket.error as msg:
+ print "Sock err (3): {}".format(msg)
finally:
sock.close()
return data
diff --git a/map.html b/map.html
index 073e1cb..2c48c39 100644
--- a/map.html
+++ b/map.html
@@ -61,7 +61,7 @@
function initMap() {
Map = new google.maps.Map(document.getElementById('map_canvas'), {
center: {lat: CenterLat, lng: CenterLon},
- zoom: 5,
+ zoom: 8,
mapTypeId: google.maps.MapTypeId.ROADMAP
});
window.setInterval(function () {
@@ -74,17 +74,17 @@
var r = 255, g = 255, b = 0;
var maxalt = 40000;
/* Max altitude in the average case */
- var invalt = maxalt - parseInt(plane['altitude']);
+ var invalt = maxalt - parseFloat(plane['altitude']);
var selected = (Selected == plane['ICAO24']);
if (invalt < 0) invalt = 0;
- b = parseInt(255 / maxalt * invalt);
+ b = parseFloat(255 / maxalt * invalt);
return {
strokeWeight: (selected ? 2 : 1),
path: google.maps.SymbolPath.FORWARD_CLOSED_ARROW,
scale: 5,
fillColor: 'rgb(' + r + ',' + g + ',' + b + ')',
fillOpacity: 0.9,
- rotation: parseInt(plane['heading'])
+ rotation: parseFloat(plane['heading'])
};
}
function selectPlane() {
@@ -133,20 +133,20 @@
var myplane = Planes[plane['ICAO24']];
marker = myplane.marker;
var icon = marker.getIcon();
- var newpos = new google.maps.LatLng(parseInt(plane['latitude']), parseInt(plane['longitude']));
+ var newpos = new google.maps.LatLng(parseFloat(plane['latitude']), parseFloat(plane['longitude']));
marker.setPosition(newpos);
marker.setIcon(getIconForPlane(plane));
- myplane.altitude = parseInt(plane['altitude']);
- myplane.speed = parseInt(plane['velocity']);
- myplane.lat = parseInt(plane['latitude']);
- myplane.lon = parseInt(plane['longitude']);
- myplane.track = parseInt(plane['heading']);
+ myplane.altitude = parseFloat(plane['altitude']);
+ myplane.speed = parseFloat(plane['velocity']);
+ myplane.lat = parseFloat(plane['latitude']);
+ myplane.lon = parseFloat(plane['longitude']);
+ myplane.track = parseFloat(plane['heading']);
myplane.flight = plane['call_sign'];
if (myplane['ICAO24'] == Selected)
refreshSelectedInfo();
} else {
marker = new google.maps.Marker({
- position: new google.maps.LatLng(parseInt(plane['latitude']), parseInt(plane['longitude'])),
+ position: new google.maps.LatLng(parseFloat(plane['latitude']), parseFloat(plane['longitude'])),
map: Map,
icon: getIconForPlane(plane)
});
diff --git a/radar.py b/radar.py
index 94e63cc..575101a 100644
--- a/radar.py
+++ b/radar.py
@@ -240,11 +240,8 @@ def _blip_add(self, msg):
else:
self.blips[icao] = {'msg': msg, 'timestamp': time.time(), 'count': 1}
- # Note, we need to decode for latitude/longitude after we have updated the message as the algorithm
- # is dependent on 2 frames; odd and even
- if not self.blips[icao]['msg'].decodeCPR():
- pass # decodeCPR_relative is not validated yet and will likely yield wrong values
- # self.blips[icao]['msg'].decodeCPR_relative()
+ if not self.blips[icao]['msg'].decodeCPR_relative():
+ self.blips[icao]['msg'].decodeCPR()
self.lock.release()
if self.cfg_verbose_logging:
diff --git a/spots_config.json b/spots_config.json
index 9425339..3dd8ddb 100644
--- a/spots_config.json
+++ b/spots_config.json
@@ -14,6 +14,6 @@
"log file": "spots.log",
"log max bytes": 1048576,
"log backup count": 10,
- "spots server address": "localhost",
+ "spots server address": "",
"spots server port": 5051
}
\ No newline at end of file
diff --git a/squitter.py b/squitter.py
index b3654ee..921eb98 100644
--- a/squitter.py
+++ b/squitter.py
@@ -1,6 +1,7 @@
import basic
import math
import logging
+import time
__author__ = 'Wolfrax'
@@ -215,9 +216,8 @@ def __init__(self):
self.odd_raw_longitude = 0
self.even_raw_latitude = 0
self.even_raw_longitude = 0
- self.odd_pos = False
- self.even_pos = False
- self.even_then_odd_order = False
+ self.even_time = 0
+ self.odd_time = 0
self.on_ground = False
self.logger = logging.getLogger('spots.squitter')
@@ -278,13 +278,11 @@ def update(self, msg):
self.odd_raw_latitude = msg.odd_raw_latitude if msg.odd_raw_latitude != 0 else self.odd_raw_latitude
self.odd_raw_longitude = msg.odd_raw_longitude if msg.odd_raw_longitude != 0 else self.odd_raw_longitude
-
self.even_raw_latitude = msg.even_raw_latitude if msg.even_raw_latitude != 0 else self.even_raw_latitude
self.even_raw_longitude = msg.even_raw_longitude if msg.even_raw_longitude != 0 else self.even_raw_longitude
- self.odd_pos = msg.odd_pos if msg.odd_pos else self.odd_pos
- self.even_pos = msg.even_pos if msg.even_pos else self.even_pos
- self.even_then_odd_order = msg.even_then_odd_order if msg.even_then_odd_order else self.even_then_odd_order
+ self.even_time = msg.even_time if msg.even_time != 0.0 else self.even_time
+ self.odd_time = msg.odd_time if msg.odd_time != 0.0 else self.odd_time
def get_downlink_format(self):
return self['downlink_format']
@@ -292,124 +290,98 @@ def get_downlink_format(self):
def _get_msg_byte(self, byte_nr):
return (self.msg >> (self.no_of_bits - (byte_nr + 1) * 8)) & 0xFF
- def decodeCPR(self):
- # Basic algorithm: http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html
-
- if self.odd_pos is False or self.even_pos is False:
- return False
-
- longitude = 0.0
-
- air_dlat = 90.0 if self.on_ground else 360.0
-
- air_dlat_0 = air_dlat / 60.0
- air_dlat_1 = air_dlat / 59.0
- lat0 = self.even_raw_latitude
- lat1 = self.odd_raw_latitude
- lon0 = self.even_raw_longitude
- lon1 = self.odd_raw_longitude
-
- j = int(math.floor(((59 * lat0 - 60 * lat1) / 131072.0) + 0.5))
-
- rlat0 = air_dlat_0 * (j % 60 + lat0 / 131072.0)
- rlat1 = air_dlat_1 * (j % 59 + lat1 / 131072.0)
-
- if self.on_ground:
- surface_rlat = self['latitude'] if self['latitude'] != 0 else self.cfg_latitude
- surface_rlon = self['longitude'] if self['longitude'] != 0 else self.cfg_longitude
- rlat0 += math.floor(surface_rlat / 90.0) * 90.0
- rlat1 += math.floor(surface_rlat / 90.0) * 90.0
- longitude = 90 + math.floor(surface_rlon / 90.0) * 90.0
- else: # Adjust if we are on southern hemisphere by substracting 360 from latitude
- if rlat0 >= 270:
- rlat0 -= 360
- if rlat1 >= 270:
- rlat1 -= 360
-
- if rlat0 < -90 or rlat0 > 90 or rlat1 < -90 or rlat1 > 90:
- return False
- if CPR_NL(rlat0) != CPR_NL(rlat1):
- return False
-
- if not self.on_ground:
- if self.even_then_odd_order:
- # Decode using odd as the latest message
- ni = max(CPR_NL(rlat1) - 1, 1)
- m = int(math.floor((((lon0 * (CPR_NL(rlat1) - 1)) - (lon1 * CPR_NL(rlat1))) / 131072.0) + 0.5))
- longitude = (360.0 / ni) * ((m % ni) + lon1 / 131072.0)
- else:
- # Decode using even as the latest message
- ni = max(CPR_NL(rlat0) - 1, 1)
- m = int(math.floor((((lon0 * (CPR_NL(rlat0) - 1)) - (lon1 * CPR_NL(rlat0))) / 131072.0) + 0.5))
- longitude = (360.0 / ni) * ((m % ni) + lon0 / 131072.0)
+ def decodeCPR_relative(self):
+ # Basic algorithms
+ # https://adsb-decode-guide.readthedocs.io/en/latest/content/airborne-position.html
- if longitude > 180:
- longitude -= 360
- latitude = rlat1
+ if self.odd_time == 0 and self.even_time == 0:
+ return False # Got to have at least one message
+
+ if self.odd_time == 0: # Even message
+ d_lat = 360.0 / 60.0
+ lat_cpr = self.even_raw_latitude / self.MAX_17_BITS
+ lon_cpr = self.even_raw_longitude / self.MAX_17_BITS
+ else: # Odd message
+ d_lat = 360.0 / 59.0
+ lat_cpr = self.odd_raw_latitude / self.MAX_17_BITS
+ lon_cpr = self.odd_raw_longitude / self.MAX_17_BITS
+
+ j = math.floor(self.cfg_latitude / d_lat) + \
+ math.floor((self.cfg_latitude % d_lat) / d_lat - lat_cpr + 0.5) # latitude index
+
+ latitude = d_lat * (j + lat_cpr)
+
+ if CPR_NL(latitude) == 0:
+ d_lon = 360.0
+ else:
+ d_lon = 360.0 / CPR_NL(latitude)
+
+ m = math.floor(self.cfg_longitude / d_lon) + math.floor((self.cfg_longitude % d_lon) / d_lon - lon_cpr + 0.5)
+
+ longitude = d_lon * (m + lon_cpr)
self.data['latitude'] = str(round(latitude, 3)) if latitude != 0.0 else ""
self.data['longitude'] = str(round(longitude, 3)) if longitude != 0.0 else ""
- # Reset all flags
- self.even_pos = False
- self.odd_pos = False
+ return True
- self.odd_raw_latitude = 0
- self.odd_raw_longitude = 0
- self.even_raw_latitude = 0
- self.even_raw_longitude = 0
+ def decodeCPR(self):
+ # Basic algorithms
+ # http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html
+ # https://adsb-decode-guide.readthedocs.io/en/latest/content/airborne-position.html
- self.even_then_odd_order = False
+ if self.odd_time == 0 or self.even_time == 0:
+ return False # we need both even + odd messages
- return True
+ cpr_lat_even = self.even_raw_latitude / self.MAX_17_BITS
+ cpr_lon_even = self.even_raw_longitude / self.MAX_17_BITS
+ cpr_lat_odd = self.odd_raw_latitude / self.MAX_17_BITS
+ cpr_lon_odd = self.odd_raw_longitude / self.MAX_17_BITS
- def decodeCPR_relative(self):
- if self.odd_pos is False or self.even_pos is False:
- return False
+ j = math.floor(59 * cpr_lat_even - 60 * cpr_lat_odd + 0.5) # latitude index
- air_dlat_1 = 360 / 59.0
+ latitude_even = float((360.0 / 60.0) * (j % 60 + cpr_lat_even))
+ latitude_odd = float((360.0 / 59.0) * (j % 59 + cpr_lat_odd))
- if self.odd_raw_latitude != 0:
- latr = self.odd_raw_latitude
- elif self.even_raw_latitude != 0:
- latr = self.even_raw_latitude
- else:
- latr = self.cfg_latitude
+ # Southern hemisphere: latitudes are between 270 to 360 degrees, make it +/- 90 degrees
+ if latitude_even >= 270:
+ latitude_even -= 360
+ if latitude_odd >= 270:
+ latitude_odd -= 360
- if self.odd_raw_longitude != 0:
- longr = self.odd_raw_longitude
- elif self.even_raw_longitude != 0:
- longr = self.even_raw_longitude
- else:
- longr = self.cfg_longitude
+ if CPR_NL(latitude_even) != CPR_NL(latitude_odd):
+ return False # Different latitude zones, not possible to compute position now
- tmp1 = math.floor(latr / air_dlat_1)
- tmp2 = (int(latr) % int(air_dlat_1))
- j = int(tmp1 + math.trunc(0.5 + tmp2 / air_dlat_1 - self.odd_raw_latitude / 131072.0))
- rlat = air_dlat_1 * (j + self.odd_raw_latitude / 131072.0)
- if rlat >= 270:
- rlat -= 360
+ # We have Ok conditions, set the latitude and continue with longitude calculation
+ latitude = latitude_even if self.even_time >= self.odd_time else latitude_odd
- if rlat < -90 or rlat > 90:
- return
+ if self.even_time >= self.odd_time:
+ ni = max(CPR_NL(latitude_even), 1)
+ dlon = 360.0 / ni
+ m = math.floor(cpr_lon_even * (CPR_NL(latitude_even) - 1) - cpr_lon_odd * CPR_NL(latitude_even) + 0.5)
+ longitude = dlon * (m % ni + cpr_lon_even)
+ else:
+ ni = max(CPR_NL(latitude_odd) - 1, 1)
+ dlon = 360.0 / ni
+ m = math.floor(cpr_lon_even * (CPR_NL(latitude_odd) - 1) - cpr_lon_odd * CPR_NL(latitude_odd) + 0.5)
+ longitude = dlon * (m % ni + cpr_lon_odd)
- if abs(rlat - latr) > (air_dlat_1 / 2):
- return
+ if longitude >= 180:
+ longitude -= 360
- air_dlon = 360 / max(CPR_NL(rlat) - 1, 1)
- m = int(math.floor(longr / air_dlon)
- + math.trunc(0.5 + (int(longr) % int(air_dlon)) / air_dlon - self.odd_raw_longitude / 131072.0))
- rlon = air_dlon * (m + self.odd_raw_longitude / 131072.0)
- if rlon > 180:
- rlon -= 360
+ self.data['latitude'] = str(round(latitude, 3)) if latitude != 0.0 else ""
+ self.data['longitude'] = str(round(longitude, 3)) if longitude != 0.0 else ""
- self.data['latitude'] = str(round(rlat, 3)) if rlat != 0.0 else ""
- self.data['longitude'] = str(round(rlon, 3)) if rlon != 0.0 else ""
- self.even_pos = False
- self.odd_pos = False
- self.even_then_odd_order = False
+ # Reset all flags
+ self.odd_time = 0
+ self.even_time = 0
- return
+ self.odd_raw_latitude = 0
+ self.odd_raw_longitude = 0
+ self.even_raw_latitude = 0
+ self.even_raw_longitude = 0
+
+ return True
def parse(self, obj):
"""
@@ -501,7 +473,8 @@ def decode_ADSB_msg(self):
if self.TC_ID_CAT_D_1 <= self.type_code <= self.TC_ID_CAT_A_4:
self['call_sign'] = callsign(hex(self.msg)[2:-1])
- elif self.type_code == self.TC_AIRBORNE_VELOCITY_19:
+
+ if self.type_code == self.TC_AIRBORNE_VELOCITY_19:
if 1 <= sub_type <= 4:
self.vertical_rate = self._get_vertical_rate()
if 1 <= sub_type <= 2:
@@ -545,29 +518,9 @@ def decode_ADSB_msg(self):
if self._get_msg_byte(5) & 0x04:
self['heading'] = str(((((self._get_msg_byte(5) & 0x03) << 8) | self._get_msg_byte(6)) * 45) >> 7)
- elif self.TC_SURFACE_POS_5 <= self.type_code <= self.TC_AIRBORNE_POS_22:
- odd = True if (self._get_msg_byte(6) & 0x04) else False
-
- if odd:
- self.odd_pos = True
- if self.even_pos:
- self.even_then_odd_order = True
- else:
- self.even_pos = True
-
- lat = ((self._get_msg_byte(6) & 0x03) << 15) | (self._get_msg_byte(7) << 7) | (self._get_msg_byte(8) >> 1)
- lon = ((self._get_msg_byte(8) & 0x01) << 16) | (self._get_msg_byte(9) << 8) | (self._get_msg_byte(10))
- if odd:
- self.odd_raw_latitude = lat
- self.odd_raw_longitude = lon
- else:
- self.even_raw_latitude = lat
- self.even_raw_longitude = lon
-
- if self.TC_AIRBORNE_POS_9 <= self.type_code <= self.TC_AIRBORNE_POS_18:
- self['altitude'] = str(self._get_altitude())
- self.on_ground = False
- elif self.TC_AIRBORNE_POS_20 <= self.type_code <= self.TC_AIRBORNE_POS_22:
+ if self.TC_SURFACE_POS_5 <= self.type_code <= self.TC_AIRBORNE_POS_22:
+ if (self.TC_AIRBORNE_POS_9 <= self.type_code <= self.TC_AIRBORNE_POS_18) \
+ or (self.TC_AIRBORNE_POS_20 <= self.type_code <= self.TC_AIRBORNE_POS_22):
self['altitude'] = str(self._get_altitude())
self.on_ground = False
elif self.TC_SURFACE_POS_5 <= self.type_code <= self.TC_SURFACE_POS_8:
@@ -575,13 +528,33 @@ def decode_ADSB_msg(self):
self['heading'] = str(self._get_heading())
self.on_ground = True
- elif self.type_code == self.TC_RESERVED_TEST_23:
+ if (self.TC_AIRBORNE_POS_9 <= self.type_code <= self.TC_AIRBORNE_POS_18) \
+ or (self.TC_AIRBORNE_POS_20 <= self.type_code <= self.TC_AIRBORNE_POS_22):
+ odd = True if (self._get_msg_byte(6) & 0x04) != 0 else False
+
+ lat = ((self._get_msg_byte(6) & 0x03) << 15) | \
+ (self._get_msg_byte(7) << 7) | \
+ (self._get_msg_byte(8) >> 1)
+ lon = ((self._get_msg_byte(8) & 0x01) << 16) | \
+ (self._get_msg_byte(9) << 8) | \
+ (self._get_msg_byte(10))
+
+ if odd:
+ self.odd_time = time.time()
+ self.odd_raw_latitude = lat
+ self.odd_raw_longitude = lon
+ else:
+ self.even_time = time.time()
+ self.even_raw_latitude = lat
+ self.even_raw_longitude = lon
+
+ if self.type_code == self.TC_RESERVED_TEST_23:
if sub_type == 7:
id_13 = (((self._get_msg_byte(5) << 8) | self._get_msg_byte(6)) & 0xFFF1) >> 3
if id_13 != 0:
self['squawk'] = str("{:=04X}".format(parse_id13(id_13)))
- elif self.type_code == self.TC_EXT_SQ_AIRCRFT_STATUS_28:
+ if self.type_code == self.TC_EXT_SQ_AIRCRFT_STATUS_28:
if sub_type == 1:
id_13 = ((self._get_msg_byte(5) << 8) | self._get_msg_byte(6)) & 0x1FFF
if id_13 != 0:
diff --git a/test_CPR.py b/test_CPR.py
index 9ee976a..1ab9e16 100644
--- a/test_CPR.py
+++ b/test_CPR.py
@@ -1,93 +1,66 @@
import basic
import math
+import squitter
+import time
-def CPR_NL(lat):
- lat = abs(lat)
- if lat >= max(basic.ADSB.NL):
- return 1
- elif lat <= min(basic.ADSB.NL):
- return len(basic.ADSB.NL) + 1
- else:
- return basic.ADSB.NL.index(max(filter(lambda x: x < lat, basic.ADSB.NL))) + 1
-
-
-def decodeCPR(odd_msg, even_msg):
- # We assume that odd_msg is consistent with even_msg on 'on_ground'
- air_dlat = 90.0 if odd_msg.on_ground else 360.0
- air_dlat_0 = air_dlat / 60.0
- air_dlat_1 = air_dlat / 59.0
- lat0 = even_msg.raw_latitude
- lat1 = odd_msg.raw_latitude
- lon0 = even_msg.raw_longitude
- lon1 = odd_msg.raw_longitude
-
- j = int(math.floor(((59 * lat0 - 60 * lat1) / 131072.0) + 0.5))
-
- rlat0 = air_dlat_0 * (j % 60 + lat0 / 131072.0) # 10.2157745361328
- rlat1 = air_dlat_1 * (j % 59 + lat1 / 131072.0) # 10.2162144547802
-
- if odd_msg.on_ground:
- surface_rlat = odd_msg.latitude if odd_msg.latitude != 0 else odd_msg.cfg_latitude
- surface_rlon = odd_msg.longitude if odd_msg.longitude != 0 else odd_msg.cfg_longitude
- rlat0 += math.floor(surface_rlat / 90.0) * 90.0
- rlat1 += math.floor(surface_rlat / 90.0) * 90.0
- else: # Adjust if we are on southern hemisphere by substracting 360 from latitude
- if rlat0 >= 270:
- rlat0 -= 360
- if rlat1 >= 270:
- rlat1 -= 360
-
- if rlat0 < -90 or rlat0 > 90 or rlat1 < -90 or rlat1 > 90:
- return None
- if CPR_NL(rlat0) != CPR_NL(rlat1):
- return None # The NL latitude value for odd and even frames must be equal
-
- ni = max(CPR_NL(rlat1) - 1, 1) # 58
- m = int(math.floor((((lon0 * (CPR_NL(rlat1) - 1)) - (lon1 * CPR_NL(rlat1))) / 131072.0) + 0.5)) # -39
- if odd_msg.on_ground:
- longitude = 90.0
- else:
- longitude = (360.0 / ni) * ((m % ni) + lon1 / 131072.0)
- latitude = rlat1
-
- if odd_msg.on_ground:
- longitude += math.floor(surface_rlon / 90.0) * 90.0
- elif longitude > 180:
- longitude -= 360
-
- return {'latitude': latitude, 'longitude': longitude}
-
-
-class CPR(basic.ADSB):
- def __init__(self):
- basic.ADSB.__init__(self)
-
- self.raw_latitude = 0
- self.raw_longitude = 0
- self.latitude = 0
- self.longitude = 0
- self.on_ground = False
- self.odd = False
- self.cfg_latitude = 0
- self.cfg_longitude = 0
-
def main():
- CPR1 = CPR()
- CPR2 = CPR()
-
- CPR1.raw_latitude = 88385
- CPR1.raw_longitude = 125818
- CPR1.odd = True
-
- CPR2.raw_latitude = 92095
- CPR2.raw_longitude = 39846
- CPR2.odd = False
-
- res = decodeCPR(CPR1, CPR2)
+ Sq = squitter.Squitter()
+
+ # Examples from
+ # https://adsb-decode-guide.readthedocs.io/en/latest/content/airborne-position.html
+ # http://www.lll.lu/~edward/edward/adsb/DecodingADSBposition.html
+
+ # For (93000, 74158, 51372, 50194) expected result is latitude = 52.257, longitude: 3.919
+ Sq.even_raw_latitude = 93000 # 92095
+ Sq.odd_raw_latitude = 74158 # 88385
+ Sq.even_raw_longitude = 51372 # 39846
+ Sq.odd_raw_longitude = 50194 # 125818
+ Sq.odd_time = 1
+ Sq.even_time = 2 # Make sure that even msg is later than odd
+
+ if Sq.decodeCPR():
+ print "latitude: {} longitude: {}".format(Sq['latitude'], Sq['longitude'])
+ print "Expected latitude: 52.257 longitude: 3.919"
+
+ # For (92095, 88385, 39846, 125818) expected result is latitude = 10.216 longitude: 123.889
+ Sq.even_raw_latitude = 92095
+ Sq.odd_raw_latitude = 88385
+ Sq.even_raw_longitude = 39846
+ Sq.odd_raw_longitude = 125818
+ Sq.odd_time = 1
+ Sq.even_time = 2 # Make sure that even msg is later than odd
+
+ if Sq.decodeCPR():
+ print "Got latitude: {} longitude: {}".format(Sq['latitude'], Sq['longitude'])
+ print "Expected latitude: 10.216 longitude: 123.889"
+
+ m1 = [0, 0x8D40621D58C386435CC412692AD6] # Odd
+ m2 = [0, 0x8D40621D58C382D690C8AC2863A7] # Even
+
+ Sq1 = squitter.Squitter()
+ Sq1.parse(m1)
+ Sq1.decode()
+
+ Sq2 = squitter.Squitter()
+ Sq2.parse(m2)
+ Sq2.decode()
+
+ Sq1.update(Sq2) # Expected result is latitude = 52.257, longitude: 3.919
+ if Sq1.decodeCPR():
+ print "latitude: {} longitude: {}".format(Sq1['latitude'], Sq1['longitude'])
+ print "Expected latitude: 52.257 longitude: 3.919"
+
+ m3 = [0, 0x8D40621D58C382D690C8AC2863A7]
+ Sq3 = squitter.Squitter()
+ Sq3.cfg_latitude = 52.258
+ Sq3.cfg_longitude = 3.918
+ Sq3.parse(m3)
+ Sq3.decode()
+ if Sq3.decodeCPR_relative():
+ print "latitude: {} longitude: {}".format(Sq3['latitude'], Sq3['longitude'])
+ print "Expected latitude: 52.25720 longitude: 3.91937"
- if res is not None:
- print "latitude: {}, longitude: {}".format(res['latitude'], res['longitude'])
if __name__ == '__main__':
- main()
\ No newline at end of file
+ main()