⚙️ How-To Batch processing Performance & scale Updated

How to Split a Large Layer into Tiles for Batch Processing

Problem statement

One file holds 40% of the country and will not process.

gdf = gpd.read_file("greater_london.gpkg")     # 9.2 GB
MemoryError

Or it loads, and then the operation takes eleven hours while fifteen worker processes sit idle, because the unit of work is the file and there is only one of it.

Splitting the layer into tiles fixes both: each tile fits in memory, and there are hundreds of them so the pool has something to do. But a naive split produces wrong answers — features that straddle a boundary get counted twice, dissolves happen per-tile instead of per-region, and the reassembled output has seams down every tile edge.

Tiling correctly means deciding three things: how big, which tile each feature belongs to, and how much extra to read.

Quick answer

Build a grid, assign each feature to exactly one tile, and read a buffered window when the operation needs neighbours:

import geopandas as gpd
from shapely.geometry import box

def make_grid(bounds, size, crs):
    minx, miny, maxx, maxy = bounds
    cells, y = [], miny
    while y < maxy:
        x = minx
        while x < maxx:
            cells.append(box(x, y, min(x + size, maxx), min(y + size, maxy)))
            x += size
        y += size
    return gpd.GeoDataFrame(
        {"tile_id": range(len(cells))}, geometry=cells, crs=crs
    )

grid = make_grid(gpd.read_file(src, rows=0).total_bounds, 10_000, 27700)
def process_tile(tile, src, buffer_m=0):
    window = tile.geometry.buffer(buffer_m) if buffer_m else tile.geometry
    gdf = gpd.read_file(src, bbox=window.bounds)      # only this window is read
    if gdf.empty:
        return None
    out = expensive_operation(gdf)
    # keep only features whose representative point is in the CORE tile
    return out[out.geometry.representative_point().within(tile.geometry)]
Decision Rule
tile size small enough that one tile fits comfortably in memory; large enough that per-tile overhead is a small share
assignment by representative_point(), so each feature lands in exactly one tile
buffer zero for local operations; larger than the operation's reach for anything needing neighbours 
one 9.2 GB file, 11 h, 1 worker busy
      ↓
612 tiles, ~40 s each, 16 workers → 26 min

Assigning features without double-counting

Panels showing a feature straddling two tiles assigned by intersection versus by representative point.
Intersects puts a straddling feature in both tiles. A representative point puts it in exactly one.

Step-by-step solution

Vertical steps from bounds through grid, assignment, buffered read, processing and reassembly.
Assign by point, read by buffered box, write by core. Three different geometries, on purpose.

1. Get the bounds without reading the data 

import fiona

with fiona.open(src) as handle:
    bounds = handle.bounds          # from the header — instant, even on 9 GB
    count = len(handle)
    crs = handle.crs
print(f"{count:,} features over {bounds}")

gpd.read_file(src, rows=0) also works and returns an empty frame with the CRS. Either way, do not read 9 GB to find out how big it is.

2. Choose the tile size from the data, not from a round number

The goal is roughly equal work per tile, and work follows feature density rather than area. Sample first:

def suggest_tile_size(src, bounds, target_features=50_000, probes=12):
    """Estimate density from a few probe windows, then size the tile to hit a target."""
    import random
    minx, miny, maxx, maxy = bounds
    probe = min(maxx - minx, maxy - miny) / 20
    counts = []
    for _ in range(probes):
        x = random.uniform(minx, maxx - probe)
        y = random.uniform(miny, maxy - probe)
        n = len(gpd.read_file(src, bbox=(x, y, x + probe, y + probe)))
        counts.append(n / (probe * probe))          # features per m²
    density = sorted(counts)[len(counts) // 2]      # median, not mean — outliers
    if density <= 0:
        return probe
    return round((target_features / density) ** 0.5, -2)

size = suggest_tile_size(src, bounds)
print(f"suggested tile size: {size:,.0f} m")     # 10,000 m

Rules of thumb that hold up:

  • 50,000–200,000 features per tile for vector work in GeoPandas.
  • Under a minute per tile so progress is meaningful and a crash loses little.
  • Never smaller than the largest feature — a river polygon spanning three tiles defeats the point.

3. Assign each feature to exactly one tile

This is the step that decides whether counts stay correct.

# WRONG — a feature crossing a boundary appears in both tiles
assigned = gpd.sjoin(gdf, grid, predicate="intersects")
len(assigned) > len(gdf)             # True: duplicates

# RIGHT — a representative point is inside exactly one tile
probes = gdf.copy()
probes["geometry"] = gdf.geometry.representative_point()
assigned = gpd.sjoin(probes, grid, predicate="within")[["tile_id"]]
gdf = gdf.join(assigned)
assert gdf["tile_id"].notna().all()
assert len(gdf) == len(original)

representative_point() rather than centroid: the centroid of a C-shaped polygon can fall outside it, and then the feature lands in the wrong tile or in none. A representative point is guaranteed to be inside the shape.

A point exactly on a tile boundary is within neither tile in some edge cases; catch it:

missing = gdf["tile_id"].isna()
if missing.any():
    fallback = gpd.sjoin_nearest(probes[missing], grid)[["tile_id"]]
    gdf.loc[missing, "tile_id"] = fallback["tile_id"]

4. Read a buffered window when the operation needs neighbours

Whether a buffer is needed depends entirely on whether the operation is local:

# local — a feature's result depends only on that feature. buffer = 0
gdf.to_crs(3857)
gdf.geometry.buffer(50)
gdf.assign(area=gdf.geometry.area)
gdf[gdf["land_use"] == "residential"]

# non-local — needs neighbours. buffer must exceed the operation's reach
gdf.dissolve(by="ward")                      # a ward can span tiles
gpd.sjoin_nearest(gdf, hospitals)            # nearest may be in the next tile
gdf.geometry.buffer(500).union_all()         # reach = 500 m
BUFFER = 500        # must be ≥ the furthest distance the operation looks

def process_tile(tile_geom, src, buffer_m=BUFFER):
    window = tile_geom.buffer(buffer_m)
    gdf = gpd.read_file(src, bbox=window.bounds)          # read wide
    if gdf.empty:
        return None
    out = do_work(gdf)
    core = out.geometry.representative_point().within(tile_geom)
    return out[core]                                       # write narrow

Read wide, write narrow. The buffer supplies context; the core filter guarantees each output feature is produced exactly once, by exactly one tile.

For a nearest-neighbour operation the buffer must exceed the largest plausible distance, which is not always knowable. When it is not, tiling is the wrong approach — use a global spatial index instead.

5. Process the tiles in parallel 

from concurrent.futures import ProcessPoolExecutor, as_completed

def run_tiles(grid, src, out_dir, workers=8):
    out_dir.mkdir(parents=True, exist_ok=True)
    results = []
    with ProcessPoolExecutor(workers) as pool:
        futures = {
            pool.submit(process_and_write, row.geometry, row.tile_id, src, out_dir): row.tile_id
            for row in grid.itertuples()
        }
        for fut in as_completed(futures):
            tile_id = futures[fut]
            try:
                results.append(fut.result())
            except Exception as exc:
                results.append({"tile_id": tile_id, "status": "failed",
                                "error": f"{type(exc).__name__}: {exc}"})
    return results

def process_and_write(tile_geom, tile_id, src, out_dir):
    out = process_tile(tile_geom, src)
    if out is None or out.empty:
        return {"tile_id": tile_id, "status": "empty", "rows": 0}
    path = out_dir / f"tile_{tile_id:05d}.gpkg"
    tmp = path.with_suffix(".tmp.gpkg")
    out.to_file(tmp, driver="GPKG")
    tmp.replace(path)
    return {"tile_id": tile_id, "status": "ok", "rows": len(out)}

One output file per tile. Writing to a shared GeoPackage from several processes corrupts it — SQLite permits one writer.

Empty tiles are normal and should be reported, not treated as failures. A grid over a coastal region is largely sea.

6. Reassemble, and check the arithmetic 

import pandas as pd

def reassemble(out_dir, final, expected_rows=None):
    parts = sorted(out_dir.glob("tile_*.gpkg"))
    frames = [gpd.read_file(p) for p in parts]
    frames = [f for f in frames if not f.empty]
    merged = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs=frames[0].crs)

    if expected_rows is not None:
        assert len(merged) == expected_rows, (
            f"expected {expected_rows:,} rows, got {len(merged):,} "
            f"({len(merged)-expected_rows:+,}) — check the core filter"
        )
    merged.to_file(final, driver="GPKG")
    return merged

That assertion is the whole safety net for tiling. If the reassembled count exceeds the input, the core filter is not excluding buffered features; if it is short, some tile produced nothing when it should have. Either way it is caught before the output is used.

Code examples

Example 1: the complete tiled run 

import fiona, geopandas as gpd, pandas as pd
from pathlib import Path
from shapely.geometry import box
from concurrent.futures import ProcessPoolExecutor, as_completed

def tiled_process(src: Path, out: Path, *, size=10_000, buffer_m=0, workers=8):
    with fiona.open(src) as h:
        bounds, total, crs = h.bounds, len(h), h.crs

    grid = make_grid(bounds, size, crs)
    print(f"{total:,} features → {len(grid)} tiles of {size:,} m")

    tmp_dir = out.parent / f".{out.stem}_tiles"
    results = run_tiles(grid, src, tmp_dir, workers=workers)

    summary = pd.DataFrame(results)
    print(summary["status"].value_counts().to_dict())
    failed = summary[summary["status"] == "failed"]
    if len(failed):
        for r in failed.itertuples():
            print(f"  ✗ tile {r.tile_id}: {r.error}")
        raise SystemExit(f"{len(failed)} tiles failed")

    merged = reassemble(tmp_dir, out, expected_rows=None)
    print(f"wrote {len(merged):,} features to {out}")
    return merged
tiled_process(Path("greater_london.gpkg"), Path("out/london.gpkg"),
              size=10_000, buffer_m=0, workers=16)
# 8,412,993 features → 612 tiles of 10,000 m
# {'ok': 441, 'empty': 171}
# wrote 8,412,993 features to out/london.gpkg

Input and output feature counts matching is the proof that assignment and the core filter are both correct.

Example 2: dissolve across tile boundaries, correctly

Dissolve is the classic operation tiling gets wrong — a ward spanning two tiles dissolves into two pieces.

def tiled_dissolve(src, grid, by="ward", buffer_m=2_000):
    """Dissolve per tile with a buffer, then dissolve the results once more."""
    partials = []
    for row in grid.itertuples():
        window = row.geometry.buffer(buffer_m)
        gdf = gpd.read_file(src, bbox=window.bounds)
        if gdf.empty:
            continue
        partials.append(gdf.dissolve(by=by, aggfunc="sum").reset_index())

    # the second dissolve fuses the pieces that spanned tiles
    combined = gpd.GeoDataFrame(pd.concat(partials, ignore_index=True), crs=partials[0].crs)
    return combined.dissolve(by=by, aggfunc="sum").reset_index()

The two-stage pattern — partial aggregate per tile, then a final aggregate over the partials — is the general answer for any associative operation: sums, unions, counts, min/max. It does not work for medians or percentiles, which cannot be combined from partials.

Note the double-counting risk: with a buffer, the same feature appears in several tiles' partials, so aggfunc="sum" over-counts. Either use buffer 0 and accept split geometries that the second dissolve fuses, or deduplicate on a feature id before the second stage.

Example 3: tiling as a resumable job 

def process_and_write(tile_geom, tile_id, src, out_dir):
    path = out_dir / f"tile_{tile_id:05d}.gpkg"
    done = out_dir / f"tile_{tile_id:05d}.done"
    if done.exists():
        return {"tile_id": tile_id, "status": "skipped"}

    out = process_tile(tile_geom, src)
    if out is not None and not out.empty:
        tmp = path.with_suffix(".tmp.gpkg")
        out.to_file(tmp, driver="GPKG")
        tmp.replace(path)
    done.touch()                       # only after the data is safely on disk
    return {"tile_id": tile_id, "status": "ok", "rows": 0 if out is None else len(out)}

A separate .done marker, written after the data, distinguishes "this tile produced nothing" from "this tile has not run". Without it, an empty tile looks identical to an unprocessed one and gets retried on every run forever.

Explanation

A tile with its buffered read window and its core write area, showing which features are kept.
Three geometries per tile: read the buffer, process everything, write only the core.

Tiling divides a spatial problem the way sharding divides a database one, and it works for the same reason: most operations are local. Reprojecting a parcel needs that parcel and nothing else, so it does not matter which tile it is in, how many tiles there are, or in what order they run.

Every difficulty comes from operations that are not local. A dissolve needs every feature with the same key; a nearest-neighbour search needs whatever is closest, wherever it is. Cutting the data arbitrarily changes those answers, because a tile boundary is not a feature of the world — it is an artefact of how the job was divided.

The buffered read is the standard remedy, and it works when the operation has a bounded reach. A 500 m buffer operation looks at most 500 m away, so reading 500 m past the tile edge supplies every feature it could possibly need. Writing only the core then guarantees each output feature is produced exactly once. Read wide, write narrow.

The remedy fails when the reach is unbounded. A nearest hospital could be 50 km away in a rural tile; no fixed buffer is safe. There the answer is not a bigger buffer but a different structure: hold the small reference layer in memory in every worker, or build a global index once and query it per tile.

Assignment by representative point is the other load-bearing detail. Using intersects puts a straddling feature in every tile it touches, which double-counts in exactly the way that is hardest to notice — the map looks right, the totals are wrong. A representative point is inside exactly one tile, and representative_point() rather than centroid because the centroid of a concave shape can fall outside it entirely.

Finally, tiling changes the unit of work, and with it every operational property: peak memory drops to one tile, parallelism becomes useful, resume granularity becomes one cell, and progress becomes honest because the tiles are comparable.

Edge cases or notes

  • read_file(bbox=…) uses the spatial index when the format has one. GeoPackage and FlatGeobuf do; shapefile needs a .qix sidecar or it scans.
  • A feature larger than a tile is read by every tile it touches and assigned to one. That is correct, but if such features are common the tiles are too small.
  • centroid on a concave polygon falls outside it. Always representative_point() for assignment.
  • Empty tiles are normal, especially over sea or unpopulated area. Report them; do not treat them as failures.
  • Concurrent writes to one GeoPackage corrupt it. One file per tile, merged afterwards.
  • pd.concat of frames with different CRS silently produces a frame with no CRS. Assert the CRS after reassembly.
  • Tile the projected CRS, not lat/lon — a degree grid has wildly varying cell areas by latitude.
  • Medians and percentiles cannot be combined from partials. Those need the whole dataset, or an approximation like t-digest.
  • gdal_retile.py and ogr2ogr -clipsrc do the same job from the command line for rasters and vectors respectively.

FAQ

What tile size should I use?

Whatever gives 50,000–200,000 features and under a minute per tile. Sample the density with a few probe windows rather than guessing from area.

How do I stop features being counted twice?

Assign each feature to exactly one tile by its representative_point(), and after processing keep only features whose representative point is inside the core tile.

When do I need a buffer?

When the operation looks at neighbouring features — dissolve, nearest, buffer-and-union. The buffer must be at least as large as the operation's reach.

Can I tile a dissolve?

Yes, in two stages: dissolve per tile, then dissolve the partial results. That works for any associative aggregate. It does not work for medians or percentiles.

Why representative_point instead of centroid?

The centroid of a C-shaped or ring-shaped polygon can fall outside the polygon, so the feature is assigned to the wrong tile or to none.

Should each tile write its own file?

Yes. Concurrent writes to one GeoPackage corrupt it. Merge the per-tile outputs at the end.

How do I know the tiling was correct?

Compare the reassembled feature count against the input count. More means the core filter is leaking; fewer means a tile silently produced nothing.