Modelling a Mobile Service Area: Why Postal Codes Lie and Drive-Time Polygons Do Not

Every business that sends a vehicle to a customer eventually writes the same function. It takes an address and returns a boolean. Somebody names it is_serviceable, somebody else names it in_zone, and the first implementation is almost always a Python set of postal code prefixes copied out of a spreadsheet that a dispatcher maintained by hand.

I write software for KMJ Tire, a small Calgary operation whose entire service catalogue is tires and oil changes. No mechanical work, no diagnostics, no suspension. That narrowness is a gift for a systems person, because it strips away most of the domain noise and leaves a handful of genuinely hard problems. The hardest one is not inventory and it is not scheduling. It is geography. A van leaves a yard, drives to wherever the vehicle is sitting, and does the work in a driveway or a parking lot. Whether we can profitably reach a given address, and how much of the day that address consumes, is the question the whole mobile service operation rests on.

This article is about getting that question right. It walks through why two common models — a postal code allowlist and a radius circle — fail in specific, instructive ways; what the correct geometric primitives are; how to pick a coordinate reference system for Alberta and why the default choice is wrong; what an isochrone actually is under the hood; how to test a point against a polygon at scale without lighting your CPU on fire; and how to store all of it so that a decision made in April can still be explained in September. Every number in the worked examples is illustrative and invented for the purpose of the example. None of them are operating statistics.

The Question Behind the Question

"Can we serve this address" is two questions wearing one coat.

The first is a feasibility question: is the location reachable by a service van within the operating constraints we have set? The second is an economics question: what does reaching it cost in van-minutes, and does that cost fit the job we are being asked to do? A twenty-minute tire rotation ninety minutes away is a loss even if the van can physically get there.

Naive systems answer only the first and pretend the second does not exist. That is how you end up with a dispatcher manually overriding the software every morning, which is the real symptom you are looking for when you audit one of these systems. If humans routinely override the boolean, the boolean is modelling the wrong thing.

The output we actually want is richer than a bit:

Field	Type	Meaning
covered	boolean	Reachable under the active policy
zone_version_id	uuid	Exactly which polygon produced the answer
drive_time_band	enum	t20, t35, t50, outside
confidence	enum	geocoded_exact, geocoded_interpolated, unresolved
decided_at	timestamptz	When the evaluation ran

Once the return type looks like that, most of the design follows. You need polygons with identity, you need versions, and you need a record of what was decided and why.

Postal Codes Are an Addressing Artifact, Not a Geography

A Canadian postal code has the form A1A 1A1. The first three characters form the Forward Sortation Area — the FSA. The last three form the Local Delivery Unit, the LDU. In T2P 1J9, T2P is downtown Calgary and 1J9 narrows it to a specific side of a specific block, or a specific large building, or a specific set of post office boxes.

The second character carries a signal that surprises people: a zero means rural. T0L is rural southwest Alberta. T2P is urban. That single digit is the difference between an FSA covering a few dozen city blocks and an FSA covering an area you could lose a small country in.

Here is the first structural problem. Postal codes are not areas. They are routes. Canada Post assigns them to sequences of delivery points so that a letter carrier or a rural driver can walk or drive them efficiently. The polygon you get when you download an "FSA boundary file" is a derived product — typically built by taking every dissemination block whose addresses carry that FSA and unioning them together. It is a reasonable approximation of where the addresses are. It is not an authoritative statement about where a line sits on the ground, and nobody at Canada Post is maintaining it as one.

The second problem is straddling. Municipal boundaries and postal routes were drawn by different organizations for different reasons and they were never reconciled. T3Z covers Springbank, which is in Rocky View County, west of the city limits. T1X spans Chestermere and reaches into county land around it. T3R includes both far northwest Calgary and acreage subdivisions well outside the city. If your allowlist contains T3Z, you have implicitly promised service to properties scattered across a corridor of foothills, some of which are a pleasant twenty-five minutes away and some of which are up a gravel range road that eats forty minutes in each direction.

The third problem is scale. A rural FSA such as T0M or T0J can cover hundreds of square kilometres. Including it in an allowlist is not a coverage decision; it is a coin flip with the geography.

The fourth problem is time. Calgary builds new communities continuously. When a new neighbourhood goes in, the addresses exist before the postal codes do, and the postal codes exist before anybody updates your spreadsheet. An allowlist written in 2023 is silently wrong about every subdivision occupied since. The failure is invisible: the request arrives, the lookup misses, and the customer is told no. Nothing logs an error, because from the software's perspective nothing went wrong.

The fifth problem is ownership. Ask who maintains the list. In every organization I have seen, the honest answer is "the person who made it, who has since moved to a different role." An unowned lookup table in a repository is technical debt that accrues interest in the form of lost work you never find out about.

A quick summary of what each failure mode costs you:

Failure mode	Symptom	How you find out
Route, not area	Boundary disputes at the edges	Dispatcher overrides
Straddling FSA	Wildly variable travel for one code	Van runs long
Rural FSA scale	Accidental hundred-kilometre promises	Angry driver
New development gap	Silent refusals	Never, usually
No owner	Slow drift from reality	Quarterly panic

The Radius Circle and the Lie of Euclidean Distance

The natural second attempt is a circle. Take the depot coordinates, pick a radius, and test whether the customer's point falls inside. It is one line of code and it feels principled because it is at least about distance.

It is wrong for a reason that has nothing to do with implementation quality: euclidean distance and network distance are different measurements, and the ratio between them is not a constant.

That ratio has a name in transport geography — the circuity ratio, or detour index. In a clean grid with no obstacles it hovers around 1.2 to 1.3. Rivers, rail corridors, and limited-access highways push it much higher in specific places. Calgary has all three. The Bow River cuts through the city from northwest to southeast and the Elbow joins it downtown. Two addresses six hundred metres apart across the Bow can be a seven-kilometre drive if the nearest crossing is inconvenient, and the number of crossings is small relative to the length of the river.

Limited-access roads invert the intuition in the other direction. Deerfoot Trail is a freeway with interchanges, not intersections. You cannot enter it from the adjacent residential street; you enter it where an interchange exists. So a point three hundred metres from a Deerfoot lane may be structurally far from Deerfoot. Meanwhile a location fifteen kilometres up the corridor with an interchange at its doorstep is functionally close.

Stoney Trail, the ring road, does the same thing at city scale. It makes perimeter movement dramatically cheaper than radial movement through the core. A circle drawn around the depot is isotropic — it assumes travel cost grows at the same rate in every direction. The actual cost surface is nothing like that. It is stretched along the freeway corridors and compressed across the river valleys and the rail yards.

There is a second, subtler defect. A circle drawn in degrees is not even a circle. If you write ST_DWithin(depot, customer, 0.35) against geometries in EPSG:4326, the unit is degrees, and a degree of longitude at Calgary's latitude is much shorter than a degree of latitude. Your "circle" is an ellipse squashed east-west by roughly 37 percent. I have found this bug in production more than once, and it never announces itself; it just quietly refuses service to the east and west while over-promising to the north and south.

Anisotropy, Deerfoot, and the 16:45 Problem

Consider two hypothetical service requests, both illustrative:

• Request A sits eight kilometres from the yard in a straight line, but on the far side of the core, and the van would have to cross the river and traverse downtown.
• Request B sits fifteen kilometres from the yard in a straight line, but almost all of that distance is ring-road running with an interchange at each end.

At 10:30 on a Tuesday these might both be roughly twenty-two minutes. At 16:45 on a Thursday in October, Request A is forty-plus minutes and Request B is still around twenty-five. A radius model ranks A as nearer than B in both scenarios and is therefore wrong exactly when being wrong is most expensive — during the afternoon peak of the autumn changeover, when the seasonal swap rush has every slot in the day already committed.

This is anisotropy: cost depends on direction, not just magnitude. And it is non-stationary: the anisotropy field itself changes through the day. Any model that collapses to a single scalar radius has thrown away both properties.

The right abstraction is a level set of the travel-time function. Fix the origin, fix the departure time, fix the vehicle profile, and ask: what is the set of all points reachable within t minutes? That set is a polygon. It is lumpy, it has fingers running along the freeways, it has bites taken out of it where the river blocks movement, and it changes shape depending on when you ask. That polygon is an isochrone, and it is the honest answer.

Geometry Primitives You Actually Need

Before writing anything spatial, fix the vocabulary. The OGC Simple Features model gives you a small set of types and you will use four of them.

A Point is a coordinate pair. In a geographic system it is (longitude, latitude) in that order — which is the reverse of how humans say it, and the source of a truly enormous number of bugs. GeoJSON, PostGIS, and most routing APIs use x-then-y. Leaflet's LatLng does not. Write a single conversion helper and never do it inline.

A LineString is an ordered sequence of points. Roads are linestrings. Your depot-to-customer path is a linestring. Nothing about it is closed.

A Polygon is one exterior ring plus zero or more interior rings, where a ring is a closed linestring whose first and last coordinates are identical. Interior rings are holes. A service zone with a genuine exclusion — a military reserve, a large restricted industrial site — is a polygon with a hole, not two polygons.

A MultiPolygon is a set of disjoint polygons treated as one geometry. This is what you actually want for a service zone, because real coverage is often disconnected. A town twenty minutes down the highway with nothing serviceable in between is a separate part of the same zone.

Winding order deserves its own paragraph because it will eventually cost you a day. RFC 7946, the GeoJSON specification, requires exterior rings to be counterclockwise and interior rings clockwise — the right-hand rule. The Esri shapefile format uses the opposite convention for exteriors. PostGIS mostly does not care, because it treats rings as planar and infers holes from containment. Spherical geometry libraries such as S2, and by extension several cloud data warehouses, care intensely: on a sphere a ring divides the surface into two finite regions, and orientation is the only thing that tells you which one you meant. Get it backwards and your Calgary service zone becomes "the entire planet except Calgary." The symptom is unmistakable once you have seen it, and baffling the first time.

Two invariants worth asserting on ingest:

• The geometry is valid under OGC rules. Self-intersecting exterior rings are the usual offender; ST_IsValid plus ST_IsValidReason will name the problem, and ST_MakeValid will usually fix it at the cost of silently changing your data.
• Ring orientation is normalized. ST_ForcePolygonCCW on the way in means downstream consumers never have to guess.

Choosing a Coordinate Reference System for Alberta

EPSG:4326 is WGS 84, a geographic coordinate system whose units are degrees of angle. It is the right storage and interchange format. It is the wrong computation format for anything involving length or area, and the reason is straightforward: a degree is not a distance.

At Calgary's latitude, roughly 51.05 degrees north, one degree of latitude is about 111.2 kilometres. One degree of longitude is that same figure scaled by the cosine of the latitude — about 0.629 — giving roughly 70.0 kilometres. So the coordinate grid is anisotropic by a factor of 1.59 before you have done anything at all. Any planar distance, area, or buffer computed on raw lat/lon inherits that distortion.

The three ways out, in ascending order of how much I like them:

Use geography types. PostGIS geography computes on the spheroid and returns metres. ST_DWithin(a::geography, b::geography, 25000) is correct and readable. The cost is that the geography operator set is narrower and the functions are slower, sometimes considerably, on complex polygons.

Use Web Mercator. Do not. EPSG:3857 is a presentation projection. Its scale factor is the secant of the latitude, so at 51 degrees north every length is inflated by about 1.59 and every area by about 2.53. It exists so that map tiles are square. It has no business in a measurement pipeline, and the fact that it is the default in several visualization stacks has produced a lot of quietly wrong analytics.

Use a local projected CRS. This is the right answer for repeated planar work. The question is which one.

The obvious candidate is UTM. Calgary sits at about 114.07 degrees west, which places it in UTM zone 11N, EPSG:32611 — but only just. The boundary between zone 11 and zone 12 runs along the 114th meridian, roughly seven kilometres east of downtown. So the moment your service area extends east toward Chestermere at 113.7 west or Strathmore at 113.4 west, you are working across a zone boundary. You can force everything into zone 11 and accept the extrapolated distortion; near the zone edge the scale factor drifts to roughly 1.001, which is a metre per kilometre. Over a forty-kilometre run that is forty metres of error. Tolerable for coverage decisions, irritating for anything else, and philosophically unsatisfying.

Alberta publishes a better option. The provincial 3TM system uses three-degree-wide transverse Mercator zones with central meridians at 111, 114, 117, and 120 degrees west, and — critically — a scale factor of exactly 1.0 at the central meridian rather than UTM's 0.9996. The 114 W zone spans 112.5 to 115.5 degrees west, which comfortably contains Cochrane, Airdrie, Calgary, Okotoks, Chestermere, and Strathmore in a single undistorted frame with the central meridian running essentially through the city.

That definition is commonly catalogued as EPSG:3401 for the NAD83 datum, with NAD83(CSRS) realizations in the neighbouring 3775–3779 block. Verify against your own PROJ installation rather than trusting a number from an article — projinfo EPSG:3401 will print the WKT and you can confirm the central meridian yourself. Getting the datum realization wrong is a sub-metre error and usually irrelevant for coverage; getting the central meridian wrong is not.

A comparison of the candidates on the axes that matter:

CRS	Units	Distortion near Calgary	Verdict
EPSG:4326	degrees	n/a — not planar	Storage and interchange
EPSG:3857	metres	~1.59x length, ~2.53x area	Tiles only
EPSG:32611	metres	~1.001 near zone edge	Workable, awkward east
Alberta 3TM 114 W	metres	~1.0 across the region	Preferred for planar math
geography	metres	spheroidal, correct	Good default when speed allows

My standing recommendation: store in 4326, index in 4326, and reproject to the 3TM 114 W frame inside any function that buffers, measures, or computes area. Keep the reprojection at one boundary in the code so it is auditable.

What an Isochrone Really Is

An isochrone is a contour of the travel-time field. The definition is clean; the construction is not, and understanding the construction tells you what to distrust.

A routing engine — Valhalla, OSRM with a custom profile, GraphHopper, pgRouting over an imported network — holds the road network as a directed graph. Nodes are intersections, edges are road segments, and each edge carries a traversal cost derived from length, classification, speed data, turn restrictions, and access rules for the vehicle profile.

To build a fifteen-minute isochrone the engine runs a single-source shortest-path expansion from the origin, typically Dijkstra with a cost ceiling instead of a target node. It stops expanding when accumulated cost exceeds the budget. What comes back is not a polygon; it is a labelled subgraph — every node with its minimum cost from the origin.

Turning that into an area requires interpolation. Valhalla's approach is representative: it rasterizes costs onto a regular grid, assigning each cell the lowest cost reachable within it, then runs marching squares over that grid to extract the iso-cost contour, then generalizes and polygonizes the result. Every step of that pipeline introduces artifacts:

• Grid resolution determines how faithfully narrow features survive. A single road threading through otherwise unreachable terrain may be widened into a lobe or dropped entirely.
• Generalization trades vertex count for fidelity. A polygon with two hundred vertices tests fast and lies a little. One with twenty thousand vertices tests slowly and lies less.
• Holes appear where an enclosed area is genuinely unreachable within the budget, and also where the grid was too coarse to notice the road that reaches it.

Time-of-day dependence is the property people forget. If the engine is running free-flow speeds derived from road classification, the polygon describes an empty city at three in the morning. Valhalla accepts a departure timestamp and will use predicted traffic if you have loaded a traffic tile set. OSRM has no native concept of time; the usual workaround is to compile separate graphs with per-bucket speed tables and route against whichever matches the hour. Either way, the polygon is a function of (origin, budget, profile, departure_bucket, graph_version), and if your cache key omits any of those, you will eventually serve an answer generated for a different world.

Versioning and Caching Drive-Time Polygons

Generating isochrones on the request path is a mistake. It is a network hop to a service that is doing a graph expansion and a contouring pass, and it puts a third-party dependency in the latency budget of something a customer is waiting on.

Generate them offline. A batch job runs per origin, per budget, per departure bucket, writes the resulting geometry into the database with full provenance, and marks it active. Coverage evaluation then reads a local polygon and never touches the routing engine.

The provenance record should be complete enough to regenerate the exact polygon:

{
  "origin": {"lon": -114.0719, "lat": 51.0447},
  "cost_seconds": 1800,
  "profile": "auto",
  "departure_bucket": "weekday_pm_peak",
  "engine": "valhalla",
  "engine_version": "3.4.0",
  "graph_extract": "alberta-2026-07-14.osm.pbf",
  "generalize_m": 40,
  "denoise": 0.2,
  "generated_at": "2026-07-15T04:12:08Z"
}

Three rules keep this honest.

First, versions are immutable. You never update a geometry in place. A new generation writes a new row with a new version number and a validity window; the old row stays exactly as it was. This is what makes historical explanation possible, and it is worth every byte.

Second, the cache key is the full provenance tuple, not the origin alone. I have watched a team spend a day debugging "the polygon changed by itself" when in fact two jobs with different departure buckets were writing to the same key.

Third, regeneration is scheduled against the graph extract, not the calendar. When the OpenStreetMap extract updates and the network topology changes — a new interchange opens, a bridge closes for construction — the polygons are stale regardless of how recently you rebuilt them. Tie the trigger to the input, not to a cron expression that somebody picked arbitrarily.

Point in Polygon: Ray Casting, Winding Numbers, and the Gap Between Them

Given a polygon and a point, decide containment. Two classical algorithms, both O(n) in the vertex count, both worth understanding before you delegate to a library.

Ray casting, also called the crossing number or even-odd rule, shoots a ray from the query point in a fixed direction — conventionally along positive x — and counts how many polygon edges it crosses. Odd means inside, even means outside. It is four lines of arithmetic per edge and it is what most naive implementations use.

Winding number computes how many times the polygon boundary wraps around the point, summing signed angles or, more practically, using orientation tests on each edge. Nonzero means inside.

For a simple, non-self-intersecting polygon the two agree everywhere. They diverge on self-intersecting geometry: a figure-eight's crossover lobe is outside under even-odd and inside under nonzero winding. Real isochrone output is not always clean, especially after generalization, so this is not purely academic. My preference is to validate and repair geometries on ingest so that the question never arises, and then use whichever predicate the library implements — but to know which one that is.

Complexity matters once you have more than a couple of zones. The naive loop is O(Z·n) for Z zones averaging n vertices. With a dozen active zone versions at a few thousand vertices each, a single evaluation touches tens of thousands of edges. At a few hundred requests a minute that is measurable, and it is entirely avoidable.

Spatial Indexes and the Bounding-Box Prefilter

The fix is the standard two-phase spatial query: cheap filter, then exact refinement.

Phase one uses an R-tree. Every geometry contributes its minimum bounding rectangle to a balanced tree of nested rectangles. A point query descends the tree and returns the small set of candidates whose MBR contains the point. That is O(log Z) expected, and it eliminates the overwhelming majority of zones with four float comparisons each.

Phase two runs the exact predicate on the survivors only. If your zones are roughly convex and well separated, phase one usually leaves one candidate. If they are long and diagonal — a corridor along a highway is the classic case — the MBR is a poor approximation and you will get more false positives, which is an argument for splitting a straggly zone into a MultiPolygon of tighter parts.

PostGIS gives you this for free. A GIST index over a geometry column is an R-tree over the bounding boxes, implemented on top of the generalized search tree framework. The && operator is bounding-box overlap and is index-accelerated. Functions like ST_Contains and ST_Intersects are defined as && composed with the exact underscore-prefixed predicate, so the planner uses the index automatically. What is not automatic: an index on geom does nothing for a query that wraps the column in ST_Transform. Either index the expression or store a second projected column.

Shapely's STRtree is a packed Sorted-Tile-Recursive R-tree — built once, queried many times, not modifiable afterwards. That fits the isochrone workload exactly, because zone versions change on a batch schedule rather than per request.

One more trick worth its weight: prepared geometries. Both GEOS and Shapely can precompute an edge index for a polygon that will be tested repeatedly. shapely.prepared.prep returns an object whose contains is dramatically faster on high-vertex polygons after the first invocation. Build them once at index construction and hold them for the process lifetime.

A Data Model That Can Explain Itself Later

Here is the schema I would defend in review. PostGIS 3.x, PostgreSQL 15 or newer.

CREATE EXTENSION IF NOT EXISTS postgis;
CREATE EXTENSION IF NOT EXISTS btree_gist;

CREATE TABLE service_zones (
    id            uuid PRIMARY KEY DEFAULT gen_random_uuid(),
    slug          text NOT NULL UNIQUE,
    display_name  text NOT NULL,
    zone_kind     text NOT NULL CHECK (zone_kind IN ('drive_time','manual','union')),
    depot_id      uuid NOT NULL REFERENCES depots(id),
    notes         text,
    created_at    timestamptz NOT NULL DEFAULT now()
);

The zone is an identity and a policy statement. It carries no geometry, because the geometry changes and the identity does not.

CREATE TABLE zone_versions (
    id              uuid PRIMARY KEY DEFAULT gen_random_uuid(),
    zone_id         uuid NOT NULL REFERENCES service_zones(id),
    revision        integer NOT NULL,
    geom            geometry(MultiPolygon, 4326) NOT NULL,
    geom_3tm        geometry(MultiPolygon, 3401)
                    GENERATED ALWAYS AS (ST_Transform(geom, 3401)) STORED,
    source_kind     text NOT NULL,
    provenance      jsonb NOT NULL,
    cost_seconds    integer,
    departure_bucket text,
    validity        tstzrange NOT NULL,
    authored_by     text NOT NULL,
    created_at      timestamptz NOT NULL DEFAULT now(),
    UNIQUE (zone_id, revision),
    CONSTRAINT geom_is_valid CHECK (ST_IsValid(geom)),
    CONSTRAINT no_overlapping_windows
        EXCLUDE USING gist (zone_id WITH =, validity WITH &&)
);

CREATE INDEX zone_versions_geom_gix ON zone_versions USING GIST (geom);
CREATE INDEX zone_versions_live_idx ON zone_versions USING GIST (validity);

Three things in that definition are doing real work. The generated projected column means planar measurement never pays a runtime transform and the index on it stays usable. The CHECK (ST_IsValid(geom)) constraint refuses malformed contour output at the door rather than three weeks later during an incident. The exclusion constraint — which needs btree_gist to mix an equality column with a range column — makes it structurally impossible for one zone to have two simultaneously active versions. That single line removes an entire category of ambiguity from the system.

CREATE TABLE coverage_checks (
    id                 bigserial PRIMARY KEY,
    requested_at       timestamptz NOT NULL DEFAULT now(),
    raw_address        text NOT NULL,
    normalized_address text,
    geocoder           text,
    geocode_quality    text NOT NULL,
    point              geometry(Point, 4326),
    zone_version_id    uuid REFERENCES zone_versions(id),
    decision           text NOT NULL
                       CHECK (decision IN ('covered','outside','undetermined')),
    decision_reason    text NOT NULL,
    boundary_margin_m  double precision,
    served_from_cache  boolean NOT NULL DEFAULT false,
    latency_ms         integer NOT NULL
);

CREATE INDEX coverage_checks_point_gix ON coverage_checks USING GIST (point);
CREATE INDEX coverage_checks_recent_idx
    ON coverage_checks (requested_at DESC)
    INCLUDE (decision, zone_version_id);

The foreign key from a check to a version is the whole point of the design. Someone reviewing a job from four months ago asks why the system said yes to an address it now says no to. Without that column the answer is an archaeology project. With it, you join two tables and print the polygon that was live at the time, alongside the provenance blob describing which graph extract and which departure bucket produced it. boundary_margin_m earns its place too — the distance to the nearest boundary tells you whether a decision was comfortable or a coin flip, and aggregating it reveals where the polygon needs attention.

Queries Worth Writing Down

Which zones cover a point right now, ordered by tightest drive-time budget:

SELECT z.slug, zv.id AS version_id, zv.cost_seconds
FROM zone_versions zv
JOIN service_zones z ON z.id = zv.zone_id
WHERE zv.validity @> now()
  AND zv.geom && ST_SetSRID(ST_MakePoint($1, $2), 4326)
  AND ST_Intersects(zv.geom, ST_SetSRID(ST_MakePoint($1, $2), 4326))
ORDER BY zv.cost_seconds ASC
LIMIT 1;

The && clause is redundant in the sense that ST_Intersects already applies it, but writing it explicitly documents the intent and survives refactors where somebody swaps the predicate for one that is not index-backed.

How far is this point from the edge of the zone, in metres:

SELECT ST_Distance(
         ST_Transform(ST_SetSRID(ST_MakePoint($1, $2), 4326), 3401),
         ST_Boundary(zv.geom_3tm)
       ) AS margin_m
FROM zone_versions zv
WHERE zv.id = $3;

Note the transform on the point and the pre-projected column on the polygon. Doing it the other way — transforming the polygon per row — is the single most common performance mistake in PostGIS code I review.

Area of a zone in square kilometres, which is a sanity check you should run on every generated version:

SELECT z.slug,
       zv.revision,
       round((ST_Area(zv.geom_3tm) / 1e6)::numeric, 1) AS area_km2
FROM zone_versions zv
JOIN service_zones z ON z.id = zv.zone_id
WHERE zv.validity @> now()
ORDER BY area_km2 DESC;

If a regeneration changes the area by more than a threshold you have chosen — twenty percent is a reasonable illustrative starting point — that is a signal, not a routine event. Something upstream moved. Hold the new version in a staging state and make a human look at it before it goes live.

Python: Prefilter, Exact Test, and a Ray Caster for Teaching

Production code should use Shapely and let GEOS do the arithmetic. Here is the structure I use, with the two-phase filter made explicit rather than hidden.

from dataclasses import dataclass
from typing import Optional, Sequence

from shapely.geometry import Point, shape
from shapely.geometry.base import BaseGeometry
from shapely.prepared import prep
from shapely.strtree import STRtree


@dataclass(frozen=True)
class ZoneVersion:
    version_id: str
    slug: str
    cost_seconds: int
    geometry: BaseGeometry


@dataclass(frozen=True)
class Verdict:
    covered: bool
    version_id: Optional[str]
    reason: str


class CoverageIndex:
    """Immutable spatial index over the zone versions live at build time."""

    def __init__(self, versions: Sequence[ZoneVersion]):
        self._versions = list(versions)
        self._geoms = [v.geometry for v in self._versions]
        self._prepared = [prep(g) for g in self._geoms]
        self._bounds = [g.bounds for g in self._geoms]
        self._tree = STRtree(self._geoms)
        self._position = {id(g): i for i, g in enumerate(self._geoms)}

    def _bbox_hit(self, idx: int, lon: float, lat: float) -> bool:
        minx, miny, maxx, maxy = self._bounds[idx]
        return minx <= lon <= maxx and miny <= lat <= maxy

    def evaluate(self, lon: float, lat: float) -> Verdict:
        probe = Point(lon, lat)
        candidates = self._tree.query(probe)
        if len(candidates) == 0:
            return Verdict(False, None, "no_bbox_candidate")

        hits = []
        for candidate in candidates:
            idx = self._resolve(candidate)
            if not self._bbox_hit(idx, lon, lat):
                continue
            if self._prepared[idx].contains(probe):
                hits.append(self._versions[idx])

        if not hits:
            return Verdict(False, None, "bbox_hit_exact_miss")

        best = min(hits, key=lambda v: v.cost_seconds)
        return Verdict(True, best.version_id, f"inside:{best.slug}")

    def _resolve(self, candidate) -> int:
        if isinstance(candidate, (int,)):
            return int(candidate)
        return self._position[id(candidate)]

The _resolve shim exists because Shapely 2.0 changed STRtree.query to return integer indices where 1.8 returned geometry objects. If you maintain a library that must span both, hide the difference in one place.

Now the pedagogical version. Nobody should ship this, but everybody who works on spatial systems should be able to write it from memory, because when a library gives you an answer you do not believe, this is how you find out who is wrong.

def point_in_ring(x: float, y: float, ring: Sequence[tuple]) -> bool:
    """Crossing-number test with a half-open edge rule.

    The ring is a closed sequence of (x, y) pairs. The half-open
    comparison (yi > y) != (yj > y) treats each edge as containing
    its lower endpoint and excluding its upper one, which makes a
    ray that passes exactly through a vertex count once, not twice
    and not zero times.
    """
    inside = False
    count = len(ring)
    j = count - 1
    for i in range(count):
        xi, yi = ring[i]
        xj, yj = ring[j]
        if (yi > y) != (yj > y):
            t = (y - yi) / (yj - yi)
            crossing_x = xi + t * (xj - xi)
            if crossing_x > x:
                inside = not inside
        j = i
    return inside


def point_in_polygon(x, y, exterior, holes=()) -> bool:
    if not point_in_ring(x, y, exterior):
        return False
    return not any(point_in_ring(x, y, hole) for hole in holes)

The half-open rule in that conditional is the entire trick. A naive implementation writes yi <= y < yj or yj <= y < yi and gets it right by accident; an implementation that writes min(yi, yj) <= y <= max(yi, yj) double-counts vertex hits and returns the wrong parity for any point whose latitude exactly matches a vertex. Since isochrone vertices come out of a regular grid, and grid latitudes are extremely round numbers, that case is far more common than a uniform-random model would predict.

Floating Point on the Boundary

Three degenerate cases, and what to do about each.

The point lies exactly on an edge. Mathematically the predicate is undefined unless you state a convention. Computationally, whether it lands on the edge depends on whether t and crossing_x round the way you hoped. Decide the policy explicitly — boundary counts as inside is the humane choice for a service area — and implement it as a separate test rather than hoping the crossing count works out. In PostGIS that means ST_Covers rather than ST_Contains; ST_Contains returns false for boundary points and ST_Covers returns true, and confusing the two produces a bug that appears roughly one time in ten thousand and is never reproducible on demand.

The ray passes through a vertex. Handled by the half-open rule above. Verify it with a fixture: a point at exactly the same y as a polygon vertex, once with the vertex to the left and once to the right.

Two computations disagree. This is the one that ruins afternoons. The same address evaluated twice returns different answers because one path went through a projection and the other did not, or because a value round-tripped through JSON at fifteen significant digits instead of seventeen. Determinism is a property you have to engineer.

Three practices buy it:

1. Snap on ingest. ST_SnapToGrid(geom, 1e-7) quantizes coordinates to about 1.1 centimetres of latitude. Every downstream computation then starts from the same numbers, and the crossing arithmetic becomes reproducible across library versions.
2. Serialize with full precision. ST_AsGeoJSON(geom, 15) and a JSON encoder that does not truncate. A double needs 17 significant decimal digits for a lossless round trip; most defaults give you fewer.
3. Refuse to answer near the boundary. Compute the margin. If the point sits within a tolerance band of the edge — five metres is a defensible illustrative choice — return undetermined and let a human resolve it. This is not a cop-out. It is an honest statement that the polygon's own accuracy is coarser than the question being asked, and it converts a silent wrong answer into a visible, countable event.

Underneath all of this, if you need genuinely exact orientation tests, the answer is adaptive-precision arithmetic — Shewchuk's robust predicates, which JTS implements with double-double arithmetic and GEOS inherits. Knowing it exists is usually enough; needing it means your tolerances are wrong.

Geocoding Failures and Failing Closed

Coverage evaluation has a dependency most designs underweight: turning text into coordinates. If geocoding is wrong, everything downstream is confidently wrong.

Calgary makes this harder than average because of quadrants. The city is divided into NW, NE, SW, and SE, and street numbering restarts in each. "123 5 Street" is four different places. A form that lets a customer omit the quadrant, or a normalizer that drops it, will hand the geocoder an ambiguous string, and the geocoder will pick one — usually the most populous match, usually silently. The van is then dispatched to a location twelve kilometres from the vehicle. Quadrant is not optional metadata; treat a missing quadrant as a validation failure and resolve it before geocoding.

The normalizer should also handle unit designators, since a driveway visit needs the building and a bay number needs the unit; the abbreviation zoo of ST, AV, DR, CR, BV, TR, WY, GA, LI, PZ; and directional suffixes that some sources write as N.W. and others as NW.

Then fail closed. Every geocoder returns a quality signal, and the tiers matter:

Quality	Meaning	Policy
Rooftop / parcel	Resolved to the property	Evaluate normally
Interpolated	Estimated along a street segment	Evaluate, widen the margin band
Street centroid	Street matched, number did not	Undetermined, route to review
Locality only	City or town matched	Undetermined
No match	Nothing	Undetermined

An interpolated result on a long rural road can be hundreds of metres off, which is exactly where the polygon edge tends to be. Treating "interpolated" as equivalent to "rooftop" is a quiet source of bad dispatches, and the correction is small: multiply the boundary tolerance by the positional uncertainty of the geocode tier.

Failing closed means an unresolved address never returns covered. It returns undetermined with a reason code, and something in the workflow — a queue, a flagged record, a note on the service request — routes it to a person. Customers near a boundary are also frequently near an edge case in the data, and rural addresses on range roads are both. For a business whose service area genuinely extends past the city limits, that queue is not an exception path; it is a normal part of Tuesday.

Caching on Rounded Coordinates

Once volume grows, evaluating the same location repeatedly is waste. The obvious cache key is the coordinate pair, and the obvious refinement is to round it so that near-identical points share an entry.

Round to four decimal places and, at Calgary's latitude, a cell is roughly 11.1 metres north-south by 7.0 metres east-west. Five decimals gives about 1.1 by 0.7 metres. The tradeoff is direct: coarser cells mean a higher hit rate and a larger maximum positional error introduced by the key itself.

The failure mode is specific. If a rounding cell straddles a zone boundary, two genuinely different addresses collapse to one key, and whichever was evaluated first determines the answer for both. One of them is wrong, and the record shows a clean, plausible, cached decision.

The mitigation is to make cacheability conditional on the margin you already compute:

CELL_DEG = 1e-4
CELL_DIAGONAL_M = 13.2  # illustrative, at ~51 deg N


def cache_key(lon: float, lat: float, policy_id: str) -> str:
    qx = round(lon / CELL_DEG)
    qy = round(lat / CELL_DEG)
    return f"cov:{policy_id}:{qx}:{qy}"


def maybe_store(result, margin_m: float) -> bool:
    if result.reason == "undetermined":
        return False
    return margin_m is not None and margin_m > CELL_DIAGONAL_M

Comfortable decisions get cached; marginal ones are recomputed every time. Since the marginal population is small, the hit rate barely moves and the correctness problem disappears. Include the policy or zone-set identifier in the key so that publishing a new zone version invalidates by construction rather than by remembering to flush.

When the Routing Provider Is Down

Because isochrones are precomputed, a routing outage cannot break coverage evaluation. That is most of the value of the offline design, and it is worth stating explicitly to whoever is reviewing the architecture.

What an outage does break is regeneration, and the failure is quiet. Polygons drift out of date, nothing alarms, and three weeks later somebody notices the new interchange is not represented. Alert on generation age, not on generation errors.

For the routing that does happen live — estimating travel time to a specific address so a dispatcher can size the day — build a ladder:

1. Primary. The routing engine with a strict timeout. Two seconds is a reasonable illustrative ceiling for an interactive path.
2. Secondary. A cached matrix of depot-to-cell travel times at a coarse grid resolution. Stale by hours, accurate within a few minutes, and available instantly.
3. Tertiary. Straight-line distance multiplied by a direction-dependent circuity factor derived from historical data. Crude, but it degrades gracefully and it is better than an error page.
4. Floor. Return the polygon-based coverage verdict with the estimate marked provisional, and let the human decide.

Each tier tags its output with the method used, and that tag lands in the service record. When you later analyze which estimates were poor, you can separate "the model is bad" from "the model was not running."

The same discipline applies to the geocoder. A second provider behind a circuit breaker costs a day to integrate and removes a single point of failure from the most upstream step in the chain. For urgent work — a flat repair on a vehicle stranded in a lot, or the kind of roadside situation where minutes are the whole story — a degraded answer delivered now beats a perfect answer delivered after the customer has given up.

Testing Coverage Without Brittleness

Spatial tests fail badly when written naively, because the natural assertion — "this coordinate is inside this polygon" — breaks whenever the polygon is regenerated, which is often. The tests then get deleted, and coverage logic ends up with the worst test suite in the codebase despite being one of its riskiest components.

Three layers that hold up.

Golden points. A fixture file of coordinates with expected verdicts and a human-written rationale. Choose points that are unambiguous by construction — deep inside the core, well outside in open country, on the far side of an obvious barrier — not points scraped from a map near an edge.

GOLDEN = [
    ("depot_yard",        -114.0719, 51.0447, "covered"),
    ("inner_city_core",   -114.0630, 51.0450, "covered"),
    ("north_ring_edge",   -114.0100, 51.1800, "covered"),
    ("far_east_prairie",  -112.4000, 50.9500, "outside"),
    ("mountain_park",     -116.0300, 51.1700, "outside"),
]

Assert the verdict and the reason code. If a regeneration flips a golden point, the build fails and a person looks at the polygon — which is the outcome you want.

Property tests. These survive polygon churn because they assert invariants rather than facts. The most valuable one is projection stability: reproject a point from 4326 into the working CRS and back, and require that the containment verdict is unchanged and the coordinate agrees to within a tight tolerance.

from hypothesis import given, strategies as st
from pyproj import Transformer

FWD = Transformer.from_crs("EPSG:4326", "EPSG:3401", always_xy=True)
INV = Transformer.from_crs("EPSG:3401", "EPSG:4326", always_xy=True)


@given(
    lon=st.floats(min_value=-115.4, max_value=-112.6),
    lat=st.floats(min_value=49.8, max_value=52.2),
)
def test_reprojection_round_trip_preserves_verdict(lon, lat):
    x, y = FWD.transform(lon, lat)
    lon2, lat2 = INV.transform(x, y)
    assert abs(lon2 - lon) < 1e-9
    assert abs(lat2 - lat) < 1e-9
    assert INDEX.evaluate(lon, lat).covered == INDEX.evaluate(lon2, lat2).covered

Other invariants that pay for themselves: every stored geometry is valid; no two active versions of one zone overlap in time; a point inside a tighter drive-time band is inside every looser band from the same origin; and the same input evaluated twice in one process yields byte-identical output.

Boundary tests without brittleness. Do not assert what happens on the line. Assert monotonicity around it. Take a boundary vertex, step fifty metres along the inward normal, and require covered. Step fifty metres outward and require not covered. Between those, require only that the answer is stable and that the margin computation reports a value below the tolerance. You are testing the shape of the decision function rather than the coordinates of the polygon, which is what you actually care about.

One more thing worth automating: an area-delta check between consecutive versions, run in CI against a fixture polygon set. Contouring parameters get tuned, somebody changes a generalization value, and the polygon quietly loses a limb. A percentage-change assertion catches it in the pipeline instead of in dispatch.

Seasonality as a Time-Windowed Zone Version

Alberta has two weeks in spring and two in fall when the entire province changes its tires at once. Demand during the changeover peak is a different distribution from the rest of the year, not merely a larger one, and it breaks any assumption that coverage is a static property.

The arithmetic is unforgiving. A van that completes a fixed number of jobs a day at forty minutes of driving between them completes meaningfully fewer than one averaging twenty. During peak weeks the marginal long-distance job does not add revenue; it removes two closer jobs from the day. So the rational operating policy is a smaller effective service area precisely when demand is highest.

The wrong way to implement that is to edit the polygon. Somebody opens the geometry, shrinks it, and in April somebody else expands it again. History is destroyed, the audit trail lies, and a service record from March cannot be explained.

The right way is a second zone version with a bounded validity window, competing for the same zone identity:

INSERT INTO zone_versions (
    zone_id, revision, geom, source_kind, provenance,
    cost_seconds, departure_bucket, validity, authored_by
)
VALUES (
    :zone_id,
    (SELECT COALESCE(MAX(revision), 0) + 1
       FROM zone_versions WHERE zone_id = :zone_id),
    ST_Multi(ST_GeomFromGeoJSON(:geojson)),
    'isochrone',
    :provenance::jsonb,
    1500,                          -- 25 minutes, tightened from 30
    'weekday_pm_peak',
    tstzrange('2026-10-05 00:00+00', '2026-10-26 00:00+00', '[)'),
    'seasonal-policy-job'
);

The exclusion constraint enforces the rest. If the peak window overlaps the standing version's window, the insert fails, and you are forced to close the standing window explicitly rather than creating a silent ambiguity. That is the constraint doing its job: converting a policy question into a schema error at the moment somebody tries to be careless.

A few consequences fall out of this design that I did not anticipate the first time:

• Quotes issued before the window opens must be honoured against the version live at issue time, which means the quote stores a zone_version_id and the fulfilment path reads it rather than re-evaluating.
• The tightened polygon is a different isochrone, not a scaled copy. Regenerate it with the peak departure bucket; a negative buffer of the standing polygon shrinks it uniformly, which is exactly the isotropic mistake the whole model exists to avoid.
• Announcing the window in advance is cheaper than declining requests during it. Whatever surfaces coverage to customers — the self-serve scheduling page, a dispatcher's screen — should read the version that will be live on the requested date, not the one live right now.

The same mechanism handles other temporal policies without new machinery: a reduced winter zone when a road closes seasonally, an expanded window when a second van joins the rotation, a temporary exclusion around a construction detour. Every one of them is a row with a validity range, and every one of them is explainable afterward.

Latency, Monitoring, and What to Alert On

Coverage evaluation sits on an interactive path, so it needs a latency budget and observability that matches.

The realistic decomposition for a cold request:

Stage	Typical share	Notes
Address normalization	small	Pure CPU, string work
Geocoding	dominant	Network, external provider
Index query and exact test	small	Microseconds with prepared geometry
Margin computation	small	One projected distance
Persisting the check record	small	Single insert, can be async

Geocoding is the entire story. Everything else is noise. That has two implications: cache geocodes aggressively and separately from coverage verdicts, since a normalized address maps to a coordinate far more stably than a coordinate maps to a verdict; and instrument the provider call with its own timer so a provider slowdown is not misattributed to your spatial code.

Metrics I would page on, or at least chart on a wall:

• p50, p95, and p99 end-to-end evaluation latency, split by cache-hit status. Averaging across both hides everything.
• Geocode failure rate by tier. A rise in "street centroid" results usually means a normalizer regression, not a provider problem.
• The undetermined rate. This is your honesty metric. A sudden climb means the geocoder degraded or a new polygon has a boundary running through a populated area.
• Boundary-margin distribution. If a large share of requests cluster within a few hundred metres of an edge, the polygon is drawn through demand rather than around it, and the budget should move.
• Decision-flip rate across version transitions. Replay a recent sample of points against the incoming version before it goes live; a flip count outside the expected range blocks publication.
• Age of the newest zone version per origin, compared against the age of the graph extract.

That last one has caught more real problems than the rest combined, because stale geometry produces answers that look completely normal.

What I Would Build Differently

Three things I got wrong on the first pass, offered as a shortcut past them.

I started with a single polygon per origin and a boolean verdict. Bands should have existed from the beginning. The moment dispatch wants to sequence a day, "inside" is not enough information; they need the drive-time tier, and retrofitting bands meant reworking every consumer. Model concentric budgets — twenty, thirty-five, fifty minutes — from the start, and let the boolean be a derived view over the innermost band that policy currently accepts.

I treated the peak-season constraint as an operational note rather than a data-model requirement. It is a data-model requirement. Anything that changes coverage over time belongs in the versioning mechanism, and discovering that in the middle of an October rush is an unpleasant way to learn it.

I under-invested in the undetermined path. The first version had two outcomes, and every ambiguous case was forced into one of them, which meant the system was confidently wrong at exactly the addresses where being wrong mattered most — acreages, new subdivisions, rural routes. A third state plus a review queue turned an invisible error rate into a visible work item, and the work item turned out to be small.

The broader lesson generalizes past tires. A service area is not a fact about your business; it is a model of a cost surface that changes with the road network, the hour, the season, and the number of vehicles you are running. Postal codes model none of that. A radius models one dimension of it badly. Versioned isochrones with explicit provenance model most of it, and — more importantly — they tell you when they are unsure.

If you want to see the shape of the real-world constraints this design is serving, the operational side is documented plainly: the fleet programs that drive most of the multi-vehicle routing, the commercial work that pushes the polygon outward toward industrial parks, and the Calgary yard that anchors every isochrone in this article. The software exists to answer one question honestly, and the geometry is how it does that.