commit 966ee378a20731d328472f5d1f32ecc748c9b7b8
Author: Pim Kunis <pim@kunis.nl>
Date:   Sun Aug 27 12:37:43 2023 +0200

    init

diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..bfb14b3
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,2 @@
+*.gpx
+*.pickle
diff --git a/create_polygon.py b/create_polygon.py
new file mode 100755
index 0000000..e2c14e4
--- /dev/null
+++ b/create_polygon.py
@@ -0,0 +1,64 @@
+#!/bin/python3
+
+import osmapi
+import shapely
+import pickle
+from os.path import exists
+
+
+OSM_KOGGENLAND_ID = 161930
+PICKLE_FILE = 'polygon.pickle'
+
+
+def chunks(lst, n):
+    """Yield successive n-sized chunks from lst."""
+    for i in range(0, len(lst), n):
+        yield lst[i:i + n]
+
+
+def dump_polygon():
+    osm_api = osmapi.OsmApi()
+    koggenland = osm_api.RelationGet(OSM_KOGGENLAND_ID)
+    ways = filter(lambda x: x['type'] == 'way', koggenland['member'])
+    shell_ids = []
+
+    for way in ways:
+        way = osm_api.WayGet(way['ref'])
+        node_ids = way['nd']
+        if not shell_ids:
+            shell_ids.append(node_ids[0])
+        for node_id in node_ids[1:]:
+            shell_ids.append(node_id)
+
+    node_data = dict()
+    for chunk in chunks(shell_ids, 100):
+        node_data |= osm_api.NodesGet(chunk)
+
+    shell = list()
+    for shell_id in shell_ids:
+        node = node_data[shell_id]
+        shell.append((node['lat'], node['lon']))
+
+    koggenland_polygon = shapely.Polygon(shell)
+    with open(PICKLE_FILE, 'wb') as f:
+        pickle.dump(koggenland_polygon, f)
+
+
+def main():
+    if not exists(PICKLE_FILE):
+        print('Using OSM API to create polygon...')
+        dump_polygon()
+        print('Dumped polygon to file.')
+    with open(PICKLE_FILE, 'rb') as f:
+        area = pickle.load(f)
+        print('Loaded polygon from file.')
+
+    point_out = shapely.Point(52.68811, 4.91267)
+    point_in = shapely.Point(52.69168, 4.91324)
+
+    print('out', area.contains(point_out))
+    print('in', area.contains(point_in))
+
+
+if __name__ == '__main__':
+    main()
diff --git a/find_centre.py b/find_centre.py
new file mode 100755
index 0000000..7c79b0d
--- /dev/null
+++ b/find_centre.py
@@ -0,0 +1,58 @@
+#!/bin/python3
+
+import osmapi
+import itertools
+import math
+# import pickle
+# import os
+
+OSM_KOGGENLAND_ID = 161930
+PICKLE_FILE = 'koggenland.pickle'
+
+
+def chunks(lst, n):
+    """Yield successive n-sized chunks from lst."""
+    for i in range(0, len(lst), n):
+        yield lst[i:i + n]
+
+
+def flatten(lst):
+    return list(itertools.chain(*lst))
+
+
+def main():
+    osm_api = osmapi.OsmApi()
+    koggenland = osm_api.RelationFullRecur(OSM_KOGGENLAND_ID)
+    node_ids = list(map(lambda a: a['data']['id'], filter(
+        lambda a: a['type'] == 'node', koggenland)))
+
+    coords = []
+    for chunk in chunks(node_ids, 100):
+        coords += map(lambda a: (a['lat'], a['lon']),
+                      osm_api.NodesGet(chunk).values())
+
+    cx, cy = 0., 0.
+    for (x, y) in coords:
+        cx += x
+        cy += y
+    cx /= len(coords)
+    cy /= len(coords)
+
+    print('centre', (cx, cy))
+
+    # 7KM radius (eye-balled)
+
+    # max_coord = None
+    # max_dist = 0
+
+    # for (x, y) in coords:
+    #     dist = math.sqrt((cx - x)**2 + (cy - y)**2)
+    #     if dist > max_dist:
+    #         max_dist = dist
+    #         max_coord = (x, y)
+
+    # print('max', max_coord)
+
+
+if __name__ == '__main__':
+    main()
diff --git a/parse_gpx.py b/parse_gpx.py
new file mode 100755
index 0000000..e4e6b93
--- /dev/null
+++ b/parse_gpx.py
@@ -0,0 +1,28 @@
+#!/bin/python3
+
+import gpxpy
+import gpxpy.gpx
+import shapely
+import pickle
+
+GPX_IN_FILE = 'caches.gpx'
+GPX_OUT_FILE = 'koggenland.gpx'
+PICKLE_FILE = 'polygon.pickle'
+
+
+def main():
+    with open(PICKLE_FILE, 'rb') as f:
+        area = pickle.load(f)
+
+    with open(GPX_IN_FILE, 'rb') as f:
+        gpx = gpxpy.parse(f)
+
+    gpx.waypoints = list(filter(lambda wp: area.contains(
+        shapely.Point(wp.latitude, wp.longitude)), gpx.waypoints))
+
+    with open(GPX_OUT_FILE, 'w') as f:
+        f.write(gpx.to_xml())
+
+
+if __name__ == '__main__':
+    main()