geo_json_to_gpx_json convert geojson to gpx

etl/etl_routemodel/etl_routemodel.py

@@ -8,6 +8,171 @@
 from geoalchemy2 import Geometry
 from geopy import distance
 from sqlalchemy import create_engine, inspect, String
+from collections import defaultdict
+from sklearn.neighbors import BallTree
+
+
+def geo_json_to_gpx_json(routes_filepath, routepoints_filepath, waypoints_filepath):
+    routeDataGJ = []
+    with open(routes_filepath, 'r') as f:
+        routeDataGJ = json.load(f)
+
+    routePointsGJ = []
+    with open(routepoints_filepath, 'r') as f:
+        routePointsGJ = json.load(f)
+        routePointsGJ = routePointsGJ["features"]
+
+    wayPointsGJ = []
+    with open(waypoints_filepath, 'r') as f:
+        wayPointsGJ = json.load(f)
+        wayPointsGJ = wayPointsGJ['features']
+
+    def flatten_route_points(rpGJ):
+        routePoints = []
+        for rp in rpGJ:
+            routePoints.append(
+                [
+                    rp["properties"]["route_fid"],
+                    rp["properties"]["route_point_id"],
+                    rp["geometry"]["coordinates"][0],
+                    rp["geometry"]["coordinates"][1],
+                ]
+            )
+        return np.array(routePoints)
+
+    allRPs = flatten_route_points(routePointsGJ)
+
+    # Could be useful but not being used
+    # def group_route_points(rpGJ):
+    #     groupedRPs = {}
+    #     currId = rpGJ[0]["properties"]["route_fid"]
+    #     points = []
+    #     for rp in rpGJ:
+    #         fid = rp["properties"]["route_fid"]
+    #         lng = rp["geometry"]["coordinates"][0]
+    #         lat = rp["geometry"]["coordinates"][1]
+    #         ele = rp["properties"]["ele"]
+
+    #         if fid != currId:
+    #             groupedRPs[currId] = points
+    #             currId = fid
+    #             points = []
+    #         points.append([lng, lat, ele])
+
+    #     return groupedRPs
+
+
+    # groupedRPs = group_route_points(routePointsGJ)
+
+    def assign_wps_to_routes(WPs, RPs):
+        coordinates = RPs[:, -2:]
+        routes = RPs[:, 0]
+        groupedWPs = defaultdict(list)
+
+        bt = BallTree(np.radians(coordinates), metric="haversine")
+        for wp in WPs:
+            lng = wp['geometry']['coordinates'][0]
+            lat = wp['geometry']['coordinates'][1]
+            # Increase k if a waypoint can be part of multiple routes, that is we can check the
+            # k closest points to see which routes they belong to (k = 1 => closest route only)
+            distances, indices = bt.query([np.radians([lng, lat])], k = 1)
+            closestRoutes = set()
+            for i in indices[0]:
+                closestRoutes.add(routes[i])
+            for r in closestRoutes:
+                groupedWPs[int(r)].append(wp)
+
+        return groupedWPs
+
+    groupedWPs = assign_wps_to_routes(wayPointsGJ, allRPs)
+
+    routesGJ = routeDataGJ['features']
+    routesJ = []
+    routeNames = []
+
+    # Iterate over each route/loop
+    for i, routeGJ in enumerate(routesGJ):
+        name = routeGJ['properties']['name']
+        routePointsGJ = routeGJ['geometry']['coordinates']
+        start = routePointsGJ[0]
+        end = routePointsGJ[-1]
+
+        points = []
+        prevP = None
+
+
+        # Iterate over each point in the current route
+        for p in routePointsGJ:
+
+            # Loops over all waypoints assigned to the current route and determines if it should be inserted
+            # at the current point
+            def findAndAddWP(prevP = None, p = None):
+                if not prevP:
+                    return False
+
+                # Determines if a wp is between p1 and p2 (wp is in the bounding box defined by p1 and p2)
+                def isInBoundingBox(wp, p1, p2):
+                    maxLng = max(p1[0], p2[0]) + 0.000045 # Might need to account for width of road in some cases?
+                    minLng = min(p1[0], p2[0]) - 0.000045
+                    maxLat = max(p1[1], p2[1]) + 0.000045
+                    minLat = min(p1[1], p2[1]) - 0.000045
+
+                    return wp[0] <= maxLng and wp[0] >= minLng and wp[1] <= maxLat and wp[1] >= minLat
+
+                for wp in groupedWPs[i]:
+                    lng = wp['geometry']['coordinates'][0]
+                    lat = wp['geometry']['coordinates'][1]
+                    if isInBoundingBox((lng, lat), prevP, p):
+                        points.append({
+                            'lng': wp['geometry']['coordinates'][0],
+                            'lat': wp['geometry']['coordinates'][1],
+                            "type": "waypoint"
+                        })
+                        points.append({
+                            'lng': p[0], 
+                            'lat': p[1]
+                        })
+                        groupedWPs[i].remove(wp)
+                        return True
+
+                return False
+
+            if not findAndAddWP(prevP, p):
+                points.append(
+                    {
+                        'lng': p[0], 
+                        'lat': p[1]
+                    }
+                )
+            prevP = p
+
+        # If groupedWPs[i] is not empty, potentialy waypoints at the start and end of loop
+        # They are not between existing route points so they would not be added in findAdnAddWP
+        for wp in groupedWPs[i]:
+            lng = wp['geometry']['coordinates'][0]
+            lat = wp['geometry']['coordinates'][1]
+            if distance.great_circle((lat, lng), (start[1], start[0])).m <= 50:
+                points.insert(0, {
+                    'lng': wp['geometry']['coordinates'][0],
+                    'lat': wp['geometry']['coordinates'][1],
+                    "type": "waypoint"
+                })
+            elif distance.great_circle((lat, lng), (end[1], end[0])).m <= 50:
+                points.append({
+                    'lng': wp['geometry']['coordinates'][0],
+                    'lat': wp['geometry']['coordinates'][1],
+                    "type": "waypoint"
+                })
+
+        route = {
+            "mode": "driving",
+            'points': points
+        }
+
+        routesJ.append(route)
+        routeNames.append(name)
+
+    return (routesJ, routeNames)
 
 
 def gpx_json_to_gdf(json_filepath):

etl/requirements.txt

@@ -28,6 +28,8 @@ shapely==2.0.2
 six==1.16.0
 soupsieve==2.5
 SQLAlchemy==2.0.22
+scikit-learn==1.4.1.post1
+scipy==1.12.0
 termcolor==2.3.0
 typing_extensions==4.8.0
 tzdata==2023.3