Spatial Epi Public Health GIS

Parallelizing Route Matrix Computation Safely

This guide, part of Batch Routing & Error Handling for Healthcare Access Analytics, covers running large origin–destination travel time matrices in parallel without losing, duplicating or corrupting results.

Problem Context & Constraints

A statewide access analysis might need travel times from 12,000 block groups to 2,400 clinics — nearly 29 million pairs. A self-hosted OSRM server answers a table request of a few hundred by a few hundred in well under a second, so the whole matrix is a matter of minutes if requests run in parallel, and hours if they run one at a time.

Parallelism brings its own failure modes. Overloading the router: too many concurrent requests exhaust its threads or memory, and it starts returning errors or timing out, which the client then retries, making things worse. Partial results: a worker that crashes halfway through leaves some chunks written and others not, and a naive restart either recomputes everything or silently assembles an incomplete matrix. Nondeterminism: results that arrive in completion order and are appended to one file end up in a different row order each run, so two runs cannot be compared by checksum. Shared state: several processes writing to the same output file or cache corrupt it.

The design principles that avoid these are simple: chunk the work into independent units with deterministic identifiers, bound concurrency to what the router can sustain, write each chunk atomically to its own file, and assemble and verify only when every chunk is present.

Chunk, Bound, Write Atomically, Verify Four stages. The origin–destination matrix is divided into a grid of chunks with deterministic identifiers derived from origin and destination block ranges. A bounded pool of workers, sized to the routing server's capacity, processes chunks. Each worker writes its chunk to a temporary file and renames it into place only when complete, so a crash never leaves a partial chunk. Assembly runs only when every chunk file exists, checks shapes and value ranges, and orders rows and columns deterministically before writing the final matrix with a checksum. Parallel, restartable, and the same answer every time Chunk origin × destination blocks, stable IDs Bounded pool workers ≤ router capacity Atomic write temp file, then rename into place Assemble all chunks present checks + checksum a restart skips chunks whose files already exist

Prerequisites

  • python 3.11
  • httpx 0.27.0, numpy 1.26.4, pandas 2.2.2
  • A self-hosted OSRM 5.27 server with the table service enabled and --max-table-size set above the chunk size
  • Origins and destinations with stable identifiers and coordinates in WGS84 longitude–latitude
  • Error handling for individual requests as in Handling API Timeouts in Batch OSM Routing

Step-by-Step Solution

Sort inputs by identifier, split into origin and destination blocks, run chunks through a thread pool bounded by a semaphore, write each chunk atomically, and assemble with verification.

# Parallel, restartable OSRM table computation with atomic chunk files.
# Pinned: httpx==0.27.0, numpy==1.26.4, pandas==2.2.2
import hashlib
import logging
import os
import time
from concurrent.futures import ThreadPoolExecutor, as_completed
from pathlib import Path

import httpx
import numpy as np
import pandas as pd

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("par_matrix")

OSRM = "http://osrm.internal:5000"
BLOCK = 250               # origins or destinations per chunk side
WORKERS = 8               # sized to the router's --threads, not the client's cores
RETRIES = 4


def chunks(origins: pd.DataFrame, dests: pd.DataFrame):
    o = origins.sort_values("id").reset_index(drop=True)
    d = dests.sort_values("id").reset_index(drop=True)
    for oi in range(0, len(o), BLOCK):
        for di in range(0, len(d), BLOCK):
            yield f"o{oi:07d}_d{di:07d}", o.iloc[oi:oi + BLOCK], d.iloc[di:di + BLOCK]


def fetch(client: httpx.Client, o: pd.DataFrame, d: pd.DataFrame) -> np.ndarray:
    coords = ";".join(f"{x:.6f},{y:.6f}" for x, y in pd.concat([o, d])[["lon", "lat"]].to_numpy())
    src = ";".join(map(str, range(len(o))))
    dst = ";".join(map(str, range(len(o), len(o) + len(d))))
    url = f"{OSRM}/table/v1/driving/{coords}?sources={src}&destinations={dst}&annotations=duration"
    for attempt in range(RETRIES):
        try:
            r = client.get(url, timeout=60)
            r.raise_for_status()
            body = r.json()
            if body.get("code") != "Ok":
                raise RuntimeError(body.get("code"))
            return np.array(body["durations"], dtype=float) / 60.0      # minutes; None -> nan
        except (httpx.HTTPError, RuntimeError) as e:
            wait = 2 ** attempt
            log.warning("retry %d after %s (%ss)", attempt + 1, e, wait)
            time.sleep(wait)
    raise RuntimeError("chunk failed after retries")


def run_chunk(client, out_dir: Path, cid: str, o, d) -> str:
    final = out_dir / f"{cid}.npz"
    if final.exists():
        return "skipped"
    m = fetch(client, o, d)
    tmp = out_dir / f".{cid}.{os.getpid()}.tmp.npz"
    np.savez(tmp, minutes=m, o=o["id"].to_numpy(), d=d["id"].to_numpy())
    os.replace(tmp, final)                       # atomic on the same filesystem
    return "done"


def compute(origins, dests, out_dir: Path) -> None:
    out_dir.mkdir(parents=True, exist_ok=True)
    work = list(chunks(origins, dests))
    limits = httpx.Limits(max_connections=WORKERS, max_keepalive_connections=WORKERS)
    with httpx.Client(limits=limits) as client, ThreadPoolExecutor(WORKERS) as pool:
        futs = {pool.submit(run_chunk, client, out_dir, cid, o, d): cid for cid, o, d in work}
        status = {"done": 0, "skipped": 0, "failed": 0}
        for f in as_completed(futs):
            try:
                status[f.result()] += 1
            except Exception as e:
                status["failed"] += 1
                log.error("chunk %s failed: %s", futs[f], e)
    log.info("chunks=%d %s", len(work), status)


def assemble(origins, dests, out_dir: Path) -> pd.DataFrame:
    expected = {cid for cid, _, _ in chunks(origins, dests)}
    present = {p.stem for p in out_dir.glob("o*_d*.npz")}
    missing = expected - present
    if missing:
        raise RuntimeError(f"{len(missing)} chunks missing; rerun compute before assembling")
    parts = []
    for cid in sorted(expected):
        z = np.load(out_dir / f"{cid}.npz", allow_pickle=True)
        parts.append(pd.DataFrame(z["minutes"], index=z["o"], columns=z["d"]).stack(dropna=False))
    M = pd.concat(parts).unstack().sort_index().sort_index(axis=1)
    if M.shape != (len(origins), len(dests)):
        raise RuntimeError(f"shape {M.shape} != expected")
    digest = hashlib.sha256(np.ascontiguousarray(M.to_numpy()).tobytes()).hexdigest()[:16]
    log.info("assembled %s, unreachable=%.2f%%, sha256_16=%s",
             M.shape, 100 * M.isna().to_numpy().mean(), digest)
    return M

The pool size is set by the router, not the client. OSRM processes each table request on one of its worker threads; sending more concurrent requests than it has threads only builds a queue on the server, raises latency and triggers client timeouts. Match WORKERS to the server’s --threads setting, then tune down if the server shares a machine with other work.

Because each chunk’s file name comes from its origin and destination offsets in sorted order, a rerun produces the same chunk boundaries and skips completed ones. The rename from a temporary file is what makes a chunk either fully present or absent: a crash mid-write leaves only a temporary file, which the next run ignores. The rename is only atomic within one filesystem, so keep temporary and final files in the same directory rather than writing temporaries to a system temp location on a different mount.

More Workers Is Not Always Faster A curve of pairs computed per second against the number of concurrent workers. Throughput rises almost linearly up to eight workers, which matches the routing server's thread count, then flattens, and beyond sixteen workers it falls as requests queue, time out and are retried. The recommended setting sits at the knee of the curve. Throughput peaks at the router's thread count 8 workers = router threads timeouts and retries pairs per second concurrent workers measure the curve once per server configuration

Validation & Edge Cases

Unreachable pairs. OSRM returns null for pairs it cannot route, such as islands without ferries or points snapped to disconnected road fragments. These become NaN, not zero. Report the unreachable share, and investigate origins with many unreachable destinations — usually a bad snap.

Snapping differences across chunks. OSRM snaps each coordinate to the nearest road. The same origin snaps identically in every chunk, so results are consistent, but a point snapped to a service road or a divided highway on the wrong side produces long detours. Snap origins once ahead of time with the nearest service, review the snap distances, and use the snapped coordinates in every chunk.

Changing data between runs. A restart after the router’s map data was updated mixes chunks from two network versions. Record the OSRM data timestamp in each chunk file and refuse to assemble chunks from different versions.

Memory at assembly. A 12,000 × 2,400 matrix of float64 values is about 230 MB, manageable, but the stacked intermediate in the code above is larger. For bigger matrices, write chunks into a preallocated memory-mapped array by offset rather than stacking.

Idempotence of reruns. Deleting one chunk file and rerunning must recompute only that chunk and produce the same assembled checksum. Make that a test: remove a random chunk, rerun, compare checksums.

Rate limits on hosted services. Against a commercial routing API, the bound is the provider’s rate limit rather than server threads. Implement a token bucket across workers rather than a pool size alone, since several workers can burst together and exceed a per-second limit.

Orphaned temporary files. A crash leaves temporary chunk files behind. They do no harm to correctness, since assembly reads only final names, but they accumulate. Clean files matching the temporary pattern that are older than the longest plausible chunk time at the start of each run.

INFO par_matrix chunks=480 {'done': 473, 'skipped': 7, 'failed': 0}
INFO par_matrix assembled (12004, 2398), unreachable=0.31%, sha256_16=5e1c0a9b44f2d7a3
What a Restart Sees A grid of chunk cells after a crash. Most cells hold completed chunk files and are skipped on restart. A few cells are empty because their workers had not started, and are computed. Two cells contain only temporary files from workers that crashed mid-write; those are ignored and their chunks recomputed. After the restart every cell holds a completed file and assembly proceeds. After a crash: skip, compute, ignore temporaries tmp tmp complete: skipped empty: computed temporary only: recomputed assembly waits until every cell is complete

Compliance Notes

Record the routing engine version and data timestamp, the profile, chunk size and worker count, retry policy, snap review results, the unreachable share, and the assembled matrix checksum. Keep chunk files until the matrix has been used and archived, so any cell can be traced to the request that produced it.

Origins are usually population-weighted block or block-group points, not addresses, and carry no personal data. If a matrix is ever computed from patient addresses, the request URLs contain those coordinates: send them only to a self-hosted router inside the secure environment, never to a public routing API, and exclude request logs from shared storage.