Chunking County-Scale Runs by Spatial Tile

Splitting a county run into tiles is the cheapest way to make it parallel, resumable and bounded in memory — and the easiest way to change the answer without noticing. Every operation in a compliance pipeline that looks at a parcel’s neighbours (setback encroachment, sliver detection, nearest-feature distance) reads across the tile edge, so a naive partition silently drops the relationships that cross it. This guide builds tiles that balance work rather than area, buffers each tile so edge parcels see their neighbours, and asserts that the partitioned run returns exactly what the whole-county run returns. It is the partitioning half of batch processing optimization.

Prerequisites

Step-by-step

Step 1: Determine the interaction distance before choosing tiles

The buffer that makes tiling safe is determined by the rules, not by the geometry. Work it out explicitly and let it fail loudly if a rule exceeds it.

def interaction_distance(rules) -> float:
    """The furthest any rule reaches from the parcel it is evaluating."""
    reaches = []
    for r in rules:
        if r["kind"] == "setback":
            reaches.append(r["threshold"])
        elif r["kind"] == "proximity":
            reaches.append(r["buffer_ft"])
        elif r["kind"] == "adjacency":
            reaches.append(r.get("search_ft", 0.0))
    return max(reaches) if reaches else 0.0

HALO_FT = interaction_distance(RULES) * 1.5     # margin for geometry that overruns

The multiplier is there because a rule’s threshold is measured from a parcel boundary, and a parcel near a tile edge may itself extend past that edge. A halo of 1.5× the largest threshold covers both effects on any normal fabric, and a check in step 5 confirms it did.

The halo comes from the rules, not from the tileEach rule contributes a reach; the halo is the largest of them with a margin, and a tile reads everything within it while evaluating only what it owns.Setback rules — reach is the thresholdtens of feetProximity rules — reach is the bufferhundreds of feetHalo = largest reach x 1.5margin for parcels overrunning the edgeAsserted at 90% during the runwarning before the partition changes answers
A larger tile does not make a rule look further. Deriving the halo from the tile size is how a partition silently starts changing answers.

Step 2: Balance the tiles by parcel count, not by area

A uniform grid over a county gives one tile with 40,000 parcels and forty tiles with none, so the run takes as long as its worst tile. Splitting recursively on the median coordinate — a KD-tree split — gives tiles of similar population and wildly different sizes, which is exactly what balances the work.

import numpy as np
from shapely.geometry import box

def balanced_tiles(gdf, max_per_tile=5000):
    """Recursive median splits on the longer axis until every tile is small enough."""
    def split(sub, bounds):
        if len(sub) <= max_per_tile:
            return [(box(*bounds), sub.index)]
        minx, miny, maxx, maxy = bounds
        pts = sub.geometry.representative_point()
        if (maxx - minx) >= (maxy - miny):
            cut = float(np.median(pts.x))
            left, right = pts.x < cut, pts.x >= cut
            return (split(sub[left], (minx, miny, cut, maxy))
                    + split(sub[right], (cut, miny, maxx, maxy)))
        cut = float(np.median(pts.y))
        low, high = pts.y < cut, pts.y >= cut
        return (split(sub[low], (minx, miny, maxx, cut))
                + split(sub[high], (minx, cut, maxx, maxy)))

    return split(gdf, tuple(gdf.total_bounds))

The tiles tessellate exactly — every point of the county is in one tile and no point is in two — which is the property the assignment rule in step 4 depends on. Where an H3 index is already on the parcels, its coarser resolutions make a serviceable alternative partition; the trade-off is that H3 cells balance by area rather than by population, so they inherit the problem this step exists to solve unless population is checked.

Step 3: Read each tile with its halo

The tile defines which parcels are evaluated; the halo defines which are available. Reading with a bounding-box filter pushes the selection into the file format, which is where GeoParquet earns its place — see converting parcel shapefiles to GeoParquet for batch runs.

import geopandas as gpd

def read_tile(path, tile_geom, halo_ft):
    haloed = tile_geom.buffer(halo_ft)
    gdf = gpd.read_parquet(path, bbox=haloed.bounds)
    gdf = gdf[gdf.intersects(haloed)]
    # Evaluate only what the tile owns; the rest is context.
    gdf["in_tile"] = gdf.geometry.representative_point().within(tile_geom)
    return gdf

Reading by bounding box rather than filtering after a full read is the difference between a tile that loads in under a second and one that loads the whole county forty times.

Step 4: Assign every parcel to exactly one tile, deterministically

The halo means a parcel appears in several tiles. Exactly one of them must own it, and ownership must not depend on which tile finished first.

def evaluate_tile(gdf, rules):
    """Evaluate everything for context; emit results only for parcels this tile owns."""
    results = run_rules(gdf, rules)               # uses the full haloed frame
    return results[results["parcel_id"].isin(gdf.loc[gdf["in_tile"], "parcel_id"])]

Ownership by representative point is deterministic because the tiles tessellate and a representative point is always inside its own parcel: every parcel’s point falls in exactly one tile regardless of evaluation order. Ownership by intersection would not be — a parcel straddling an edge intersects two tiles and would be emitted twice.

Which tile owns a parcel that straddles an edgeEvery parcel appears in several tiles because of the halo. Exactly one must emit a result for it, and not the one that happened to finish first.Does the parcel’srepresentative point fallinside this tile?noContext only — evaluate, emit nothingthe owning tile will emit ityesThis tile owns it — emit the resultevaluated with the full haloed neighbourhoodTiles tessellate, so every parcel has exactly one ownerdeterministic regardless of which tile ran first
Ownership by intersection emits the parcel twice; ownership by representative point emits it exactly once, whatever the execution order.

Step 5: Assert the halo was big enough

The halo is an assumption about the rules, and assumptions about rules should be checked at runtime rather than believed.

def check_halo(gdf, results, tile_geom, halo_ft):
    """Nothing a tile relied on may have come from beyond the halo."""
    used = gdf[gdf["parcel_id"].isin(results["neighbour_id"].dropna())]
    max_reach = used.geometry.distance(tile_geom.boundary).max()
    if max_reach > halo_ft * 0.9:
        raise RuntimeError(
            f"a rule reached {max_reach:.1f} ft from the tile edge; halo is {halo_ft:.1f} ft")

Failing at 90% of the halo rather than at 100% gives warning before the partition actually starts changing answers. This check is what makes the tiling defensible: without it, the correctness of a county run rests on somebody having remembered every rule’s reach.

Step 6: Write per-tile output that survives a retry

A tile that fails must be re-runnable without corrupting what succeeded, which means the output has to be addressed by tile rather than appended to a shared file.

import os

def write_tile_result(results, out_dir, tile_id, run_id):
    path = f"{out_dir}/run={run_id}/tile={tile_id:05d}/results.parquet"
    os.makedirs(os.path.dirname(path), exist_ok=True)
    tmp = path + ".tmp"
    results.to_parquet(tmp, index=False)
    os.replace(tmp, path)      # atomic: a reader never sees a partial file
    return path
Output addressed by tile is output that survives a retryEach tile writes to its own path under the run, via a temporary name and an atomic rename, so a failed tile is re-runnable without touching what succeeded.Write to a temporary namea reader never sees a partial fileRename atomically into placerun id and tile id in the pathA failed tile leaves no filepresence is the completion recordResume by listing what is missingno bookkeeping to keep in sync
Re-running only the tiles with no result file is correct by construction — which is what makes a county run resumable rather than restartable.

Writing to a temporary name and renaming is what makes a tile’s output atomic. It is the same discipline the queueing in building async rule queues for batch zoning validation relies on, and it turns a failed run into a resumable one: re-running only the tiles with no result file is correct by construction.

Verification

The partitioned run must return exactly the unpartitioned run. That is testable on a sample area and should be a standing test.

whole = run_rules(gpd.read_parquet(PARCELS), RULES).sort_values("parcel_id")
tiled = pd.concat([pd.read_parquet(p) for p in tile_outputs]).sort_values("parcel_id")

assert len(tiled) == len(whole), f"tiled emitted {len(tiled)} vs {len(whole)}"
assert tiled["parcel_id"].is_unique, "a parcel was emitted by two tiles"
pd.testing.assert_frame_equal(
    tiled.reset_index(drop=True), whole.reset_index(drop=True), check_like=True)

# Balance: the slowest tile bounds the run.
sizes = [len(pd.read_parquet(p)) for p in tile_outputs]
print(f"{len(sizes)} tiles, {min(sizes)}–{max(sizes)} parcels, "
      f"imbalance {max(sizes) / (sum(sizes) / len(sizes)):.2f}x")

An imbalance ratio near 1.2 or below means the median splits did their job. Above 2, the run is waiting on one tile and the max_per_tile threshold wants lowering.

Common Pitfalls

  • Tiling without a halo. Every neighbour relationship crossing a tile edge disappears, and the result looks complete because every parcel still got a verdict.
  • Assigning parcels by intersection. A parcel straddling an edge is emitted by both tiles, so it appears twice with possibly different results.
  • Uniform grids over a county. Population per tile varies by orders of magnitude and the run takes as long as its densest tile.
  • Appending all tiles to one file. A retried tile then duplicates its rows, and a failed run cannot be resumed without redoing everything.
  • Deriving the halo from the tile size. It comes from the rules’ reach. A larger tile does not make a rule look further.
  • Filtering after a full read. Reading the whole county per tile turns a partitioning optimisation into a slowdown proportional to the tile count.

Frequently Asked Questions

How large should a tile be?

Small enough that a tile plus its halo fits comfortably in one worker’s memory, large enough that the halo is a small fraction of the tile. A few thousand parcels usually satisfies both; below about a thousand, the halo starts dominating the read and the overhead grows.

Does the halo make the run substantially more expensive?

It makes each tile read more parcels than it evaluates, and the overhead is roughly the ratio of the halo’s ring area to the tile’s area. For a 5,000-parcel tile and a 200-foot halo, that is a few percent — negligible against the parallelism gained. It becomes significant only when tiles get small or a rule has a very long reach.

Can I tile by administrative boundary instead?

You can, and it is attractive because the outputs align with how results are consumed. The costs are that municipalities vary in size by orders of magnitude, which reintroduces the balance problem, and that boundary-crossing relationships are exactly the ones a jurisdiction-based partition hides — see jurisdictional boundary and precedence resolution.

How does this interact with Dask?

Well, but the halo has to be explicit — Dask-GeoPandas partitions spatially and will happily give you partitions with no overlap. The pattern that works is to compute tiles as above, then map the per-tile function over them rather than relying on the framework’s own partitioning. Parallelizing parcel validation with Dask-GeoPandas covers the execution side.

Should tiles be recomputed each run?

Only when the fabric changes materially. Stable tile boundaries make run-to-run comparison much easier, because a result diff is then per-tile and a changed tile localises the change. Recomputing the partition every run means every result moves, and the diff tells you nothing — which matters for detecting rule drift between pipeline releases.

What if a rule genuinely has unbounded reach?

Then it cannot be tiled and should be run as a separate whole-county pass. Nearest-feature queries are the usual case: the nearest protected feature could in principle be arbitrarily far, though in practice a generous cut-off plus a fallback for the parcels that find nothing works well — the approach in measuring distance to the nearest protected feature. Trying to force such a rule into the tiled pass is how a halo assumption gets quietly violated.

Part of: Batch processing optimization