Merge branch 'cpr'
This commit is contained in:
commit
fdaa496b8f
194
python/cpr.py
194
python/cpr.py
@ -26,13 +26,9 @@ from air_modes.exceptions import *
|
|||||||
#the decoder is implemented as a class, cpr_decoder, which keeps state for local decoding.
|
#the decoder is implemented as a class, cpr_decoder, which keeps state for local decoding.
|
||||||
#the encoder is cpr_encode([lat, lon], type (even=0, odd=1), and surface (0 for surface, 1 for airborne))
|
#the encoder is cpr_encode([lat, lon], type (even=0, odd=1), and surface (0 for surface, 1 for airborne))
|
||||||
|
|
||||||
latz = 15
|
#TODO: remove range/bearing calc from CPR decoder class. you can do this outside of the decoder.
|
||||||
|
|
||||||
def nbits(surface):
|
latz = 15
|
||||||
if surface == 1:
|
|
||||||
return 19
|
|
||||||
else:
|
|
||||||
return 17
|
|
||||||
|
|
||||||
def nz(ctype):
|
def nz(ctype):
|
||||||
return 4 * latz - ctype
|
return 4 * latz - ctype
|
||||||
@ -49,38 +45,30 @@ def dlat(ctype, surface):
|
|||||||
else:
|
else:
|
||||||
return tmp / nzcalc
|
return tmp / nzcalc
|
||||||
|
|
||||||
def nl_eo(declat_in, ctype):
|
|
||||||
return nl(declat_in) - ctype
|
|
||||||
|
|
||||||
def nl(declat_in):
|
def nl(declat_in):
|
||||||
if abs(declat_in) >= 87.0:
|
if abs(declat_in) >= 87.0:
|
||||||
return 1.0
|
return 1.0
|
||||||
return math.floor( (2.0*math.pi) * pow(math.acos(1.0- (1.0-math.cos(math.pi/(2.0*latz))) / pow( math.cos( (math.pi/180.0)*abs(declat_in) ) ,2.0) ),-1.0))
|
return math.floor( (2.0*math.pi) * math.acos(1.0- (1.0-math.cos(math.pi/(2.0*latz))) / math.cos( (math.pi/180.0)*abs(declat_in) )**2 )**-1)
|
||||||
|
|
||||||
def dlon(declat_in, ctype, surface):
|
def dlon(declat_in, ctype, surface):
|
||||||
if surface == 1:
|
if surface:
|
||||||
tmp = 90.0
|
tmp = 90.0
|
||||||
else:
|
else:
|
||||||
tmp = 360.0
|
tmp = 360.0
|
||||||
nlcalc = nl_eo(declat_in, ctype)
|
nlcalc = max(nl(declat_in)-ctype, 1)
|
||||||
if nlcalc == 0:
|
return tmp / nlcalc
|
||||||
return tmp
|
|
||||||
else:
|
|
||||||
return tmp / max(nlcalc, 1.0)
|
|
||||||
|
|
||||||
def decode_lat(enclat, ctype, my_lat, surface):
|
def decode_lat(enclat, ctype, my_lat, surface):
|
||||||
tmp1 = dlat(ctype, surface)
|
tmp1 = dlat(ctype, surface)
|
||||||
tmp2 = float(enclat) / (2**nbits(surface))
|
tmp2 = float(enclat) / (2**17)
|
||||||
j = math.floor(my_lat/tmp1) + math.floor(0.5 + ((my_lat % tmp1) / tmp1) - tmp2)
|
j = math.floor(my_lat/tmp1) + math.floor(0.5 + ((my_lat % tmp1) / tmp1) - tmp2)
|
||||||
# print "dlat gives " + "%.6f " % tmp1 + "with j = " + "%.6f " % j + " and tmp2 = " + "%.6f" % tmp2 + " given enclat " + "%x" % enclat
|
|
||||||
|
|
||||||
return tmp1 * (j + tmp2)
|
return tmp1 * (j + tmp2)
|
||||||
|
|
||||||
def decode_lon(declat, enclon, ctype, my_lon, surface):
|
def decode_lon(declat, enclon, ctype, my_lon, surface):
|
||||||
tmp1 = dlon(declat, ctype, surface)
|
tmp1 = dlon(declat, ctype, surface)
|
||||||
tmp2 = float(enclon) / (2.0**nbits(surface))
|
tmp2 = float(enclon) / (2**17)
|
||||||
m = math.floor(my_lon / tmp1) + math.floor(0.5 + ((my_lon % tmp1) / tmp1) - tmp2)
|
m = math.floor(my_lon / tmp1) + math.floor(0.5 + ((my_lon % tmp1) / tmp1) - tmp2)
|
||||||
# print "dlon gives " + "%.6f " % tmp1 + "with m = " + "%.6f " % m + " and tmp2 = " + "%.6f" % tmp2 + " given enclon " + "%x" % enclon
|
|
||||||
|
|
||||||
return tmp1 * (m + tmp2)
|
return tmp1 * (m + tmp2)
|
||||||
|
|
||||||
@ -93,7 +81,11 @@ def cpr_resolve_local(my_location, encoded_location, ctype, surface):
|
|||||||
|
|
||||||
return [decoded_lat, decoded_lon]
|
return [decoded_lat, decoded_lon]
|
||||||
|
|
||||||
def cpr_resolve_global(evenpos, oddpos, mostrecent, surface):
|
def cpr_resolve_global(evenpos, oddpos, mypos, mostrecent, surface):
|
||||||
|
#cannot resolve surface positions unambiguously without knowing receiver position
|
||||||
|
if surface and mypos is None:
|
||||||
|
raise CPRNoPositionError
|
||||||
|
|
||||||
dlateven = dlat(0, surface)
|
dlateven = dlat(0, surface)
|
||||||
dlatodd = dlat(1, surface)
|
dlatodd = dlat(1, surface)
|
||||||
|
|
||||||
@ -114,7 +106,6 @@ def cpr_resolve_global(evenpos, oddpos, mostrecent, surface):
|
|||||||
#This checks to see if the latitudes of the reports straddle a transition boundary
|
#This checks to see if the latitudes of the reports straddle a transition boundary
|
||||||
#If so, you can't get a globally-resolvable location.
|
#If so, you can't get a globally-resolvable location.
|
||||||
if nl(rlateven) != nl(rlatodd):
|
if nl(rlateven) != nl(rlatodd):
|
||||||
#print "Boundary straddle!"
|
|
||||||
raise CPRBoundaryStraddleError
|
raise CPRBoundaryStraddleError
|
||||||
|
|
||||||
if mostrecent == 0:
|
if mostrecent == 0:
|
||||||
@ -122,21 +113,41 @@ def cpr_resolve_global(evenpos, oddpos, mostrecent, surface):
|
|||||||
else:
|
else:
|
||||||
rlat = rlatodd
|
rlat = rlatodd
|
||||||
|
|
||||||
dl = dlon(rlat, mostrecent, surface)
|
#disambiguate latitude
|
||||||
nlthing = nl(rlat)
|
if surface:
|
||||||
ni = max(nlthing - mostrecent, 1)
|
if mypos[0] < 0:
|
||||||
|
rlat -= 90
|
||||||
|
|
||||||
m = math.floor(((evenpos[1]*(nlthing-1)-oddpos[1]*(nlthing))/2**17)+0.5) #longitude index
|
dl = dlon(rlat, mostrecent, surface)
|
||||||
|
nl_rlat = nl(rlat)
|
||||||
|
|
||||||
|
m = math.floor(((evenpos[1]*(nl_rlat-1)-oddpos[1]*nl_rlat)/2**17)+0.5) #longitude index
|
||||||
|
|
||||||
|
#when surface positions straddle a disambiguation boundary (90 degrees),
|
||||||
|
#surface decoding will fail. this might never be a problem in real life, but it'll fail in the
|
||||||
|
#test case. the documentation doesn't mention it.
|
||||||
|
|
||||||
if mostrecent == 0:
|
if mostrecent == 0:
|
||||||
enclon = evenpos[1]
|
enclon = evenpos[1]
|
||||||
else:
|
else:
|
||||||
enclon = oddpos[1]
|
enclon = oddpos[1]
|
||||||
|
|
||||||
rlon = dl * (((ni+m) % ni)+enclon/2**17)
|
rlon = dl * ((m % max(nl_rlat-mostrecent,1)) + enclon/2.**17)
|
||||||
|
|
||||||
|
#print "DL: %f nl: %f m: %f rlon: %f" % (dl, nl_rlat, m, rlon)
|
||||||
|
#print "evenpos: %x, oddpos: %x, mostrecent: %i" % (evenpos[1], oddpos[1], mostrecent)
|
||||||
|
|
||||||
|
if surface:
|
||||||
|
#longitudes need to be resolved to the nearest 90 degree segment to the receiver.
|
||||||
|
wat = mypos[1]
|
||||||
|
if wat < 0:
|
||||||
|
wat += 360
|
||||||
|
zone = lambda lon: 90 * (int(lon) / 90)
|
||||||
|
rlon += (zone(wat) - zone(rlon))
|
||||||
|
|
||||||
|
#limit to (-180, 180)
|
||||||
if rlon > 180:
|
if rlon > 180:
|
||||||
rlon = rlon - 360.0
|
rlon -= 360.0
|
||||||
|
|
||||||
return [rlat, rlon]
|
return [rlat, rlon]
|
||||||
|
|
||||||
@ -172,52 +183,47 @@ def range_bearing(loc_a, loc_b):
|
|||||||
class cpr_decoder:
|
class cpr_decoder:
|
||||||
def __init__(self, my_location):
|
def __init__(self, my_location):
|
||||||
self.my_location = my_location
|
self.my_location = my_location
|
||||||
self.lkplist = {}
|
|
||||||
self.evenlist = {}
|
self.evenlist = {}
|
||||||
self.oddlist = {}
|
self.oddlist = {}
|
||||||
|
self.evenlist_sfc = {}
|
||||||
|
self.oddlist_sfc = {}
|
||||||
|
|
||||||
def set_location(new_location):
|
def set_location(new_location):
|
||||||
self.my_location = new_location
|
self.my_location = new_location
|
||||||
|
|
||||||
def weed_poslists(self):
|
def weed_poslists(self):
|
||||||
for poslist in [self.lkplist, self.evenlist, self.oddlist]:
|
for poslist in [self.evenlist, self.oddlist]:
|
||||||
for key, item in poslist.items():
|
for key, item in poslist.items():
|
||||||
if time.time() - item[2] > 900:
|
if time.time() - item[2] > 10:
|
||||||
|
del poslist[key]
|
||||||
|
for poslist in [self.evenlist_sfc, self.oddlist_sfc]:
|
||||||
|
for key, item in poslist.items():
|
||||||
|
if time.time() - item[2] > 25:
|
||||||
del poslist[key]
|
del poslist[key]
|
||||||
|
|
||||||
def decode(self, icao24, encoded_lat, encoded_lon, cpr_format, surface):
|
def decode(self, icao24, encoded_lat, encoded_lon, cpr_format, surface):
|
||||||
|
if surface:
|
||||||
|
oddlist = self.oddlist_sfc
|
||||||
|
evenlist = self.evenlist_sfc
|
||||||
|
else:
|
||||||
|
oddlist = self.oddlist
|
||||||
|
evenlist = self.evenlist
|
||||||
|
|
||||||
#add the info to the position reports list for global decoding
|
#add the info to the position reports list for global decoding
|
||||||
if cpr_format==1:
|
if cpr_format==1:
|
||||||
self.oddlist[icao24] = [encoded_lat, encoded_lon, time.time()]
|
oddlist[icao24] = [encoded_lat, encoded_lon, time.time()]
|
||||||
else:
|
else:
|
||||||
self.evenlist[icao24] = [encoded_lat, encoded_lon, time.time()]
|
evenlist[icao24] = [encoded_lat, encoded_lon, time.time()]
|
||||||
|
|
||||||
[decoded_lat, decoded_lon] = [None, None]
|
[decoded_lat, decoded_lon] = [None, None]
|
||||||
|
|
||||||
#okay, let's traverse the lists and weed out those entries that are older than 15 minutes, as they're unlikely to be useful.
|
#okay, let's traverse the lists and weed out those entries that are older than 10 seconds
|
||||||
self.weed_poslists()
|
self.weed_poslists()
|
||||||
|
|
||||||
if icao24 in self.lkplist:
|
if (icao24 in evenlist) \
|
||||||
#do emitter-centered local decoding
|
and (icao24 in oddlist):
|
||||||
[decoded_lat, decoded_lon] = cpr_resolve_local(self.lkplist[icao24][0:2], [encoded_lat, encoded_lon], cpr_format, surface)
|
newer = (oddlist[icao24][2] - evenlist[icao24][2]) > 0 #figure out which report is newer
|
||||||
self.lkplist[icao24] = [decoded_lat, decoded_lon, time.time()] #update the local position for next time
|
[decoded_lat, decoded_lon] = cpr_resolve_global(evenlist[icao24][0:2], oddlist[icao24][0:2], self.my_location, newer, surface) #do a global decode
|
||||||
|
|
||||||
elif (icao24 in self.evenlist) \
|
|
||||||
and (icao24 in self.oddlist) \
|
|
||||||
and (abs(self.evenlist[icao24][2] - self.oddlist[icao24][2]) < 10) \
|
|
||||||
and (surface == 0):
|
|
||||||
newer = (self.oddlist[icao24][2] - self.evenlist[icao24][2]) > 0 #figure out which report is newer
|
|
||||||
[decoded_lat, decoded_lon] = cpr_resolve_global(self.evenlist[icao24][0:2], self.oddlist[icao24][0:2], newer, surface) #do a global decode
|
|
||||||
self.lkplist[icao24] = [decoded_lat, decoded_lon, time.time()]
|
|
||||||
|
|
||||||
#so we really can't guarantee that local decoding will work unless you are POSITIVE that you can't hear more than 180nm out.
|
|
||||||
#this will USUALLY work, but you can't guarantee it!
|
|
||||||
#elif surface == 1 and self.my_location is not None:
|
|
||||||
# [local_lat, local_lon] = cpr_resolve_local(self.my_location, [encoded_lat, encoded_lon], cpr_format, surface) #try local decoding
|
|
||||||
# [rnge, bearing] = range_bearing(self.my_location, [local_lat, local_lon])
|
|
||||||
# if rnge < validrange: #if the local decoding can be guaranteed valid
|
|
||||||
# self.lkplist[icao24] = [local_lat, local_lon, time.time()] #update the local position for next time
|
|
||||||
# [decoded_lat, decoded_lon] = [local_lat, local_lon]
|
|
||||||
else:
|
else:
|
||||||
raise CPRNoPositionError
|
raise CPRNoPositionError
|
||||||
|
|
||||||
@ -232,83 +238,95 @@ class cpr_decoder:
|
|||||||
#encode CPR position
|
#encode CPR position
|
||||||
def cpr_encode(lat, lon, ctype, surface):
|
def cpr_encode(lat, lon, ctype, surface):
|
||||||
if surface is True:
|
if surface is True:
|
||||||
scalar = float(2**19)
|
scalar = 2.**19
|
||||||
else:
|
else:
|
||||||
scalar = float(2**17)
|
scalar = 2.**17
|
||||||
|
|
||||||
dlati = float(dlat(ctype, False))
|
#encode using 360 constant for segment size.
|
||||||
|
dlati = dlat(ctype, False)
|
||||||
yz = math.floor(scalar * ((lat % dlati)/dlati) + 0.5)
|
yz = math.floor(scalar * ((lat % dlati)/dlati) + 0.5)
|
||||||
rlat = dlati * ((yz / scalar) + math.floor(lat / dlati))
|
rlat = dlati * ((yz / scalar) + math.floor(lat / dlati))
|
||||||
|
|
||||||
nleo = nl_eo(rlat, ctype)
|
#encode using 360 constant for segment size.
|
||||||
if nleo == 0:
|
dloni = dlon(lat, ctype, False)
|
||||||
dloni = 360.0
|
|
||||||
else:
|
|
||||||
dloni = 360.0 / nl_eo(rlat, ctype)
|
|
||||||
|
|
||||||
xz = math.floor(scalar * ((lon % dloni)/dloni) + 0.5)
|
xz = math.floor(scalar * ((lon % dloni)/dloni) + 0.5)
|
||||||
|
|
||||||
yz = int(yz % scalar)
|
yz = int(yz) & (2**17-1)
|
||||||
xz = int(xz % scalar)
|
xz = int(xz) & (2**17-1)
|
||||||
|
|
||||||
return (yz, xz) #lat, lon
|
return (yz, xz) #lat, lon
|
||||||
|
|
||||||
if __name__ == '__main__':
|
if __name__ == '__main__':
|
||||||
import sys, random
|
import sys, random
|
||||||
|
|
||||||
rounds = 10000
|
rounds = 10001
|
||||||
threshold = 1e-3 #0.001 deg lat/lon
|
threshold = 1e-3 #0.001 deg lat/lon
|
||||||
#this accuracy is highly dependent on latitude, since at high
|
#this accuracy is highly dependent on latitude, since at high
|
||||||
#latitudes the corresponding error in longitude is greater
|
#latitudes the corresponding error in longitude is greater
|
||||||
|
|
||||||
bs = 0
|
bs = 0
|
||||||
|
surface = False
|
||||||
|
|
||||||
|
lats = [i/(rounds/170.)-85 for i in range(0,rounds)]
|
||||||
|
lons = [i/(rounds/360.)-180 for i in range(0,rounds)]
|
||||||
|
|
||||||
for i in range(0, rounds):
|
for i in range(0, rounds):
|
||||||
decoder = cpr_decoder(None)
|
even_lat = lats[i]
|
||||||
even_lat = random.uniform(-85, 85)
|
#even_lat = random.uniform(-85, 85)
|
||||||
even_lon = random.uniform(-180,180)
|
even_lon = lons[i]
|
||||||
odd_lat = even_lat + 1e-2
|
#even_lon = random.uniform(-180, 180)
|
||||||
odd_lon = min(even_lon + 1e-2, 180)
|
odd_lat = even_lat + 1e-3
|
||||||
|
odd_lon = min(even_lon + 1e-3, 180)
|
||||||
|
decoder = cpr_decoder([odd_lat, odd_lon])
|
||||||
|
|
||||||
#encode that position
|
#encode that position
|
||||||
(evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, False)
|
(evenenclat, evenenclon) = cpr_encode(even_lat, even_lon, False, surface)
|
||||||
(oddenclat, oddenclon) = cpr_encode(odd_lat, odd_lon, True, False)
|
(oddenclat, oddenclon) = cpr_encode(odd_lat, odd_lon, True, surface)
|
||||||
|
|
||||||
#perform a global decode
|
#try to perform a global decode -- this should fail since the decoder
|
||||||
|
#only has heard one position. need two for global decoding.
|
||||||
icao = random.randint(0, 0xffffff)
|
icao = random.randint(0, 0xffffff)
|
||||||
try:
|
try:
|
||||||
evenpos = decoder.decode(icao, evenenclat, evenenclon, False, False)
|
evenpos = decoder.decode(icao, evenenclat, evenenclon, False, surface)
|
||||||
#print "CPR global decode with only one report: %f %f" % (evenpos[0], evenpos[1])
|
|
||||||
raise Exception("CPR test failure: global decode with only one report")
|
raise Exception("CPR test failure: global decode with only one report")
|
||||||
except CPRNoPositionError:
|
except CPRNoPositionError:
|
||||||
pass
|
pass
|
||||||
|
|
||||||
|
#now try to do a real decode with the last packet's odd complement
|
||||||
|
#watch for a boundary straddle -- this isn't fatal, it just indicates
|
||||||
|
#that the even and odd reports lie on either side of a longitudinal boundary
|
||||||
|
#and so you can't get a position
|
||||||
try:
|
try:
|
||||||
(odddeclat, odddeclon, rng, brg) = decoder.decode(icao, oddenclat, oddenclon, True, False)
|
(odddeclat, odddeclon, rng, brg) = decoder.decode(icao, oddenclat, oddenclon, True, surface)
|
||||||
except CPRBoundaryStraddleError:
|
except CPRBoundaryStraddleError:
|
||||||
bs += 1
|
bs += 1
|
||||||
continue
|
continue
|
||||||
except CPRNoPositionError:
|
except CPRNoPositionError:
|
||||||
raise Exception("CPR test failure: no decode after even/odd inputs")
|
raise Exception("CPR test failure: no decode after even/odd inputs")
|
||||||
|
|
||||||
#print "Lat: %f Lon: %f" % (ac_lat, ac_lon)
|
|
||||||
|
|
||||||
if abs(odddeclat - odd_lat) > threshold or abs(odddeclon - odd_lon) > threshold:
|
if abs(odddeclat - odd_lat) > threshold or abs(odddeclon - odd_lon) > threshold:
|
||||||
print "odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat)
|
print "F odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat)
|
||||||
print "odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon)
|
print "F odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon)
|
||||||
raise Exception("CPR test failure: global decode error greater than threshold")
|
raise Exception("CPR test failure: global decode error greater than threshold")
|
||||||
|
# else:
|
||||||
|
# print "S odddeclat: %f odd_lat: %f" % (odddeclat, odd_lat)
|
||||||
|
# print "S odddeclon: %f odd_lon: %f" % (odddeclon, odd_lon)
|
||||||
|
|
||||||
nexteven_lat = odd_lat + 1e-2
|
nexteven_lat = odd_lat + 1e-3
|
||||||
nexteven_lon = min(odd_lon + 1e-2, 180)
|
nexteven_lon = min(odd_lon + 1e-3, 180)
|
||||||
|
|
||||||
(nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, False)
|
(nexteven_enclat, nexteven_enclon) = cpr_encode(nexteven_lat, nexteven_lon, False, surface)
|
||||||
|
|
||||||
|
#try a locally-referenced decode
|
||||||
try:
|
try:
|
||||||
(evendeclat, evendeclon) = cpr_resolve_local([even_lat, even_lon], [nexteven_enclat, nexteven_enclon], False, False)
|
(evendeclat, evendeclon) = cpr_resolve_local([even_lat, even_lon], [nexteven_enclat, nexteven_enclon], False, surface)
|
||||||
except CPRNoPositionError:
|
except CPRNoPositionError:
|
||||||
raise Exception("CPR test failure: local decode failure to resolve")
|
raise Exception("CPR test failure: local decode failure to resolve")
|
||||||
|
|
||||||
|
#check to see if the positions were valid
|
||||||
if abs(evendeclat - nexteven_lat) > threshold or abs(evendeclon - nexteven_lon) > threshold:
|
if abs(evendeclat - nexteven_lat) > threshold or abs(evendeclon - nexteven_lon) > threshold:
|
||||||
|
print "F evendeclat: %f nexteven_lat: %f evenlat: %f" % (evendeclat, nexteven_lat, even_lat)
|
||||||
|
print "F evendeclon: %f nexteven_lon: %f evenlon: %f" % (evendeclon, nexteven_lon, even_lon)
|
||||||
raise Exception("CPR test failure: local decode error greater than threshold")
|
raise Exception("CPR test failure: local decode error greater than threshold")
|
||||||
|
|
||||||
print "CPR test successful. There were %i boundary straddles over %i rounds." % (bs, rounds)
|
print "CPR test successful. There were %i boundary straddles over %i rounds." % (bs, rounds)
|
||||||
|
@ -69,7 +69,9 @@ class output_print(air_modes.parse):
|
|||||||
pass
|
pass
|
||||||
|
|
||||||
def output(self, msg):
|
def output(self, msg):
|
||||||
print self.parse(msg)
|
parsed = self.parse(msg)
|
||||||
|
if parsed is not None:
|
||||||
|
print self.parse(msg)
|
||||||
|
|
||||||
def print0(self, shortdata, ecc):
|
def print0(self, shortdata, ecc):
|
||||||
[vs, cc, sl, ri, altitude] = self.parse0(shortdata)
|
[vs, cc, sl, ri, altitude] = self.parse0(shortdata)
|
||||||
|
Loading…
Reference in New Issue
Block a user