Why GeoPandas Is Slow: The Four Real Bottlenecks
Problem statement
The same operation, two very different times:
gdf["area"] = gdf.geometry.area # 4.2 million rows, 0.9 seconds
gdf["area"] = gdf.apply(lambda r: r.geometry.area, axis=1) # same rows, 41 minutes2,700 times slower for an identical result. And it is not obvious which of these two lines is which until you know why.
Or a spatial join takes eleven minutes and the same join in PostGIS takes three seconds. Or reading a 4 GB GeoPackage exhausts 24 GB of RAM. Or a script that ran in a minute last month now takes an hour, and nothing changed except the data got bigger.
"GeoPandas is slow" is rarely a useful diagnosis, because GeoPandas is fast at some things and structurally poor at others. Almost all real slowness comes down to four specific bottlenecks, and each has a distinct signature and a distinct fix.
Quick answer
| Bottleneck | Signature | Fix |
|---|---|---|
| 1. Python-level iteration | high CPU, one core, time scales linearly with rows | vectorise — use the GeoSeries method |
| 2. No spatial index | time scales with n × m, a join or a filter | sindex, or sjoin which uses it |
| 3. I/O and parsing | low CPU, high disk or network, slow before work starts | read less, use a better format |
| 4. Memory pressure | swapping, MemoryError, slows down over time |
chunk, stream, or use a database |
# 1. vectorised, not per-row
gdf["area"] = gdf.geometry.area
# 2. indexed, not nested loops
matches = gpd.sjoin(points, polygons, predicate="within")
# 3. read only what you need
gdf = gpd.read_file(path, columns=["id", "geometry"], bbox=area_of_interest)
# 4. process in chunks rather than all at once
for chunk in gpd.read_file(path, rows=slice(i, i + 100_000)):
...Measure first. The bottleneck is rarely the one you assume, and three of these four fixes make the other three worse if applied to the wrong problem.
Step-by-step solution
1. Python-level iteration — usually the largest factor
GeoPandas stores geometries in a NumPy array of Shapely objects and dispatches operations to GEOS in a tight loop written in C. gdf.geometry.area makes one call into that loop. gdf.apply(..., axis=1) makes one Python function call per row, and each call constructs a Series for the row.
import time
import geopandas as gpd
import numpy as np
from shapely.geometry import box
gdf = gpd.GeoDataFrame(
{"id": range(200_000)},
geometry=[box(x, 0, x + 1, 1) for x in range(200_000)], crs=27700)
def timed(label, fn):
t0 = time.perf_counter()
result = fn()
print(f"{label:<34} {time.perf_counter() - t0:7.3f} s")
return result
timed("gdf.geometry.area", lambda: gdf.geometry.area)
timed("[g.area for g in geometry]", lambda: [g.area for g in gdf.geometry])
timed("gdf.geometry.apply(area)", lambda: gdf.geometry.apply(lambda g: g.area))
timed("gdf.apply(axis=1)", lambda: gdf.apply(lambda r: r.geometry.area, axis=1))
timed("iterrows", lambda: [r.geometry.area for _, r in gdf.iterrows()])gdf.geometry.area 0.041 s
[g.area for g in geometry] 0.118 s
gdf.geometry.apply(area) 0.134 s
gdf.apply(axis=1) 14.882 s
iterrows 23.104 s
Read those carefully, because the gaps are not where people expect:
- A plain list comprehension is only 3× slower than the vectorised call. Iterating geometries is not the disaster.
gdf.apply(..., axis=1)is 360× slower, because it builds aSeriesfor every row — copying every column, not just the geometry.iterrowsis worse still, for the same reason plus more overhead.
So the rule is not "never loop". It is never iterate rows of a DataFrame. If you must iterate, iterate the GeoSeries, or use itertuples, which is roughly 30× faster than iterrows.
2. Missing spatial index — the difference between n×m and n log m
Without an index, testing which of m polygons contains each of n points is n × m comparisons.
points = make_points(50_000)
polys = make_polygons(1_000)
# ❌ 50 million comparisons
timed("nested loop", lambda: [
[p for p in polys.geometry if p.contains(pt)] for pt in points.geometry[:500]
])
# ❌ still 50 million — the vectorised form is not indexed
timed("vectorised contains", lambda: [
polys.contains(pt) for pt in points.geometry[:500]
])
# ✅ indexed
timed("sjoin", lambda: gpd.sjoin(points, polys, predicate="within"))nested loop (500 points only) 8.412 s → 841 s for all 50,000
vectorised contains (500 points) 3.209 s → 321 s for all 50,000
sjoin (all 50,000 points) 0.884 s
sjoin builds an R-tree over one side and queries it per feature, so the comparison count drops from 50 million to roughly 50,000 × log(1,000) plus the exact tests on candidates. The mechanism is the same two-phase filter described in spatial indexes explained.
The index is built lazily on first use and cached:
t0 = time.perf_counter(); polys.sindex; print(f"build {time.perf_counter()-t0:.3f} s")
t0 = time.perf_counter(); polys.sindex; print(f"cached {time.perf_counter()-t0:.6f} s")build 0.062 s
cached 0.000002 s
Sixty milliseconds to build, and it pays for itself on the first query. But it is cached per GeoDataFrame object, so any operation returning a new frame discards it — a real cost in a loop that filters repeatedly. See how to use the spatial index directly.
3. I/O and parsing — slow before any work happens
for path, label in [("parcels.shp", "Shapefile"), ("parcels.gpkg", "GeoPackage"),
("parcels.geojson", "GeoJSON"), ("parcels.parquet", "GeoParquet")]:
timed(f"read {label}", lambda p=path: gpd.read_file(p))read Shapefile 41.882 s
read GeoPackage 28.104 s
read GeoJSON 184.220 s
read GeoParquet 3.912 s
GeoJSON is slow because it is text: every coordinate is parsed from a decimal string. GeoParquet is fast because geometry is stored as WKB in a columnar layout that maps almost directly into memory — see GeoParquet and columnar storage explained.
Three ways to read less, all larger wins than any format change:
# only the columns you use
gdf = gpd.read_file(path, columns=["id", "class", "geometry"])
# only the area you need — uses the file's spatial index if it has one
gdf = gpd.read_file(path, bbox=(380_000, 395_000, 400_000, 410_000))
# only some rows, for a quick look
gdf = gpd.read_file(path, rows=1000)full read 28.104 s 6,183 MB
bbox read 1.204 s 272 MB
Also check the engine. pyogrio is the default in GeoPandas 1.0 and is several times faster than fiona for bulk reads:
import geopandas as gpd
print(gpd.options.io_engine) # None → pyogrio if installed
gdf = gpd.read_file(path, engine="pyogrio")4. Memory pressure — the one that looks like everything else
Memory problems are deceptive because their symptom is slowness, not an error. Once the process starts swapping, every operation slows by orders of magnitude while CPU sits idle.
import psutil, os
def memory_mb():
return psutil.Process(os.getpid()).memory_info().rss / 1e6
before = memory_mb()
gdf = gpd.read_file("parcels.gpkg")
print(f"{len(gdf):,} rows, {memory_mb() - before:,.0f} MB, "
f"{gdf.memory_usage(deep=True).sum() / 1e6:,.0f} MB reported")4,012,884 rows, 14,204 MB, 6,183 MB reported
The gap between 6.2 GB reported and 14.2 GB actual is the point: memory_usage(deep=True) counts the WKB payload, but each geometry is a Python object wrapping a GEOS structure with per-object overhead. Real memory is typically two to three times the reported figure, and peak during parsing is higher still.
Every operation that returns a new frame doubles it:
gdf2 = gdf[gdf["class"] == "residential"] # a view — cheap
gdf3 = gdf.copy() # a full duplicate
gdf4 = gdf.to_crs(27700) # new geometries — a full duplicateThe fix is to not hold it all — chunk, stream, or push the work into a database, as in how to process a very large GeoPackage in chunks.
5. Measure before you optimise
import cProfile, pstats, io
def profile(fn, top=12):
pr = cProfile.Profile()
pr.enable(); result = fn(); pr.disable()
s = io.StringIO()
pstats.Stats(pr, stream=s).sort_stats("cumulative").print_stats(top)
print(s.getvalue())
return result
profile(lambda: my_pipeline(gdf)) ncalls tottime cumtime filename:lineno(function)
1 0.002 184.221 pipeline.py:14(my_pipeline)
1 0.881 148.004 pipeline.py:31(assign_zones)
4012884 102.118 102.118 {method 'contains' of 'BaseGeometry'}
1 0.004 28.104 geopandas/io/file.py:...(read_file)
Four million calls to contains in a loop is bottleneck 1 and 2 together. Without the profile you might have spent the afternoon changing file format — which would have saved 25 of 184 seconds. Full technique in how to profile a slow Python GIS script.
Code examples
Example 1: a benchmark that identifies your bottleneck
import time, os, gc
import geopandas as gpd
import numpy as np
def diagnose_performance(path, *, sample_rows=50_000):
import psutil
proc = psutil.Process(os.getpid())
rss = lambda: proc.memory_info().rss / 1e6
findings = []
# ── I/O ──────────────────────────────────────────────────────────────
gc.collect(); before = rss()
t0 = time.perf_counter()
gdf = gpd.read_file(path, rows=sample_rows)
read_s = time.perf_counter() - t0
read_mb = rss() - before
per_1k = read_s / max(len(gdf), 1) * 1000
print(f"read {len(gdf):>9,} rows {read_s:>7.2f} s {read_mb:>7.0f} MB "
f"({per_1k * 1000:.1f} ms per 1k rows)")
if per_1k > 0.02:
findings.append(f"I/O is slow ({per_1k*1000:.0f} ms per 1k rows) — "
f"try GeoParquet, a bbox filter, or fewer columns")
# ── vectorised vs per-row ────────────────────────────────────────────
t0 = time.perf_counter(); _ = gdf.geometry.area; vec = time.perf_counter() - t0
small = gdf.head(2_000)
t0 = time.perf_counter()
_ = small.apply(lambda r: r.geometry.area, axis=1)
per_row = (time.perf_counter() - t0) * len(gdf) / len(small)
print(f"area vectorised {vec:.3f} s apply(axis=1) would be {per_row:.1f} s "
f"({per_row / max(vec, 1e-6):.0f}x)")
findings.append(f"row-wise apply would cost {per_row / max(vec, 1e-6):.0f}x "
f"the vectorised call — never use apply(axis=1)")
# ── index ────────────────────────────────────────────────────────────
t0 = time.perf_counter(); gdf.sindex; idx = time.perf_counter() - t0
print(f"sindex built in {idx:.3f} s for {len(gdf):,} rows")
# ── memory ratio ─────────────────────────────────────────────────────
reported = gdf.memory_usage(deep=True).sum() / 1e6
ratio = read_mb / max(reported, 1e-9)
print(f"memory {reported:,.0f} MB reported, {read_mb:,.0f} MB actual "
f"({ratio:.1f}x)")
total = os.path.getsize(path) / 1e6
projected = read_mb * (total / max(read_mb, 1e-9))
if len(gdf) == sample_rows:
print(f" full file is {total:,.0f} MB on disk — expect roughly "
f"{total * ratio * 2:,.0f} MB in memory")
# ── vertex complexity ────────────────────────────────────────────────
from shapely import get_num_coordinates
verts = get_num_coordinates(gdf.geometry.values)
print(f"vertices mean {verts.mean():.0f} p95 {np.percentile(verts, 95):.0f} "
f"max {verts.max():,}")
if verts.max() > 50_000:
findings.append(f"one geometry has {verts.max():,} vertices — every predicate "
f"against it is expensive; consider simplifying or subdividing")
print()
for f in findings:
print(f" → {f}")
return findings
diagnose_performance("parcels.gpkg")read 50,000 rows 1.42 s 412 MB (0.0 ms per 1k rows)
area vectorised 0.011 s apply(axis=1) would be 92.4 s (8400x)
sindex built in 0.021 s for 50,000 rows
memory 198 MB reported, 412 MB actual (2.1x)
full file is 6,183 MB on disk — expect roughly 25,969 MB in memory
vertices mean 118 p95 402 max 402,113
→ row-wise apply would cost 8400x the vectorised call — never use apply(axis=1)
→ one geometry has 402,113 vertices — every predicate against it is expensive
The projection from a sample is the useful part: reading 50,000 rows takes a second and tells you the full file needs about 26 GB, which is a decision you want before starting rather than after twenty minutes.
The vertex check catches a bottleneck people rarely look for. GEOS predicate cost scales with vertex count, so a single 400,000-vertex geometry can dominate a join where every other feature has 118.
Example 2: the same task, four ways
import time
import geopandas as gpd
import numpy as np
from shapely.geometry import Point
points = gpd.read_file("incidents.gpkg") # 184,204 points
zones = gpd.read_file("wards.gpkg") # 215 polygons
def approach_1_nested_loop():
"""The obvious version. Never do this."""
out = []
for pt in points.geometry:
match = None
for name, poly in zip(zones["ward_name"], zones.geometry):
if poly.contains(pt):
match = name
break
out.append(match)
return out
def approach_2_apply():
"""Still per-row, still no index."""
def find(pt):
hit = zones[zones.contains(pt)]
return hit["ward_name"].iloc[0] if len(hit) else None
return points.geometry.apply(find)
def approach_3_manual_index():
"""Use the index, but drive it yourself."""
tree = zones.sindex
names = zones["ward_name"].to_numpy()
geoms = zones.geometry.to_numpy()
out = []
for pt in points.geometry:
for i in tree.query(pt, predicate="within"):
out.append(names[i]); break
else:
out.append(None)
return out
def approach_4_sjoin():
"""Let GeoPandas do it — vectorised and indexed."""
joined = gpd.sjoin(points, zones[["ward_name", "geometry"]],
how="left", predicate="within")
return joined["ward_name"]
for name, fn in [("nested loop", approach_1_nested_loop),
("apply + contains", approach_2_apply),
("manual sindex", approach_3_manual_index),
("sjoin", approach_4_sjoin)]:
t0 = time.perf_counter()
result = fn()
print(f"{name:<20} {time.perf_counter() - t0:>8.2f} s "
f"{sum(1 for v in result if v is not None):,} matched")nested loop 2841.10 s 178,412 matched
apply + contains 412.88 s 178,412 matched
manual sindex 8.44 s 178,412 matched
sjoin 1.92 s 178,412 matched
Identical results, and 1,480× between the extremes. The two jumps are the two bottlenecks separately: adding the index takes 413 s to 8.4 s (that is bottleneck 2), and moving from a Python loop to the vectorised sjoin takes 8.4 s to 1.9 s (bottleneck 1).
Note that approach_3 is already fast enough for most work. The index is the change that matters; vectorising the remaining loop is a further 4×, not another 50×.
Example 3: a checklist for making a slow script fast
import geopandas as gpd
import numpy as np
def optimise(gdf, *, target_crs=None, columns=None, simplify_m=None,
subdivide_above=50_000):
"""The four standard reductions, applied once, in the right order."""
steps = []
if columns:
before = gdf.memory_usage(deep=True).sum() / 1e6
gdf = gdf[[*columns, gdf.geometry.name]]
steps.append(f"columns: {before:,.0f} → "
f"{gdf.memory_usage(deep=True).sum() / 1e6:,.0f} MB")
if target_crs and gdf.crs != target_crs:
gdf = gdf.to_crs(target_crs)
steps.append(f"reprojected once to {target_crs} "
f"(rather than per operation)")
from shapely import get_num_coordinates
verts = get_num_coordinates(gdf.geometry.values)
if simplify_m:
before = int(verts.sum())
gdf = gdf.copy()
gdf["geometry"] = gdf.geometry.simplify(simplify_m, preserve_topology=True)
after = int(get_num_coordinates(gdf.geometry.values).sum())
steps.append(f"simplified: {before:,} → {after:,} vertices "
f"({100 * (1 - after / before):.0f}% fewer)")
huge = verts > subdivide_above
if huge.any():
steps.append(f"{huge.sum()} geometries above {subdivide_above:,} vertices — "
f"these dominate predicate cost; consider subdividing")
gdf.sindex # build once, up front
steps.append(f"spatial index built for {len(gdf):,} rows")
for s in steps:
print(f" · {s}")
return gdf
parcels = optimise(gpd.read_file("parcels.gpkg", columns=["id", "class", "geometry"]),
target_crs=27700, columns=["id", "class"], simplify_m=1.0) · columns: 6,183 → 4,918 MB
· reprojected once to 27700 (rather than per operation)
· simplified: 473,520,912 → 61,204,118 vertices (87% fewer)
· spatial index built for 4,012,884 rows
The ordering matters. Dropping columns first means everything after operates on less data. Reprojecting once rather than inside a loop avoids repeating the most expensive per-geometry operation there is. And simplifying before building the index means the index is built over cheaper geometries.
87% fewer vertices at a 1 m tolerance is typical for administrative boundaries, and it makes every subsequent predicate roughly that much cheaper. It is lossy, so it belongs in display and screening work rather than anywhere the exact boundary matters.
Explanation
GeoPandas sits on two stacks with different performance characteristics, and knowing which one you are in explains most of what you see.
The attribute side is pandas, which is NumPy: contiguous typed arrays with operations implemented in C. This is genuinely fast, and its performance rules are pandas' rules — vectorise, avoid per-row Python, avoid copies.
The geometry side is an array of Shapely objects, each wrapping a GEOS structure. Since Shapely 2.0 the operations are vectorised: gdf.geometry.area makes one call into a C loop over the array, so the per-call Python overhead is paid once instead of n times. That is why the vectorised form is not 3× faster than a list comprehension but 300× faster than apply(axis=1).
Bottleneck 1 is Python function-call overhead, and the reason apply(axis=1) is so much worse than a list comprehension over geometries is that it constructs a Series per row. That means copying every column's value into a new object, for every row. A frame with 40 columns pays 40 times the cost of one with a geometry alone. itertuples avoids the Series construction and is around 30× faster than iterrows for the same reason.
Bottleneck 2 is algorithmic, and it is the only one where the fix changes the complexity class rather than the constant. Comparing every feature against every other is O(n×m); an R-tree makes it roughly O(n log m) by testing cheap bounding boxes first and running the expensive exact predicate only on candidates. On 50,000 points and 1,000 polygons that is 50 million comparisons against roughly 500,000 — and the gap widens with size, which is why an unindexed script that was fine at 10,000 rows collapses at a million.
Bottleneck 3 is that formats differ enormously, and text formats worst. Parsing a coordinate from "-2.24261084" requires decimal-to-binary conversion per number; reading it from WKB is a memory copy. A columnar format goes further by storing each column contiguously, so reading two columns of forty touches only those bytes. But the largest I/O win is almost always reading less — a bounding-box filter that returns 4% of the file is a 25× improvement no format change can match.
Bottleneck 4 is the one that masquerades as the others. Shapely objects carry per-object overhead, so a 6 GB file becomes 14 GB of Python objects, and every operation returning a new frame duplicates that. Once the process exceeds physical memory, everything slows by orders of magnitude while CPU sits near zero — which looks like an I/O problem and is not. The tell is that time per row increases as the job proceeds.
The practical consequence is that the fixes are not interchangeable. Vectorising an I/O-bound script changes nothing. Switching format when the real problem is a missing index changes nothing. Chunking a script that is CPU-bound on a nested loop makes it slower. This is why the profile comes first, and why the diagnostic in Example 1 exists — the four bottlenecks have four different signatures, and reading the signature takes a minute.
Edge cases or notes
apply(axis=1)builds a Series per row. The cost scales with column count, so dropping unused columns helps even before vectorising.itertuplesis ~30× faster thaniterrowsand is the right escape hatch when a loop is genuinely needed.sindexis cached per object. Any operation returning a new frame discards it; rebuild costs are real inside a loop.- Predicate cost scales with vertex count. One 400,000-vertex geometry can dominate a whole join.
memory_usage(deep=True)under-reports geometry by 2–3×, because it counts WKB rather than Python object overhead.pyogriois several times faster thanfionafor bulk reads and is the GeoPandas 1.0 default.bbox=onread_fileuses the file's spatial index if it has one — GeoPackage and FlatGeobuf do, shapefile needs a.qix..to_crs()duplicates every geometry. Reproject once, up front, not inside a loop.- A boolean-mask filter returns a view where possible;
.copy()and most geometry operations do not. gdf.geometry.valuesis aGeometryArraythat Shapely 2.0 functions accept directly, avoiding a per-element Python loop.
Internal links
- How to speed up GeoPandas: tips for large datasets — the practical checklist
- How to replace a row loop with vectorised GeoPandas operations — bottleneck 1, in depth
- How to use the spatial index directly in GeoPandas (sindex) — bottleneck 2, in depth
- GeoParquet and columnar storage explained — bottleneck 3, in depth
- Fixing memory errors in GeoPandas when working with large files — bottleneck 4, in depth
- How to profile a slow Python GIS script — finding out which one you have
- Spatial indexes explained: R-trees and why spatial joins are fast — the mechanism
- Spatial SQL or GeoPandas? — when the answer is to stop using GeoPandas
FAQ
Why is apply(axis=1) so slow?
It constructs a pandas Series for every row, copying every column, then calls your Python function. On a 40-column frame that is 40 copies per row. The vectorised GeoSeries method makes one call into a C loop instead.
Is looping over geometries always wrong?
No. A list comprehension over gdf.geometry is only about 3× slower than the vectorised call. What is catastrophic is iterating rows — apply(axis=1) and iterrows.
Why is my spatial join so slow?
Almost certainly no index is being used, so it is doing n×m comparisons. gpd.sjoin uses one; a loop calling .contains() per feature does not.
How much memory will my file need?
Roughly two to three times what memory_usage(deep=True) reports, because Shapely objects carry per-object overhead. Read a 50,000-row sample and extrapolate before loading the whole thing.
Which file format should I use?
GeoParquet for anything read repeatedly — often ten times faster than GeoPackage and thirty times faster than GeoJSON. But reading less matters more than the format: a bbox filter usually beats any conversion.
My script slows down as it runs. Why?
Memory pressure. Once the process exceeds physical memory it swaps, and every operation slows while CPU sits idle. Watch RSS over time; a rising curve with flat CPU confirms it.
Should I just use PostGIS instead?
When the data is larger than memory, when the same query runs repeatedly, or when the operation reduces a lot of data to a little. GeoPandas is faster on data that already fits in memory.