Road map route finder using Dijkstra's algorithm and Blender

dijkstra_roadmap.py

import bpy
import random
import math
import collections
import json
import queue

westernmost_longitude = -130

def distance(a, b):
    ax, ay = a
    bx, by = b
    return math.sqrt((ax - bx)**2 + (ay - by)**2)

indices = {}
reverse_indices = {}
coords = []
edges = []

# https://hub.arcgis.com/datasets/esri::usa-freeway-system/explore?layer=1
with open('USA_Freeway_System.json') as f:
    data = json.load(f)

for feature in data['features']:
    feature_lines = []
    if feature['geometry']['type'] == 'LineString':
        feature_lines = [feature['geometry']['coordinates']]
    elif feature['geometry']['type'] == 'MultiLineString':
        feature_lines = feature['geometry']['coordinates']

    for line in feature_lines:
       for [x, y] in line:
           # ignore alaska and hawaii
           if x > westernmost_longitude:
               indices[(x, y)] = len(coords)
               reverse_indices[len(coords)] = (x, y)
               coords.append((x, y))
       for a, b in zip(line, line[1:]):
           ax, ay = a
           bx, by = b
           if ax > westernmost_longitude and bx > westernmost_longitude:
               edges.append((indices[tuple(a)], indices[tuple(b)]))

print('num coords:', len(coords))

start = (-118.226138727415, 34.0285571022828)
end = (-67.7811927273625, 46.1339890992705)

neighbors = collections.defaultdict(set)
for a, b in edges:
    neighbors[reverse_indices[a]].add(reverse_indices[b])
    neighbors[reverse_indices[b]].add(reverse_indices[a])

# dijkstra's algorithm
dist = {}
prev = {}
remaining = queue.PriorityQueue()

for coord in coords:
    dist[coord] = math.inf
    prev[coord] = None

dist[start] = 0
remaining.put((dist[start], start))

while not remaining.empty():
    _, u = remaining.get(False)

    if u == end:
        break

    for v in neighbors[u]:
        alt = dist[u] + distance(u, v)
        if alt < dist[v]:
            dist[v] = alt
            prev[v] = u
            remaining.put((alt, v))

s = []
u = end
while prev[u]:
    if u == start:
        break
    s = [u] + s
    u = prev[u]

# [coord] -> [edge]
print('path length:', len(s))
pathedges = []
for a, b in zip(s, s[1:]):
    pathedges.append((indices[a], indices[b]))

def generate_mesh(name, edges):
    mesh = bpy.data.meshes.new(name + "mesh")
    obj = bpy.data.objects.new(name, mesh)
    mesh.from_pydata([(x, y, 0) for (x, y) in coords], edges, [])
    mesh.update()
    bpy.context.collection.objects.link(obj)

generate_mesh('basemap', edges)
generate_mesh('path', pathedges)