added aprs plotting draft
[wikirepo] / notes / plotting_aprs_coverage / plot_aprs_log.py
diff --git a/notes/plotting_aprs_coverage/plot_aprs_log.py b/notes/plotting_aprs_coverage/plot_aprs_log.py
new file mode 100644 (file)
index 0000000..24d81ce
--- /dev/null
@@ -0,0 +1,209 @@
+#!/usr/bin/env python3
+
+import six
+import re
+import os
+from PIL import Image
+import cartopy.io.img_tiles
+from urllib.request import urlopen, Request
+import cartopy.crs as ccrs
+import matplotlib.pyplot as plt
+from collections import Counter
+import numpy as np
+
+## bodge our way around some changes in how cartopy interacts with OSM servers
+## see https://github.com/SciTools/cartopy/issues/1341
+def get_image_inner(self, tile):
+    url = self._image_url(tile) 
+    req = Request(url)
+    req.add_header('User-agent', 'your bot 0.1')
+    fh = urlopen(req)
+    im_data = six.BytesIO(fh.read())
+    fh.close()
+    img = Image.open(im_data)
+
+    img = img.convert(self.desired_tile_form)
+
+    return img, self.tileextent(tile), 'lower'
+
+cartopy.io.img_tiles.GoogleWTS.get_image = get_image_inner
+
+## Wrap a quick and dirty cache around the tile server to keep things
+## quick
+class OSMTilesCached(cartopy.io.img_tiles.OSM):
+    def __init__(self, *args, cachepath=None, **kwargs):
+        if cachepath is None:
+            self.cachepath = os.path.join('.', 'tile_cache')
+        else:
+            self.cachepath = cachepath
+
+        if not os.path.exists(self.cachepath):
+            os.mkdir(self.cachepath)
+
+        super().__init__(*args, **kwargs)
+
+    def get_image(self, tile):
+        cachefile = os.path.join(self.cachepath, str(tile) + '.png')
+
+        if not os.path.exists(cachefile):
+            (img, extent, origin) = super().get_image(tile)
+            img.save(cachefile)
+        else:
+            img = Image.open(cachefile)
+            origin = 'lower'
+            extent = self.tileextent(tile)
+        return img, extent, origin
+
+## interpret the simplist forms of lat and lon from aprs messages
+def lonlat_from_offset(msg, offset):
+    lonlat = (None, None)
+    if msg[offset].isnumeric():
+        signs = {'N': 1, 'E': 1, 'S':-1, 'W':-1}
+        lonlat = (msg[offset+9:offset+18], msg[offset:offset+8])
+        lonlat = ((int(lonlat[0][:3]) + 1/60*(float(lonlat[0][3:8])))*signs[lonlat[0][-1]],
+                  (int(lonlat[1][:2]) + 1/60*(float(lonlat[1][2:7])))*signs[lonlat[1][-1]])
+
+    elif msg[offset] == '/' or msg[offset] == 'S':
+        clon = [ord(x) - 33 for x in msg[offset+5:offset+9]]
+        clat = [ord(x) - 33 for x in msg[offset+1:offset+1+4]]
+        lon = -180 + (clon[0]*(91**3) + clon[1]*(91**2) + clon[2]*91 + clon[3])/190463
+        lat = 90 - (clat[0]*(91**3) + clat[1]*(91**2) + clat[2]*91 + clat[3])/380926
+        lonlat = (lon, lat)
+
+    else:
+        raise Exception('unknown lat lon format: +%d %s' % (offset, msg))
+
+    return lonlat
+            
+## parse a log generated by kissutil (from direwolf) to generate a list of
+## dicts represeting packets
+def load_packets_from_kisslog(logfilename):
+    packets = []
+
+    matchre = re.compile('^\[\d+ (\d{4}-\d{2}-\d{2} \d{2}:\d{2})\] ([^>]*)\>([^,]*),([^:]*):(.*)')
+    matchkeys = ['timestamp', 'call', 'type', 'path', 'msg']
+    packettypes = Counter()
+    with open(logfilename, 'r') as logfh:
+        for lineno, line in enumerate(logfh):
+            match = re.match(matchre, line)
+            if match:
+                packet = dict(zip(matchkeys, match.groups()))
+                packet['line'] = line
+                packet['path'] = [c.replace('*', '') for c in packet['path'].split(',') if not c.startswith('WIDE')]
+
+                packettypes[packet['msg'][0]] += 1
+
+                if packet['msg'].startswith('!') or \
+                   packet['msg'].startswith('='):
+                       #position without timestamp - no msg
+                       #position without timestamp - no aprs msg
+                       packet['lonlat'] = lonlat_from_offset(packet['msg'], 1)
+                elif packet['msg'].startswith('/') or \
+                   packet['msg'].startswith('@'):     
+                       #position with timestamp - no aprs msg
+                       #position with timestamp - aprs msg
+                       packet['lonlat'] = lonlat_from_offset(packet['msg'], 8)
+                elif packet['msg'].startswith('}'): #3rd party traffic
+                    pass
+                elif packet['msg'][0] in ['>', 'T', '<', ':', ';', '`', '\'']:
+                    pass
+                else:
+                    raise Exception('Unhandled packet type \'%s\'' % packet['msg'][0])
+
+
+                packets.append(packet)
+    return packets
+
+if __name__ == "__main__":
+    packets = load_packets_from_kisslog('aprs.log')
+    
+    print('%d packets loaded' % len(packets))
+    packets = [r for r in packets if 'lonlat' in r and r['lonlat'] != (None, None)]
+    print('%d packets with position loaded' % len(packets))
+
+    ## pre-populate with stations I know I don't have off-air locations for
+    stations = dict(M0ZRN=dict(call='M0ZRN', lonlat=(-0.01,52.14), src='file'),
+                             MB7UG=dict(call='MB7UG', lonlat=(0.969833,51.189000), src='file'),
+                             MB7UJ=dict(call='MB7UJ', lonlat=(-0.742000, 51.451500), src='file'),
+                             MB7UJD=dict(call='MB7UJD', lonlat=(1.375000, 52.815500), src='file'),
+                             MB7UNB=dict(call='MB7UNB', lonlat=(1.414833, 52.657667), src='file'),
+                             MB7UP=dict(call='MB7UP', lonlat=(-1.068000, 50.854333), src='file'),
+                             MB7USW=dict(call='MB7USW', lonlat=(-1.700333, 51.486833), src='file'),
+                             MB7UTG=dict(call='MB7UTG', lonlat=(0.130667, 52.736167), src='file'),
+                             MB7UW=dict(call='MB7UW', lonlat=( -1.357500, 51.065000), src='file'))
+                    
+    edges = Counter()
+    path_lengths = Counter()
+    ## process the packets to keep a list of stations heard, and a list edges between stations
+    for packet in packets:
+        if packet['call'] not in stations or stations[packet['call']]['src'] == 'file':
+            stations[packet['call']] = dict(call=packet['call'],
+                                             lonlat=packet['lonlat'],
+                                             src='radio',
+                                             count=1)
+        else:
+            stations[packet['call']]['count'] += 1
+
+        path = [packet['call']] + packet['path'] + ['M0ZRN']
+        path_lengths[len(path)] += 1
+        edges.update((path[i], path[i+1]) for i in range(0, len(path)-1))
+
+
+    ## now pop everything on a map
+    tiles = OSMTilesCached()
+    f = plt.figure()
+    f.set_size_inches(12, 12)
+    ax = plt.axes(projection=tiles.crs)
+    extent = [-2.5, 1.5, 50.5, 53]
+    ax.set_extent([-2.5, 1.5, 50.5, 53], crs=ccrs.PlateCarree())
+
+    def points_in_extent(points, extent):
+        def point_in_extent(point, extent):
+            return point[0] >= extent[0] and point[0] <= extent[1] \
+                    and point[1] >= extent[2] and point[1] <= extent[3]
+        return all(point_in_extent(point, extent) for point in points)
+
+    for start, stop in edges:
+        startp = stations[start]['lonlat']
+        stopp = stations[stop]['lonlat']
+
+        if not points_in_extent([startp, stopp], extent):
+            continue
+
+        ax.plot([startp[0], stopp[0]], [startp[1], stopp[1]],
+                transform=ccrs.PlateCarree(),
+                c='w',
+                lw=2*(0.5+np.log10(edges[(start, stop)])),
+                zorder=1)
+        ax.plot([startp[0], stopp[0]], [startp[1], stopp[1]],
+                transform=ccrs.PlateCarree(),
+                c='k',
+                lw=0.5+np.log10(edges[(start, stop)]),
+                zorder=2)
+
+    for station in stations:
+        color = {'file': 'r', 'radio': 'g'}[stations[station]['src']]
+
+        if not points_in_extent([stations[station]['lonlat']], extent):
+            continue
+
+        ax.scatter(stations[station]['lonlat'][0],
+                stations[station]['lonlat'][1],
+                c=color,
+                s=50,
+                transform=ccrs.PlateCarree(), zorder=3)
+        ax.annotate(station,
+                (stations[station]['lonlat'][0], stations[station]['lonlat'][1]),
+                bbox=dict(ec=color, fc='w', boxstyle='round,pad=0.1',),
+                xycoords=ccrs.PlateCarree()._as_mpl_transform(ax))
+
+    ax.add_image(tiles, 8)
+    ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
+
+    f.suptitle('APRS packets heard at M0ZRN between %s and %s' % (packets[0]['timestamp'], packets[-1]['timestamp']))
+    f.savefig('map.png', dpi=100)
+    plt.show()
+                    
+
+    
+