A spatial grid may shrink the candidates, never the lookup's answer

11 November 2026 · 11 min read

Delivery coverage is a point-in-polygon question you ask over and over against a catalogue that changes rarely, so any implementation that survives contact with traffic grows an index sooner or later. The danger is that the index stops being an optimisation and becomes part of the answer. Assign a polygon to one cell instead of every cell it touches, or hand candidates back in cell order instead of catalogue order, and the index has quietly changed which store serves an address. In our delivery-zones library (shipmindlabs/delivery-zones) the grid is held to a single invariant (a lookup returns what a linear scan would return, in the same order) and everything else about it is negotiable.

That constraint is not a style preference. Coverage polygons are a business artefact: someone drew a ring around a district, someone else set a fee tier on it, and a dispatcher reads the result. If the index can change the result, every disagreement about coverage turns into an argument about geometry internals instead of an argument about the catalogue. So we allow the index exactly one effect, fewer polygons tested, and no others.

The invariant, written down as a test#

The invariant is checkable at all because the slow path is a real, public, supported operation. point_in_polygon tests one polygon once, and Zone exposes its polygons, so a linear scan is a few lines of caller code rather than a reimplementation of the library's internals:

from delivery_zones import ZoneIndex, point_in_polygon


def linear_scan(zones, point):
    matches = []
    for zone in zones:
        if any(point_in_polygon(polygon, point) for polygon in zone.polygons):
            matches.append(zone)
    return tuple(matches)


index = ZoneIndex(zones)
assert index.zones_containing(point) == linear_scan(index.zones, point)

That is what tests/test_index_oracle.py is for. We do not test the grid by asserting that a particular polygon landed in a particular cell, since that would pin the implementation and prove nothing about correctness. We test it against the scan: for a catalogue and a spread of positions, the two must agree, element for element, in order. A bucketing bug then shows up as a disagreement about zones, which is the language the rest of the system speaks.

Having an oracle also changes what the grid is permitted to be. Because the answer is defined elsewhere, the cell layout is free to be crude. It does not need to be a quadtree, it does not need to balance, and it does not need to handle degenerate catalogues gracefully in any clever way. It only has to avoid losing a candidate.

The grid filters; the same geometry decides#

There is one containment implementation, and both entry points go through it. _PreparedPolygon keeps the polygon's bounding box and turns the exterior and each hole into a _PreparedRing, which stores edges as (lat_start, lat_end, lon_start, inverse_slope) and drops horizontal edges outright. A horizontal ray never crosses them, which is also what makes the slope division safe to precompute:

    @classmethod
    def prepare(cls, points: Sequence[Point]) -> "_PreparedRing":
        edges: list[tuple[float, float, float, float]] = []
        previous = points[-1]
        for current in points:
            if previous.lat != current.lat:
                slope = (previous.lon - current.lon) / (previous.lat - current.lat)
                edges.append((current.lat, previous.lat, current.lon, slope))
            previous = current
        return cls(tuple(edges))

    def wraps(self, lon: float, lat: float) -> bool:
        inside = False
        for lat_start, lat_end, lon_start, slope in self.edges:
            if (lat_start > lat) != (lat_end > lat):
                if lon < lon_start + (lat - lat_start) * slope:
                    inside = not inside
        return inside

Containment is then bounding box, exterior, holes, in that order:

    def contains(self, point: Point) -> bool:
        if not self.bbox.contains(point):
            return False
        if not self.exterior.wraps(point.lon, point.lat):
            return False
        return not any(hole.wraps(point.lon, point.lat) for hole in self.holes)

point_in_polygon prepares one _PreparedPolygon and calls that method. A lookup through ZoneIndex calls the same method on polygons prepared at construction time. The grid sits in front of it and has no vote on the outcome. That is the whole architecture: preparation is shared, and the index contributes only a shorter list of things to prepare once and then test.

It also means the candidate set is deliberately sloppy. Cells are indexed by bounding box, so a candidate is a polygon whose envelope overlaps the query cell, not a polygon that contains the point. False positives are normal and expected. contains throws them out, starting with the cheap bounding box check that the grid has already made likely to pass.

Every polygon in every cell its envelope touches#

Bucketing is where indexes usually acquire their bugs, because the tempting representation is one cell per polygon. By centroid, say, which keeps the cell lists tidy and is wrong at every border. A district whose envelope straddles four cells is only findable from one of them, and the addresses that break are exactly the addresses near the edges of a zone, which is where coverage disputes come from in the first place.

So we register a polygon in every cell its envelope spans: from the column of box.min_lon to the column of box.max_lon, inclusive, and likewise for rows. A wide zone is stored many times over, and that redundancy is deliberate, because it is how the filter stays conservative. The axis mapping clamps, so a position on or beyond the envelope boundary resolves to a real cell rather than an out-of-range one:

def _cell_axis(value: float, low: float, high: float, side: int) -> int:
    span = high - low
    if span <= 0.0:
        return 0
    return min(side - 1, max(0, int((value - low) / span * side)))

The zero-span branch is not defensive padding. A catalogue can legitimately have no extent on an axis (one polygon collapsed to a line, or a set of zones sharing a longitude), and division would be the only thing standing between that catalogue and a lookup. Returning column zero folds the whole axis into one band, which is slower and still correct.

The envelope of the catalogue itself is the union of the polygon boxes, and it is degenerate when there is nothing to bound:

def _envelope(boxes: Sequence[BoundingBox]) -> BoundingBox:
    if not boxes:
        return BoundingBox(0.0, 0.0, 0.0, 0.0)
    return BoundingBox(
        min(box.min_lon for box in boxes),
        min(box.min_lat for box in boxes),
        max(box.max_lon for box in boxes),
        max(box.max_lat for box in boxes),
    )

An empty index has no cells either, so the degenerate box off the coast of Africa is never consulted. The cell lookup finds nothing to test regardless.

Outside the envelope is an answer, not a search#

The one place we let the index answer without geometry is outside its own envelope. A position beyond _bounds cannot be inside any polygon whose box contributed to that envelope, so candidate selection yields nothing and no crossing count is ever run. covers is the same candidate walk under any, so it short-circuits on the first containing polygon rather than collecting the rest.

That early exit is not a micro-optimisation for the request path alone. The coverage tools in coverage.py hammer the same entry points: scan_coverage walks a lattice of positions and asks index.covers(point) at each one, and scan_ambiguity runs zones_containing with a policy over the same lattice to find the positions where the policy elects nobody. Those scans are bounded by max_samples for a reason (the sample count grows with the square of the inverse spacing), and they are only usable as a routine check because the per-sample cost is a short cell list rather than the whole catalogue.

Declaration order survives the grid#

We flatten the polygons once, at construction, into pairs of (position, _PreparedPolygon) in zone order. The position is the zone's index in the catalogue, and the cells hold slots into that flattened tuple. Registration walks the polygons in order and candidates come back in slot order, so the candidates for any cell emerge in declaration order, not cell order, not insertion-into-cell order, not whatever order a hash happened to produce.

That is what makes the file the tie-breaker of last resort. A zone with several polygons is deduplicated by a seen set of positions, so it appears once before the matches reach the policy, and the default policy keeps all of them. Only then does a business rule get to choose:

serving = index.zones_containing(address, policy=by_priority(TieBreak.FIRST))

TieBreak.FIRST means the first zone in the catalogue among the tied ones. If the grid reordered candidates, that would silently become the first zone in whichever cell the address happened to fall in: a rule that is stable per address and arbitrary across addresses, which is probably the worst kind of wrong, because it looks deterministic in every test you write by hand.

Sizing by the square root, capped#

We pick the grid side from the data rather than from configuration:

        self._side = max(1, min(_MAX_GRID_SIDE, isqrt(len(self._polygons))))

With side taken as the integer square root of the polygon count, the grid has roughly as many cells as there are polygons, which keeps average cell lists short without making the cell table larger than the catalogue it indexes. _MAX_GRID_SIDE is 64, so the table never exceeds a few thousand cells. Past that point the duplication from wide-envelope zones dominates anyway, and more cells buy memory and build time instead of speed. The max(1, ...) floor matters at the other end: a one-polygon catalogue gets a single cell, so a small catalogue pays nothing for the grid beyond one bounding box test.

There is no tuning knob for this, deliberately. A knob would be a way for a caller to change lookup cost, and the only honest thing a caller could do with it is measure, which the library cannot do for them and they will not do either.

A snapshot, not a container#

ZoneIndex has no add, no remove, and no reindex. It prepares every ring once in the constructor and holds the zone tuple, the flattened polygons, the envelope, the side and the cells behind __slots__. When the catalogue changes, you build a new index and swap the reference.

This is the structural half of the invariant. The claim a lookup returns what a linear scan over index.zones returns is only verifiable if the zones and the prepared geometry cannot drift apart, and a mutable index is exactly a machine for making them drift. A zone appended without rebucketing, or rebucketed without reflattening, produces a lookup that is wrong in a way no type checker will catch. Immutability removes the state where the two disagree.

What it costs to run#

Memory is the edge lists for every ring, plus one slot reference per cell each polygon's envelope spans. Build time is one pass over the rings (a division per non-horizontal edge, one bounding box per polygon), paid at load or reload, not per request. A lookup is one envelope check, one cell list, and then the same bounding box plus crossing count that the linear scan would have run, on fewer polygons.

The cost that is easy to miss is the rebuild. A catalogue that changes often turns index construction into a recurring cost, and the shape of the library pushes that cost to where it belongs: a deployment, a configuration reload, a scheduled refresh. In exchange, the request path never sees a partially updated index, and a coverage scan over a lattice reflects one consistent catalogue rather than whatever was true as the walk progressed.

The library is pre-alpha and the public API is still moving, but this part is settled, because it is the part that makes the rest reviewable. A grid that only filters can be checked against a scan. A grid that answers has to be trusted.