Spatial Epi Public Health GIS

Streaming Large Boundary Files Without Exhausting Memory

This guide, part of Spatial Data Types & Formats in Production Epidemiology Pipelines, covers the practical problem of joining case data to reference geography that is too large to load: census blocks for a large state, national parcel data, a building footprint layer.

Problem Context & Constraints

The standard pattern — read the whole boundary file into a GeoDataFrame, then spatially join — works until it does not. A state’s census blocks run to hundreds of thousands of polygons with detailed coastlines; a county parcel layer can exceed a million. Loaded fully, with geometry objects, spatial index and the join’s intermediate results, the working memory can exceed what an analysis server has, and the operating system kills the process partway through a nightly run.

The pipeline’s actual need is almost always smaller than the file. A day’s cases touch a few hundred blocks, not all of them. An exposure join for one county needs that county’s parcels. The waste comes from reading everything to use a fraction.

Three constraints shape a streaming approach. Filter at the source — read only the features that can matter, by bounding box or attribute, using the format’s own filtering so unwanted geometry is never deserialised. Process in bounded chunks when the needed subset is still large, holding only one chunk’s geometry at a time. And handle chunk edges correctly: a point near the edge of one tile may belong to a polygon that mostly lies in the next, so tiles need a halo, and each point must be assigned exactly once.

Memory Over Time for Three Read Strategies Memory use plotted over the course of a join for three strategies. Reading the whole file climbs steadily until it exceeds the server's memory limit and the process is killed. Reading only the bounding box of the day's cases reaches a much lower plateau and completes. Reading that bounding box in tiles produces a sawtooth pattern that never rises above one tile's worth of geometry plus the running result, and completes regardless of how large the file is. Read what the job needs, a tile at a time server memory limit whole file: killed bounding box only tiled with halo time through the join memory the tiled line's slow rise is the accumulated result, not the boundary geometry

Prerequisites

  • python 3.11
  • geopandas 0.14.4, pyogrio 0.8.0, shapely 2.0.4, pandas 2.2.2
  • Boundary data in a format supporting spatial filtering on read — GeoPackage, FlatGeobuf, or GeoParquet with bounding-box columns
  • Case points in the same projected coordinate reference system as the boundaries
  • A known maximum polygon size, which sets the halo width

Step-by-Step Solution

Read only the bounding box of the day’s points; if that is still too large, tile the box, read each tile with a halo, join the points whose home is that tile, and write each tile’s result before reading the next.

# Memory-bounded point-in-polygon against a very large boundary file.
# Pinned: geopandas==0.14.4, pyogrio==0.8.0, shapely==2.0.4, pandas==2.2.2
import logging
import math

import numpy as np
import pandas as pd
import geopandas as gpd
import pyogrio

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


def read_bbox(path: str, bbox: tuple, columns: list[str]) -> gpd.GeoDataFrame:
    """Read only features intersecting bbox, and only the named columns."""
    gdf = pyogrio.read_dataframe(path, bbox=bbox, columns=columns)
    log.info("read %d features in bbox %s", len(gdf), tuple(round(v) for v in bbox))
    return gdf


def tiled_join(points: gpd.GeoDataFrame, boundary_path: str, id_col: str,
               tile_m: float = 20_000.0, halo_m: float = 5_000.0,
               out_path: str | None = None) -> pd.DataFrame:
    """Assign each point to a polygon, one tile of boundary geometry at a time.

    Each point belongs to exactly one tile (its 'home' tile by floor division), so
    it is joined exactly once. Boundaries are read for the tile PLUS a halo, so a
    point near a tile edge still finds a polygon that mostly lies next door.
    """
    points = points.sort_values("record_key").reset_index(drop=True)
    x = points.geometry.x.to_numpy()
    y = points.geometry.y.to_numpy()
    tx = np.floor(x / tile_m).astype(int)
    ty = np.floor(y / tile_m).astype(int)
    points = points.assign(_tx=tx, _ty=ty)
    results, n_tiles, peak = [], 0, 0

    for (i, j), grp in points.groupby(["_tx", "_ty"], sort=True):
        minx, miny = i * tile_m - halo_m, j * tile_m - halo_m
        maxx, maxy = (i + 1) * tile_m + halo_m, (j + 1) * tile_m + halo_m
        polys = read_bbox(boundary_path, (minx, miny, maxx, maxy), [id_col])
        peak = max(peak, len(polys))
        joined = gpd.sjoin(grp.drop(columns=["_tx", "_ty"]), polys[[id_col, "geometry"]],
                           how="left", predicate="within")
        # A point on a shared edge can match two polygons; keep one deterministically.
        joined = (joined.sort_values(["record_key", id_col])
                        .drop_duplicates("record_key", keep="first"))
        res = joined[["record_key", id_col]]
        if out_path:
            res.to_csv(out_path, mode="a", header=(n_tiles == 0), index=False)
        else:
            results.append(res)
        n_tiles += 1
        del polys, joined                   # release the tile before the next read

    log.info("tiled join: %d points over %d tiles, peak %d polygons in memory",
             len(points), n_tiles, peak)
    return pd.concat(results, ignore_index=True) if results else pd.DataFrame()

Assigning each point to a single home tile by floor division is what guarantees it is joined exactly once. The halo is what guarantees it is joined correctly: a point in the corner of its home tile may lie inside a large rural block whose bulk is in the neighbouring tile, and reading the tile’s boundaries without a halo would miss that polygon entirely and leave the point unassigned.

The halo must be at least as large as the biggest polygon’s extent beyond any tile edge, which for census blocks in sparsely populated areas can be several kilometres. Setting it from the maximum polygon extent in the file, rather than guessing, makes the correctness guarantee real.

Why Each Tile Needs a Halo Two adjacent tiles with a large rural polygon whose body lies mostly in the right tile but extends a little into the left one. A case point sits in the left tile, inside the small part of the polygon on that side. Reading only the left tile's exact bounding box would return that polygon because it intersects the box, but a polygon filter using centroids or a clipped read would miss it. Reading the tile plus a halo of the maximum polygon extent guarantees the polygon is present and the point is assigned correctly, while its home tile ensures it is joined only once. Home tile decides once; the halo decides correctly home tile neighbour tile large rural polygon case point dashed: home tile plus halo, the area actually read

Validation & Edge Cases

A halo smaller than the largest polygon. Points near tile edges in sparse areas go unassigned. Compute the maximum polygon extent once from the file’s metadata or a single pass over bounding boxes, and set the halo from it.

Formats without spatial filtering. A shapefile read with a bounding box still scans the whole file on some drivers, which defeats the purpose. Convert large reference layers once to a format with an index — GeoPackage or FlatGeobuf — as part of registering the vintage.

Points on shared edges. A point exactly on the boundary between two polygons matches both. Keep the one with the smallest identifier, consistently, as the code does, and count how many needed the rule.

Unassigned points after the join. Points that match no polygon are either outside the boundary layer’s coverage or in a halo that was too small. Distinguish them by re-reading a larger halo for the unassigned points only.

Accumulated results growing too large. For very large point sets, even the results can exhaust memory. Write each tile’s result to disk as it completes, as the out_path option does, and assemble afterwards if needed.

INFO stream.join read 14208 features in bbox (1540000, 1860000, 1570000, 1890000)
INFO stream.join tiled join: 48900 points over 212 tiles, peak 18114 polygons in memory

A useful reconciliation after any tiled join is to count results: the number of distinct record keys in the output must equal the number of input points. More means a point was joined in two tiles and the home-tile rule failed; fewer means a tile’s result was lost or a point fell outside every tile. Both are cheap to check and both indicate a defect that would otherwise surface only as a slightly wrong count somewhere downstream.

The output reconciliation described above has three possible outcomes, and each points to a different defect:

Count the Output Rows Against the Input Points Three outcomes of comparing distinct record keys in the tiled join output with the number of input points. Equal counts mean every point was joined exactly once. More output rows than input points mean the home-tile rule failed and some points were joined in two tiles. Fewer output rows mean a tile's result was lost or some points fell outside every tile. Each outcome names the part of the procedure to inspect. One count comparison names the defect output = input each point joined once correct output > input points joined in two tiles check the home-tile rule output < input a tile lost, or points outside check writes and tiling run it after every tiled join; it costs one count

Compliance Notes

Record the boundary file and vintage, the format and whether it supports spatial filtering, the tile size and halo width with how the halo was derived, the shared-edge rule and the number of points it resolved, and the count of points left unassigned.

Streaming changes nothing about disclosure — the points are protected data either way and the join runs inside the protected environment — but writing intermediate results to disk tile by tile creates files that did not exist in the in-memory version. Those files contain record keys paired with geography and must be written inside the protected environment, deleted when the job completes, and covered by the same retention rules as any other intermediate.