The Engineer’s Field Guide to Geospatial Analytics, Route Optimization, and GIS

From coordinate systems to contraction hierarchies: the data structures, algorithms, and tooling that power modern location intelligence.

Geospatial engineering sits at an unusual intersection: computational geometry, graph theory, combinatorial optimization, distributed systems, and a 200-year-old body of geodesy that refuses to be abstracted away. Most software engineers first encounter it when someone asks “how far apart are these two points?” and discover that the honest answer starts with “well, it depends on your ellipsoid.”

This article is a deep, practical tour of the field: the data models, the coordinate math, spatial indexing, the routing algorithm hierarchy from Dijkstra to Customizable Route Planning, vehicle routing optimization, map matching, and the open-source stack that ties it all together. Code is in Python and SQL, with pointers to production-grade engines where a from-scratch implementation would be malpractice.

1. The Two Data Models: Vector and Raster

Everything in GIS reduces to one of two representations.

Vector data models the world as discrete geometries: points, linestrings, polygons, and their multi-part variants, as formalized in the OGC Simple Features specification (ISO 19125). A road is a LineString, a delivery stop is a Point, a service territory is a MultiPolygon with holes. Vector data carries attributes (speed limits, population, zoning class) and supports exact topological predicates.

Raster data models the world as a regular grid of cells, each holding one or more band values: elevation (DEMs), satellite imagery, land cover classification, population density surfaces like WorldPop. Rasters trade topological precision for uniform sampling and array-based math — slope, hillshade, NDVI, and cost surfaces are all raster-native operations.

The serialization formats you’ll actually touch:

Format Model Notes WKT / WKB Vector Text/binary geometry encoding; the lingua franca of PostGIS GeoJSON (RFC 7946) Vector Always WGS84; web-friendly, verbose, no CRS flexibility per spec Shapefile Vector Legacy, ubiquitous, terrible (2GB limit, 10-char field names, multi-file) GeoPackage Vector+Raster SQLite-based, the sane modern replacement for shapefiles GeoParquet Vector Columnar, cloud-native, the analytics workhorse (works with DuckDB, Spark) FlatGeobuf Vector Streamable, includes a packed Hilbert R-tree for HTTP range-request reads Cloud-Optimized GeoTIFF (COG) Raster Internal tiling + overviews; range-request friendly MBTiles / PMTiles Tiles Tile pyramids for map rendering; PMTiles is single-file, serverless

If you’re building an analytics pipeline in 2026 and you’re not defaulting to GeoParquet + COG + PMTiles, you should have a specific reason.

2. Coordinate Reference Systems: Where “It Depends” Begins

A coordinate pair is meaningless without a CRS. A CRS is (1) a datum — an ellipsoid model of the Earth plus its anchoring — and (2) optionally a projection mapping the ellipsoid to a plane.

The ones that matter daily:

EPSG:4326 (WGS84) — geographic coordinates (longitude, latitude in degrees) on the WGS84 ellipsoid. GPS-native. Not a planar system: Euclidean math on degrees is wrong, and the error is latitude-dependent.
EPSG:3857 (Web Mercator) — the projection behind every slippy map. Conformal-ish, preserves angles locally, catastrophically distorts area at high latitudes (Greenland ≠ Africa). Never measure distance or area in 3857.
UTM zones (EPSG:326xx / 327xx) — 6°-wide transverse Mercator strips. Meters as units, <0.1% scale distortion within a zone. The default answer to “I need to buffer by 500 meters.”
Local/national grids — State Plane (US), OSGB (UK), etc. Use them when your data provider does.

The axis-order trap

WGS84’s official axis order is latitude, longitude (EPSG registry), but GeoJSON, WKB, PostGIS, and virtually all software use longitude, latitude (x, y). PROJ ≥ 6 honors authority order unless told otherwise, which is the source of an entire genre of “my points are in the ocean off Somalia” bugs:

from pyproj import Transformer

# always_xy=True forces lon/lat (x/y) order regardless of EPSG authority order
transformer = Transformer.from_crs("EPSG:4326", "EPSG:32633", always_xy=True)
x, y = transformer.transform(13.4050, 52.5200)  # Berlin → UTM 33N

Distance: haversine vs. geodesic

Three tiers of correctness for point-to-point distance on the Earth:

Equirectangular approximation — treat lat/lon as planar with a cos(lat) correction on longitude. Fast, fine under ~10 km.
Haversine — great-circle distance on a sphere. Error up to ~0.5% because Earth is oblate (flattening f ≈ 1/298.257).
Geodesic (Karney’s algorithm) — exact on the ellipsoid to nanometer precision; supersedes Vincenty, which fails to converge for near-antipodal points. This is what pyproj.Geod and GeographicLib implement.

import numpy as np

def haversine(lon1, lat1, lon2, lat2, R=6371008.8):
    """Vectorized haversine; inputs in degrees, output meters."""
    lon1, lat1, lon2, lat2 = map(np.radians, (lon1, lat1, lon2, lat2))
    dlon, dlat = lon2 - lon1, lat2 - lat1
    a = np.sin(dlat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dlon/2)**2
    return 2 * R * np.arcsin(np.sqrt(a))

from pyproj import Geod
geod = Geod(ellps="WGS84")
_, _, dist_m = geod.inv(13.4050, 52.5200, 2.3522, 48.8566)  # Berlin→Paris, exact

Rule of thumb: haversine for candidate filtering and ML features; geodesic when the number goes on an invoice; project to UTM when you need buffers, areas, or many pairwise operations in a local extent.

3. Spatial Indexing: Making O(n²) Disappear

Naïve spatial queries — “which of these 10M points fall in this polygon?” — are quadratic disasters. Spatial indexes fix this, and they come in two philosophical families.

3.1 Data-driven trees: the R-tree family

An R-tree groups geometries into hierarchies of minimum bounding rectangles (MBRs). Queries descend only into nodes whose MBRs intersect the query window. Variants: R*-tree (better split heuristics, the PostGIS GiST default behavior), STR-packed R-tree (bulk-loaded, read-only, what Shapely’s STRtree builds), Hilbert R-tree (orders entries along a Hilbert curve before packing — this is what's embedded in FlatGeobuf).

The critical mental model: the index only filters by bounding box. Every spatial query is a two-phase filter → refine pipeline. ST_Intersects(a, b) in PostGIS expands to "MBRs overlap (index) AND exact geometry test (refine)." This is why a query on long diagonal linestrings can be index-friendly yet still slow — the refine phase eats you.

import geopandas as gpd
from shapely import STRtree

points = gpd.read_file("deliveries.geoparquet")
zones  = gpd.read_file("zones.geoparquet").to_crs(points.crs)

# GeoPandas sjoin builds an STRtree under the hood
joined = gpd.sjoin(points, zones, how="inner", predicate="within")

3.2 Space-driven grids: geohash, S2, H3

Instead of adapting to the data, discrete global grid systems (DGGS) partition the planet into hierarchical cells and turn geometry problems into string/integer problems — which makes them shardable, joinable, and GROUP BY-able in ordinary data warehouses.

Geohash — Z-order (Morton) curve interleaving of lat/lon bits into base-32 strings. Prefix = containment. Weaknesses: cells are rectangles that vary wildly in shape with latitude, and Z-order has bad locality discontinuities at cell boundaries (two adjacent points can have no common prefix).
Google S2 — projects the sphere onto a cube, subdivides each face as a quadtree, orders cells by Hilbert curve. 64-bit cell IDs, ~1cm resolution at level 30. Excellent for covering arbitrary regions with unions of cells (S2RegionCoverer) — this is how "give me everything inside this polygon" becomes a set of integer range scans in Bigtable-class stores.
Uber H3 — hexagonal hierarchical grid (icosahedral gnomonic projection, aperture-7). Hexagons have a property rectangles don’t: all six neighbors are equidistant, with no ambiguous corner adjacency. This makes H3 the right choice for flow modeling, k-ring smoothing, gradient analysis, and ML feature grids. The cost: hexagons don’t perfectly nest across resolutions (the aperture-7 subdivision is rotational), so hierarchical rollups are approximate.

import h3

cell = h3.latlng_to_cell(52.5200, 13.4050, res=9)   # ~174m edge hexagon
ring = h3.grid_disk(cell, k=2)                       # 19 cells: self + 2 rings

# Aggregating a pickup dataset into an H3 heat surface:
df["h3_9"] = df.apply(lambda r: h3.latlng_to_cell(r.lat, r.lon, 9), axis=1)
surface = df.groupby("h3_9").size().rename("pickups").reset_index()

Decision heuristic: R-trees inside a database, H3/S2 across a data warehouse or between systems. They’re complementary, not competing.

3.3 PostGIS in anger

PostGIS remains the most complete spatial engine ever shipped. The operational essentials:

CREATE TABLE stops (
  id      bigint PRIMARY KEY,
  geom    geometry(Point, 4326)
);
CREATE INDEX stops_geom_gist ON stops USING gist (geom);

-- K-nearest-neighbor via the index-assisted distance operator <->
-- (order by <-> uses the GiST index; ST_Distance in ORDER BY does not)
SELECT id
FROM stops
ORDER BY geom <-> ST_SetSRID(ST_MakePoint(13.405, 52.52), 4326)
LIMIT 5;

-- Filter/refine join with a geography cast for meter-true distance
SELECT s.id, z.name
FROM stops s
JOIN zones z
  ON ST_Intersects(z.geom, s.geom);   -- && bbox filter is implicit

-- Fix the classic perf bug: ST_DWithin, not ST_Distance < x
SELECT count(*)
FROM stops
WHERE ST_DWithin(geom::geography,
                 ST_MakePoint(13.405, 52.52)::geography,
                 500);  -- meters, index-accelerated

Three PostGIS habits that separate professionals from tourists: ST_DWithin instead of ST_Distance comparisons (the former is sargable), ST_Subdivide on huge polygons before joins (caps refine-phase vertex counts), and geography type only at the edges of the system (its function coverage is narrower and math is slower than projected geometry).

The topological predicates themselves (ST_Intersects, ST_Contains, ST_Touches, ST_Crosses, …) are all projections of the DE-9IM (Dimensionally Extended 9-Intersection Model) matrix — a 3×3 matrix recording the dimensionality of intersections between the interiors, boundaries, and exteriors of two geometries. When the named predicates aren't precise enough, ST_Relate(a, b, 'T*T***T**')-style patterns give you the full algebra.

4. From Map to Graph: Building Routable Networks

Routing requires a graph, and real road networks are messier than textbook graphs in specific, important ways.

The source. OpenStreetMap is the default substrate. Ways carry tags (highway=motorway, maxspeed=100, oneway=yes, access=private); nodes shared between ways define topology. Commercial alternatives (HERE, TomTom) add lane-level detail, truck restrictions, and traffic history.

The modeling decisions that actually matter:

Node-based vs. edge-based graphs. The intuitive model — intersections as vertices, road segments as edges — cannot represent turn restrictions (“no left turn 7am–7pm”) or turn costs. The standard fix is the edge-based (line) graph: each directed road segment becomes a vertex, and each permitted turn becomes an edge. Vertex count roughly doubles-to-triples, but turn penalties become first-class edge weights. OSRM builds an edge-based graph; Valhalla keeps a node-based graph with explicit turn-cost lookup.
Weights are a function, not a number. “Shortest” (meters), “fastest” (seconds under a speed model), “cheapest” (tolls + fuel), truck-legal (height/weight constraints prune edges entirely), or time-dependent (weight = f(edge, departure_time), which breaks several speedup techniques — more below).
Strongly connected components. Raw OSM extracts contain thousands of tiny disconnected fragments (parking lots, private roads, mapping errors). Route between components → no path. Every serious pipeline keeps only the largest SCC (Tarjan’s or Kosaraju’s algorithm) per travel mode.

OSMnx makes the whole pipeline a few lines and is the right tool for research and city-scale analytics:

import osmnx as ox
import networkx as nx

G = ox.graph_from_place("Berlin, Germany", network_type="drive")
G = ox.add_edge_speeds(G)          # impute speeds from highway tags
G = ox.add_edge_travel_times(G)    # weight = length / speed

# largest strongly connected component
G = ox.truncate.largest_component(G, strongly=True)

orig = ox.nearest_nodes(G, X=13.3777, Y=52.5163)  # Brandenburg Gate
dest = ox.nearest_nodes(G, X=13.4132, Y=52.5219)  # Alexanderplatz
route = nx.shortest_path(G, orig, dest, weight="travel_time")

Note nearest_nodes — snapping is its own subproblem. Production systems snap to the nearest point on an edge (not just the nearest vertex), check side-of-street and access restrictions, and may return multiple candidates. Bad snapping is the #1 source of absurd routes.

5. The Shortest-Path Algorithm Hierarchy

There’s a clean intellectual ladder here, each rung trading preprocessing time and memory for query speed. On a continental graph (~50M edges), the spread between rungs is five orders of magnitude.

5.1 Dijkstra — the baseline

Settles vertices in increasing distance order using a priority queue. With a binary heap: O((V+E) log V). Correct for non-negative weights, no preprocessing, explores a roughly circular ball around the source — which is exactly the problem. A 700 km query settles millions of vertices, ~1s on a continental graph.

Practical notes that textbooks skip: use lazy deletion (push duplicates, skip stale pops) rather than decrease-key — it’s faster in practice despite worse asymptotics; a 4-ary heap beats binary on real hardware (cache lines); and early termination on settling the target is mandatory for point-to-point queries.

import heapq

def dijkstra(adj, source, target):
    """adj: dict[node] -> list[(neighbor, weight)]"""
    dist = {source: 0.0}
    parent = {source: None}
    pq = [(0.0, source)]
    while pq:
        d, u = heapq.heappop(pq)
        if u == target:
            break
        if d > dist.get(u, float("inf")):
            continue                      # stale entry (lazy deletion)
        for v, w in adj[u]:
            nd = d + w
            if nd < dist.get(v, float("inf")):
                dist[v], parent[v] = nd, u
                heapq.heappush(pq, (nd, v))
    path, u = [], target
    while u is not None:
        path.append(u); u = parent[u]
    return path[::-1], dist.get(target, float("inf"))

5.2 A* — goal direction

A* reorders Dijkstra’s exploration with a priority of g(v) + h(v), where h is a heuristic lower bound on remaining distance. If h is admissible (never overestimates) and consistent, A* is optimal and settles a subset of Dijkstra's vertices — the search ball becomes an ellipse pointed at the goal.

The catch on road networks: the obvious heuristic (haversine distance / max speed) is weak, because roads wiggle and speeds vary; typical speedup over Dijkstra is only 2–5×. A* earns its keep on grids and game maps; on road networks it’s a stepping stone.

5.3 Bidirectional search

Run forward search from s and backward search from t (on the reversed graph) simultaneously; stop when the frontiers meet (careful: the stopping condition is top(fwd) + top(bwd) ≥ best-seen-path, not "first common vertex"). Two balls of radius d/2 have half the total area of one ball of radius d — a reliable ~2× speedup, and the structural prerequisite for everything below.

5.4 ALT — A* with Landmarks and the Triangle inequality

Pick ~16–64 landmark vertices (farthest-point sampling works well), precompute exact distances from every vertex to every landmark. The triangle inequality then yields a much tighter admissible heuristic:

h(v) = max over landmarks L of | dist(L, t) − dist(L, v) |

Speedups of ~20–50×, preprocessing is cheap, and — crucially — ALT tolerates changing edge weights as long as they only increase (traffic!), because the landmark bounds remain admissible. This is why ALT variants survive in time-dependent settings where fancier techniques choke.

5.5 Contraction Hierarchies — the workhorse

CH (Geisberger et al., 2008) is the technique that made “instant” routing normal.

Preprocessing: impose a total order on vertices (by “importance” — a heuristic mix of edge difference, deleted neighbors, etc.). Contract vertices one by one from least to most important: when removing vertex v, for every pair of neighbors (u, w) whose shortest u→w path runs through v, insert a shortcut edge u→w with the combined weight. (A witness search — a small local Dijkstra — checks whether a path avoiding v is just as good, in which case no shortcut is needed.)

Query: bidirectional Dijkstra where the forward search only relaxes edges to more important vertices, and the backward search likewise. Both searches climb the hierarchy; they meet at some high-importance vertex. Search spaces shrink to a few hundred vertices; queries run in ~100 microseconds on continental graphs. Shortcuts are recursively unpacked to recover the real path.

The intuition for why this works is beautiful: road networks have low highway dimension — long shortest paths funnel through a tiny set of important vertices (motorways, bridges). CH is that structure made explicit.

Trade-offs: preprocessing takes minutes-to-hours and is metric-specific — change the cost function (or live traffic) and you must re-contract. Which leads to:

5.6 CRP / MLD — customizable route planning

Multi-Level Dijkstra / CRP (Delling et al., the algorithm behind OSRM’s MLD mode and originally developed for Bing Maps) splits preprocessing in two:

Metric-independent phase: partition the graph into nested cells (PUNCH, KaHIP, or Inertial Flow partitioners), computing a multilevel overlay topology. Slow, done rarely.
Customization phase: given a cost function, compute shortcut weights inside each cell — seconds to minutes, parallelizable.

Queries (~1ms) are slower than CH, but you can swap in fresh traffic every few minutes. The industry pattern: CH when the metric is static, CRP/MLD when it isn’t.

5.7 Many-to-many and the distance-matrix problem

VRP solvers (next section) need full N×N cost matrices. Running N² point-to-point queries is wasteful; the CH-based many-to-many algorithm (Knopp et al.) computes forward search spaces once per source and backward spaces once per target, then joins them in buckets — this is what OSRM’s /table endpoint and Valhalla's /sources_to_targets do. A 1,000×1,000 matrix takes ~1s. For truly huge matrices (100k+ points), switch to bucket-based batched CH or approximate with H3-cell centroid caching.

One more special case worth knowing: isochrones (reachability polygons — “everywhere reachable within 15 minutes”). Algorithmically it’s just Dijkstra with a cost cutoff; the hard part is turning the settled-node cloud into a polygon (alpha shapes / concave hulls, or Valhalla’s grid-marching approach). Isochrones power site selection, catchment analysis, and the entire “15-minute city” discourse.

6. Vehicle Routing: Where NP-Hard Meets the Loading Dock

Point-to-point routing is polynomial. The moment you ask “in what order should one or more vehicles visit these stops?” you’re in NP-hard territory.

6.1 The problem zoo

TSP — one vehicle, visit all stops, minimize tour cost. Solvable exactly to ~85,000 cities (Concorde, branch-and-cut over the subtour-elimination LP), but exact solving is rarely what industry needs.
CVRP — capacitated VRP: a fleet with capacity limits, minimize total distance/vehicles.
VRPTW — plus time windows per stop (“deliver between 9:00–11:00”). Feasibility alone becomes NP-complete.
PDPTW — pickup-and-delivery pairing (ride-sharing, courier).
The real-world monster — heterogeneous fleets, driver shifts and breaks (EU tachograph rules), multiple depots, split deliveries, time-dependent travel times, stochastic service durations, and re-optimization as new orders stream in. Every constraint you add interacts with every other one.

6.2 The algorithmic toolbox

Construction heuristics build a feasible solution fast. The classic is Clarke–Wright savings (1964): start with one out-and-back route per customer, then repeatedly merge the route pair with the largest saving s(i,j) = d(depot,i) + d(depot,j) − d(i,j), subject to feasibility. Still a respectable baseline sixty years later. Sweep and cheapest-insertion are the other standards.

Local search / improvement operators define a neighborhood over solutions:

2-opt / 3-opt — remove 2 (3) edges, reconnect. 2-opt uncrosses crossing edges.
Or-opt — relocate chains of 1–3 consecutive stops.
Lin–Kernighan / LKH — adaptive k-opt; LKH-3 is the strongest open TSP/VRP-variant heuristic in existence and routinely finds optima on huge instances.
Inter-route moves — relocate, swap, CROSS-exchange, 2-opt* between routes.

Metaheuristics escape local optima. The dominant one in modern VRP practice is Adaptive Large Neighborhood Search (ALNS): repeatedly destroy part of the solution (remove 10–40% of stops via random / worst-cost / Shaw-relatedness removal) and repair it (greedy or regret-k insertion), with operator selection probabilities adapted by historical success, wrapped in a simulated-annealing acceptance criterion. Nearly every commercial VRP engine is some flavor of ALNS. Hybrid Genetic Search (HGS) — population-based with local-search education — currently tops academic CVRP benchmarks (see Vidal’s HGS-CVRP).

Exact methods — branch-cut-and-price over set-partitioning formulations, with the pricing subproblem solved as an Elementary Shortest Path Problem with Resource Constraints. State of the art (VRPSolver) proves optimality to a few hundred customers. Use it for benchmarking heuristics, not dispatch.

6.3 OR-Tools in practice

Google OR-Tools’ RoutingModel is the pragmatic default: a constraint-programming core with ALNS-style large-neighborhood search on top, and a dimension abstraction that expresses capacities, time windows, and breaks uniformly.

from ortools.constraint_solver import pywrapcp, routing_enums_pb2

def solve_vrptw(time_matrix, time_windows, num_vehicles, depot=0,
                service_time=300, horizon=8*3600):
    n = len(time_matrix)
    manager = pywrapcp.RoutingIndexManager(n, num_vehicles, depot)
    routing = pywrapcp.RoutingModel(manager)

    def transit_cb(from_idx, to_idx):
        i = manager.IndexToNode(from_idx)
        j = manager.IndexToNode(to_idx)
        svc = service_time if i != depot else 0
        return time_matrix[i][j] + svc

    cb = routing.RegisterTransitCallback(transit_cb)
    routing.SetArcCostEvaluatorOfAllVehicles(cb)

    # Time dimension: accumulates transit, allows waiting (slack) at stops
    routing.AddDimension(cb, slack_max=3600, capacity=horizon,
                         fix_start_cumul_to_zero=False, name="Time")
    time_dim = routing.GetDimensionOrDie("Time")
    for node, (early, late) in enumerate(time_windows):
        if node == depot:
            continue
        time_dim.CumulVar(manager.NodeToIndex(node)).SetRange(early, late)

    params = pywrapcp.DefaultRoutingSearchParameters()
    params.first_solution_strategy = (
        routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC)
    params.local_search_metaheuristic = (
        routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
    params.time_limit.FromSeconds(30)

    solution = routing.SolveWithParameters(params)
    routes = []
    for v in range(num_vehicles):
        idx, route = routing.Start(v), []
        while not routing.IsEnd(idx):
            route.append(manager.IndexToNode(idx))
            idx = solution.Value(routing.NextVar(idx))
        routes.append(route)
    return routes

Field notes: (1) feed it real drive-time matrices from OSRM/Valhalla, never haversine — urban haversine underestimates by 30–60% and the “optimal” plan collapses on contact with one-way streets; (2) integers only — scale seconds/meters, no floats; (3) use AddDisjunction(node, penalty) to make stops droppable at a cost, otherwise one infeasible time window makes the entire problem infeasible; (4) Guided Local Search with a 10–60s budget is the sweet spot for daily dispatch of a few hundred stops. For heavier or research use, look at VROOM (open-source, ties directly into OSRM/Valhalla, very fast), PyVRP (HGS-based, excellent quality), and Timefold (JVM, rich constraint modeling).

6.4 Dynamic and stochastic dispatch

Real operations re-plan continuously: orders arrive, drivers run late, traffic shifts. The standard architecture is rolling-horizon re-optimization — freeze the next k stops per vehicle (drivers hate churn), re-solve the remainder every few minutes from the current solution (warm start), and price insertion of new orders via cheapest-insertion deltas against every active route. Research frontiers here involve reinforcement-learning policies for anticipatory dispatch and neural construction heuristics (attention models à la Kool et al.), which are genuinely competitive on small instances but still lose to ALNS/HGS at industrial scale.

7. Map Matching: From GPS Noise to Path Truth

GPS traces are noisy (5–20 m urban error, worse in canyons) and sparse (1–60 s sampling). Snapping each point to the nearest road independently produces garbage — it zig-zags across parallel streets. The canonical solution is the Hidden Markov Model formulation of Newson & Krumm (2009), which is what OSRM’s /match, Valhalla's Meili, and the Fast Map Matching (FMM) library all implement:

Hidden states = candidate road segments near each GPS point (from a spatial index query).
Emission probability = Gaussian in the perpendicular distance from point to candidate segment (σ ≈ GPS error, ~5 m).
Transition probability = decays exponentially in |route distance between consecutive candidates − haversine distance between the points|. Consecutive points shouldn’t require driving much farther than they moved.
Decode with the Viterbi algorithm → the jointly most probable path.

# Sketch of the Viterbi core; candidates[t] = list of (segment, emission_logp)
# route_cost(a, b) = network distance between candidate projections (via CH!)
import math

BETA = 1 / 20.0   # transition scale, meters^-1

def viterbi(candidates, gps_gaps, route_cost):
    V = [{c: e for c, e in candidates[0]}]
    back = [{}]
    for t in range(1, len(candidates)):
        V.append({}); back.append({})
        for c, e in candidates[t]:
            best, argbest = -math.inf, None
            for p in V[t-1]:
                trans = -BETA * abs(route_cost(p, c) - gps_gaps[t-1])
                score = V[t-1][p] + trans
                if score > best:
                    best, argbest = score, p
            V[t][c] = best + e
            back[t][c] = argbest
    # backtrack from argmax of V[-1] ...

The expensive part is the inner route_cost — thousands of small shortest-path queries per trace — which is why production matchers sit on top of a CH-accelerated engine. Map matching is upstream of almost everything interesting in mobility analytics: speed profiles per road segment, ETA training data, mileage-accurate billing, and detecting where the map itself is wrong.

8. The Engine Layer: What to Deploy

Engine Core technique Sweet spot OSRM CH or MLD, C++ Fastest raw queries; huge matrices; car-centric Valhalla Tiled A*-family, runtime costing Multimodal, truck attrs, isochrones, matrix, map matching, dynamic costing without recontraction GraphHopper CH + landmark hybrid, JVM Flexible custom models (declarative costing), strong bike/foot pgRouting Dijkstra/A* in PostGIS Small networks, SQL-native what-if analysis; not for serving traffic Concorde / LKH-3 Exact / k-opt TSP research, benchmark optima VROOM / PyVRP / OR-Tools ALNS / HGS / CP+LNS The VRP layer on top of any of the above

Architecture pattern that covers 90% of logistics products: Valhalla or OSRM for matrices and legs → PyVRP/VROOM/OR-Tools for sequencing → OSRM /match for after-the-fact telemetry → PostGIS + H3 for the analytics warehouse.

9. Scaling the Analytics Side

When the data outgrows one Postgres box:

DuckDB spatial — shockingly capable single-node engine; reads GeoParquet natively, does spatial joins over hundreds of millions of rows on a laptop. Try this before distributed anything.
Apache Sedona (Spark/Flink) — distributed spatial SQL with built-in spatial partitioning (KDB-tree / quad-tree partitioners) and partition-aware joins. The trick to distributed spatial joins is co-partitioning both sides by the same grid so intersecting geometries land on the same executor.
H3 as the join key — the warehouse-native pattern: assign every event an H3 cell at ingestion, and spatial joins degenerate into string-equality joins that BigQuery/Snowflake/ClickHouse eat for breakfast. Resolution choice is a bias/variance dial: res 8 (~460 m edge) for city analytics, res 9–10 for micro-mobility.
Vector tiles (MVT) for serving: ST_AsMVT in PostGIS, or tippecanoe → PMTiles for static pyramids. Never ship raw GeoJSON of a million features to a browser; tile it and let deck.gl / MapLibre do the rendering.

10. A Compact Mental Map

If you retain nothing else:

Geometry is meaningless without a CRS; measure in geodesic or projected meters, never in degrees or Web Mercator.
Every spatial query is filter-then-refine; the index prunes boxes, the geometry test decides. Choose R-trees inside engines, H3/S2 across systems.
Routing speed is bought with preprocessing: Dijkstra → A* → ALT → CH → CRP is a ladder of that trade, and road networks’ low highway dimension is why it works.
VRP is heuristics territory: construction + local search + ALNS/HGS, fed by real travel-time matrices, with droppable-stop penalties for graceful infeasibility.
Map matching is an HMM, and it converts noisy pings into the ground truth that trains your ETAs.
The open stack is genuinely excellent — PostGIS, GEOS/Shapely, GeoPandas, OSMnx, OSRM, Valhalla, PyVRP, DuckDB, H3 — and composing it well beats building any of it yourself.

The field rewards a particular kind of engineer: one who is comfortable jumping from an ellipsoid parameter to a heap invariant to an LP relaxation before lunch. If that sounds like fun, the roads — all 64 million kilometers of them — are waiting to be graphed.

Further reading: Bast et al., “Route Planning in Transportation Networks” (2016) — the definitive survey of the speedup-technique hierarchy; Toth & Vigo (eds.), “Vehicle Routing: Problems, Methods, and Applications”; Newson & Krumm (2009) for map matching; the PostGIS workshop (postgis.net) and the H3 documentation for hands-on grounding.

Brian James Curry is a Kansas City–based data scientist and applied Machine Learning and AI researcher working at the intersection of economic modeling, causal inference, and machine learning — with a standing preference for models where the economics is load-bearing. He is the founder of Vector1 Research, an independent lab exploring marketing economics, cognitive AI architecture, and AI-driven enterprise transformation, with research contributions including Papilon, PyCausalSim, and Memory-Node Encapsulation. His twenty-plus years of enterprise work span Koch Industries, Tractor Supply, Vail Resorts, Hallmark, Garmin, and AT&T.