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
Bounds without reading, tile size, one tile per feature
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
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.
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
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.qixsidecar 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.
centroidon a concave polygon falls outside it. Alwaysrepresentative_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.concatof 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.pyandogr2ogr -clipsrcdo the same job from the command line for rasters and vectors respectively.
Internal links
- Choosing the unit of work in a batch job β why tile rather than file
- How to speed up batch GIS jobs with parallel processing β what tiling makes possible
- How to process a very large GeoPackage in chunks β the non-spatial alternative
- Fixing memory errors in GeoPandas when working with large files β the symptom that sends people here
- How to build a resumable batch GIS job in Python β the
.donemarker pattern - Batch GIS job slows down and runs out of memory β per-tile memory as the fix
- How to create a fishnet grid in Python with GeoPandas β the grid, in more detail
- Spatial indexes explained: R-trees and why spatial joins are fast β what makes a bbox read cheap
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.