How to Split a Large Layer into Tiles for Batch Processing

On this page

Problem statement

One file holds a whole region and will not process: reading it raises MemoryError, or it loads and the operation runs for 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
    )

src = "parcels.gpkg"                                    # Manchester's 121,064 land parcels
grid = make_grid(gpd.read_file(src, rows=0).total_bounds, 2_000, 27700)
print(len(grid), "tiles of 2 km")
84 tiles of 2 km
def expensive_operation(gdf):
    """Stands in for the real work: here, a 5 m buffer and its area."""
    out = gdf.assign(geometry=gdf.geometry.buffer(5))
    return out.assign(area=out.geometry.area)


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)]


parts = [process_tile(row, src) for _, row in grid.iloc[36:44].iterrows()]
print("parcels from tiles 36-43:", [0 if part is None else len(part) for part in parts])
parcels from tiles 36-43: [5011, 4827, 8209, 3867, 0, 0, 5, 4149]
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

Measured on the Manchester parcels in Example 1: the whole file took 3.0 s in one process, and 84 tiles on 8 workers took 2.3 s. That is a small gain, because the file fits in memory and a 5 m buffer is cheap. What tiling guaranteed was that all 121,064 parcels came out exactly once β€” and the same code runs unchanged on a layer that does not fit.

Assigning features without double-counting

A feature straddling two tiles assigned by intersection versus by representative point. a feature crossing a tile boundarysjoin intersectsβœ•appears in both tilesβœ•processed twiceβœ•output count exceeds inputβœ•every total inflatedβœ•the map still looks rightrepresentative_point withinβœ“one point, one tileβœ“processed onceβœ“output count equals inputβœ“totals correctβœ“assignment is deterministiccentroid can fall outside a concave shape β€” use representative_point spatialworkflow.io
Intersects puts a straddling feature in both tiles. A representative point puts it in exactly one.

Bounds without reading, tile size, one tile per feature

Steps for splitting a layer into tiles and reassembling the result. three geometries per tile, on purpose1bounds from the headerdo not read 9 GB2size the grid from density50k–200k features3assign by representative pointexactly one tile4read buffered, write coreif the op needs neighbours5reassemble and check the countin must equal outOne output file per tile β€” concurrent GeoPackage writes corrupt it spatialworkflow.io
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}")
121,064 features over (378638.25, 381998.739, 390812.47, 405638.237)

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
    random.seed(0)
    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 for 50,000 features: {size:,.0f} m")
print(f"for 2,000 features: {suggest_tile_size(src, bounds, target_features=2_000):,.0f} m")
suggested tile size for 50,000 features: 20,100 m
for 2,000 features: 4,000 m

Manchester's whole parcel layer is under three times the 50,000-feature target, so the suggestion is a tile bigger than the city; the smaller target is what this page uses to get enough tiles to share out.

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.

gdf = gpd.read_file(src)
original = gdf.copy()

# WRONG β€” a feature crossing a boundary appears in both tiles
assigned = gpd.sjoin(gdf, grid, predicate="intersects")
print(f"intersects: {len(assigned):,} rows for {len(gdf):,} parcels")

# 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)
print(f"representative point: {gdf['tile_id'].notna().sum():,} assigned, {gdf['tile_id'].isna().sum()} on a boundary")
assert len(gdf) == len(original)
intersects: 124,545 rows for 121,064 parcels
representative point: 121,064 assigned, 0 on a boundary

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"]
print(gdf["tile_id"].notna().all())
True
Rows produced by assigning 120,989 Manchester parcels to 2 km tiles by intersection and by representative point. 120,989 parcels, 84 tilesthe parcels120,989sjoin intersects124,467 Β· 3,478 duplicatesrepresentative point120,989 Β· each exactly once0130kA 2.9% excess looks like nothing on a map and is wrong in every total spatialworkflow.io
Measured in step 3.

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

These, with the helpers above, are the module Example 1 runs:

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.

Tiling Manchester parcels, with a dissolve across tiles

The data used on this page

The layer is HM Land Registry's INSPIRE Index Polygons for Manchester β€” 121,064 freehold parcels in British National Grid β€” saved as parcels.gpkg. For the dissolve, each parcel also carries the ONS ward (December 2024) its representative point falls in, and a count column n, in parcels_wards.gpkg.

import os

import geopandas as gpd
import requests

WARDS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
         "Wards_December_2024_Boundaries_UK_BGC/FeatureServer/0/query"
         "?where=LAD24NM%3D%27Manchester%27&outFields=WD24CD,WD24NM&outSR=27700&f=geojson")
if not os.path.exists("parcels.gpkg"):
    if not os.path.exists("Manchester_City_Council.zip"):
        with requests.Session() as session:                # the service sets a cookie, then redirects
            open("Manchester_City_Council.zip", "wb").write(session.get(
                "https://use-land-property-data.service.gov.uk/datasets/inspire/download/"
                "Manchester_City_Council.zip", timeout=600).content)
    parcels = gpd.read_file("/vsizip/Manchester_City_Council.zip/Land_Registry_Cadastral_Parcels.gml")
    parcels = parcels[["INSPIREID", "geometry"]].rename(columns={"INSPIREID": "inspire_id"})
    parcels.set_crs(27700, allow_override=True).to_file("parcels.gpkg")
if not os.path.exists("parcels_wards.gpkg"):
    parcels = gpd.read_file("parcels.gpkg")
    wards = gpd.read_file(WARDS)
    points = parcels.assign(geometry=parcels.representative_point())
    parcels["ward"] = gpd.sjoin(points, wards[["WD24NM", "geometry"]], how="left",
                                predicate="within")["WD24NM"].values
    parcels[parcels["ward"].notna()].assign(n=1).to_file("parcels_wards.gpkg")
for path in ("parcels.gpkg", "parcels_wards.gpkg"):
    layer = gpd.read_file(path)
    print(f"{path:20} {len(layer):,} features" + (f", {layer['ward'].nunique()} wards" if "ward" in layer else ""))
parcels.gpkg         121,064 features
parcels_wards.gpkg   120,369 features, 32 wards

Example 1: the complete tiled run

The functions from the steps, as one module β€” worker processes have to be able to import what they run:

%%writefile tiling.py
from concurrent.futures import ProcessPoolExecutor, as_completed
from pathlib import Path

import fiona
import geopandas as gpd
import pandas as pd
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)


def expensive_operation(gdf):
    out = gdf.assign(geometry=gdf.geometry.buffer(5))
    return out.assign(area=out.geometry.area)


def process_tile(tile_geom, src, buffer_m=0):
    window = tile_geom.buffer(buffer_m) if buffer_m else tile_geom
    gdf = gpd.read_file(src, bbox=window.bounds)
    if gdf.empty:
        return None
    out = expensive_operation(gdf)
    return out[out.geometry.representative_point().within(tile_geom)]


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)}


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):
            try:
                results.append(fut.result())
            except Exception as exc:
                results.append({"tile_id": futures[fut], "status": "failed",
                                "error": f"{type(exc).__name__}: {exc}"})
    return results


def reassemble(out_dir, final, expected_rows=None):
    frames = [gpd.read_file(p) for p in sorted(out_dir.glob("tile_*.gpkg"))]
    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


def tiled_process(src: Path, out: Path, *, size=10_000, 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):
        raise SystemExit(f"{len(failed)} tiles failed: {failed['error'].tolist()[:3]}")
    merged = reassemble(tmp_dir, out, expected_rows=total)
    print(f"wrote {len(merged):,} features to {out}")
    return merged
import time
from pathlib import Path

import tiling

started = time.perf_counter()
tiling.expensive_operation(gpd.read_file("parcels.gpkg")).to_file("whole_file.gpkg")
print(f"whole file, one process: {time.perf_counter() - started:.1f} s")

started = time.perf_counter()
tiling.tiled_process(Path("parcels.gpkg"), Path("out/parcels_buffered.gpkg"), size=2_000, workers=8)
print(f"84 tiles, 8 workers: {time.perf_counter() - started:.1f} s")
whole file, one process: 3.0 s
121,064 features β†’ 84 tiles of 2,000 m
{'ok': 49, 'empty': 35}
wrote 121,064 features to out/parcels_buffered.gpkg
84 tiles, 8 workers: 2.3 s

Input and output feature counts matching is the proof that assignment and the core filter are both correct β€” reassemble asserts it rather than printing it.

Example 2: dissolve across tile boundaries, correctly

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

import pandas as pd


def tiled_dissolve(src, grid, by="ward", buffer_m=0, core_only=True):
    """Dissolve per tile, then dissolve the partial results once more."""
    partials = []
    for row in grid.itertuples():
        window = row.geometry.buffer(buffer_m) if buffer_m else row.geometry
        gdf = gpd.read_file(src, bbox=window.bounds)
        if core_only:                                   # each feature in exactly one tile
            gdf = gdf[gdf.geometry.representative_point().within(row.geometry)]
        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()
layer = gpd.read_file("parcels_wards.gpkg")
direct = layer.dissolve(by="ward", aggfunc="sum")
grid = tiling.make_grid(layer.total_bounds, 2_000, 27700)
for label, kwargs in (("buffer 2,000 m, no core filter", {"buffer_m": 2_000, "core_only": False}),
                      ("buffer 0, no core filter", {"buffer_m": 0, "core_only": False}),
                      ("buffer 0, core filter", {"buffer_m": 0, "core_only": True})):
    result = tiled_dissolve("parcels_wards.gpkg", grid, **kwargs).set_index("ward")
    area_gap = (result.geometry.area - direct.geometry.area).abs().max()
    print(f"{label:32} parcels counted {int(result['n'].sum()):9,} of {len(layer):,}; "
          f"{len(result)} wards; largest area difference {area_gap:.2f} mΒ²")
buffer 2,000 m, no core filter   parcels counted 1,054,608 of 120,369; 32 wards; largest area difference 0.01 mΒ²
buffer 0, no core filter         parcels counted   123,904 of 120,369; 32 wards; largest area difference 0.00 mΒ²
buffer 0, core filter            parcels counted   120,369 of 120,369; 32 wards; largest area difference 0.00 mΒ²

The union geometry was right in all three β€” dissolving twice fuses the pieces β€” and the sums were not. A bbox read returns every feature that touches the tile's box, so a parcel straddling a tile edge is read by every tile it touches: without a core filter even a zero buffer counted some parcels twice, and the 2 km buffer counted each parcel many times over. Assigning each feature to one tile before the first dissolve is what makes the sum a sum.

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.

Parcel counts from a tiled two-stage dissolve of Manchester parcels by ward, with and without a buffer and a core filter. 120,294 parcels dissolved into 32 wards, in 84 tilesper-tile readparcels countedgeometry2 km buffer, no core filter1,053,948rightno buffer, no core filter123,827rightno buffer, core filter120,294rightThe union comes out right every time; the sum is right only when each parcel belongs to one tile spatialworkflow.io
Measured in Example 2.

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)}

Run over the first ten tiles twice, with the second run standing in for a restart after a crash:

out_dir = Path("resumable")
out_dir.mkdir(exist_ok=True)
grid = tiling.make_grid(gpd.read_file("parcels.gpkg", rows=0).total_bounds, 2_000, 27700)
process_tile = lambda geom, src: tiling.process_tile(geom, src)
for attempt in (1, 2):
    statuses = [process_and_write(row.geometry, row.tile_id, "parcels.gpkg", out_dir)["status"]
                for row in grid.head(10).itertuples()]
    print(f"run {attempt}: {pd.Series(statuses).value_counts().to_dict()}")
run 1: {'ok': 10}
run 2: {'skipped': 10}

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.

Most operations are local

The three geometries involved in processing one tile. read wide, process everything, write narrowgeometry usedwhyread windowtile buffered by 500 mneighbours the operation needsprocesseverything readso edge features are correctwrite coreβœ“the tile itselfβœ“each feature written onceThe buffer must exceed the operation’s reach, or edges are wrong spatialworkflow.io
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.

bbox reads, oversized features and empty tiles

  • 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 a tile that fits in memory and runs in under a minute β€” 50,000–200,000 features is a common guide. 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 β€” after assigning each feature to one tile, or the sums over-count. Measured on Manchester's parcels, skipping the assignment counted 2.9% of parcels twice with no buffer, and far more with one. It does not work for medians or percentiles.

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.