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.12 seconds
gdf["area"] = gdf.apply(lambda r: r.geometry.area, axis=1)   # same rows, 20.6 seconds

167 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 a minute when it could take a tenth of a second. Or a one-million-row GeoPackage that memory_usage says is 97 MB takes more than a gigabyte of RAM. Or a script that ran in a minute last month now takes an hour, and nothing changed except the data got bigger.

The timings on this page were measured on an AMD Ryzen 9 7900 with Python 3.14, pandas 3.0.6, GeoPandas 1.1.4, Shapely 2.1.2 and pyogrio, on synthetic data the code builds itself, so you can run it and get your own.

"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

# 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.

Python loops, no index, parsing, memory pressure

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.006 s
[g.area for g in geometry]           0.321 s
gdf.geometry.apply(area)             0.267 s
gdf.apply(axis=1)                    0.969 s
iterrows                             2.700 s

Read those carefully, because the gaps are not where people expect:

• A plain list comprehension is 53× slower than the vectorised call, but 3× faster than apply(axis=1). Iterating geometries is not the disaster.
• gdf.apply(..., axis=1) is 160× slower than the vectorised call, because it builds a Series for every row — copying every column, not just the geometry.
• iterrows is 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 was 6.8× faster than iterrows on the same frame (0.37 s against 2.48 s; see why apply is slow).

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.

import shapely


def make_points(n, seed=0):
    """n random points in a 100 km square."""
    rng = np.random.default_rng(seed)
    xy = rng.uniform(0, 100_000, (n, 2))
    return gpd.GeoDataFrame(geometry=gpd.points_from_xy(xy[:, 0], xy[:, 1]), crs=27700)


def make_polygons(n, seed=1):
    """n Voronoi cells tiling the same square, like zones or wards."""
    rng = np.random.default_rng(seed)
    square = shapely.box(0, 0, 100_000, 100_000)
    cells = shapely.get_parts(shapely.voronoi_polygons(shapely.multipoints(rng.uniform(0, 100_000, (n, 2))),
                                                       extend_to=square))
    cells = shapely.intersection(cells, square)
    return gpd.GeoDataFrame({"name": [f"zone {i}" for i in range(len(cells))]}, geometry=cells, crs=27700)


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                          1.018 s
vectorised contains                  0.032 s
sjoin                                0.022 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:

fresh = make_polygons(1_000)          # sjoin above already built the index on `polys`
t0 = time.perf_counter(); fresh.sindex; print(f"build {time.perf_counter()-t0:.6f} s")
t0 = time.perf_counter(); fresh.sindex; print(f"cached {time.perf_counter()-t0:.6f} s")

build 0.000238 s
cached 0.000031 s

A quarter of a millisecond to build over 1,000 polygons, and 11 ms over the 50,000 parcels in Example 1, so 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 a large frame repeatedly. See how to use the spatial index directly.

3. I/O and parsing — slow before any work happens

To see this on your own machine, write one set of test parcels — a million irregular 24-vertex polygons with ten attribute columns — in each format:

rng = np.random.default_rng(7)
n = 1_000_000
cx, cy = rng.uniform(380_000, 420_000, n), rng.uniform(390_000, 430_000, n)
angle = np.linspace(0, 2 * np.pi, 24, endpoint=False)
radius = rng.uniform(8, 25, (n, 1)) * rng.uniform(0.85, 1.15, (n, 24))
rings = np.stack([cx[:, None] + radius * np.cos(angle), cy[:, None] + radius * np.sin(angle)], axis=-1)
test = gpd.GeoDataFrame({"id": np.arange(n), "class": rng.choice(["residential", "commercial", "industrial", "green"], n),
                         **{f"attr_{i}": rng.random(n) for i in range(8)}},
                        geometry=shapely.polygons(rings), crs=27700)
test.to_file("parcels.shp")
test.to_file("parcels.gpkg")
test.to_file("parcels.geojson")
test.to_parquet("parcels.parquet")
del test

for path, label in [("parcels.shp", "Shapefile"), ("parcels.gpkg", "GeoPackage"), ("parcels.geojson", "GeoJSON")]:
    timed(f"read {label}", lambda p=path: gpd.read_file(p))
timed("read GeoParquet", lambda: gpd.read_parquet("parcels.parquet"))    # read_file is for GDAL formats

read Shapefile                       2.867 s
read GeoPackage                      2.479 s
read GeoJSON                        19.806 s
read GeoParquet                      0.997 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.

Reading less helps with any format. Here a bounding box covering a fifth of the parcels cut the GeoPackage read from 2.5 s to 1.4 s, and two columns instead of eleven cut its memory from 97 to 33 MB:

reads = {
    "full read": {},
    "two columns": {"columns": ["id", "class"]},               # only the columns you use
    "bbox read": {"bbox": (380_000, 395_000, 400_000, 410_000)},  # only the area, via the file's spatial index
    "1,000 rows": {"rows": 1000},                             # only some rows, for a quick look
}
for label, kwargs in reads.items():
    t0 = time.perf_counter()
    part = gpd.read_file("parcels.gpkg", **kwargs)
    print(f"{label:<12} {time.perf_counter() - t0:7.3f} s  {len(part):>9,} rows  "
          f"{part.memory_usage(deep=True).sum() / 1e6:6.0f} MB")
del part

full read      2.461 s  1,000,000 rows      97 MB
two columns    2.181 s  1,000,000 rows      33 MB
bbox read      1.421 s    187,368 rows      18 MB
1,000 rows     0.046 s      1,000 rows       0 MB

Also check the engine. pyogrio is the default in GeoPandas 1.0, and on the million parcels it was twelve times faster than fiona:

print(gpd.options.io_engine)         # None → pyogrio if installed
for engine in ("pyogrio", "fiona"):
    t0 = time.perf_counter()
    gpd.read_file("parcels.gpkg", engine=engine)
    print(f"{engine:<8} {time.perf_counter() - t0:7.2f} s")

None
pyogrio     2.82 s
fiona      35.08 s

4. Memory pressure — the one that looks like everything else

The symptom is slowness, not an error: once the process swaps, every operation slows by orders of magnitude while the CPU sits idle.

import subprocess
import sys

# a fresh interpreter: one that has already read files reuses that memory and under-reports
measure = """
import os, psutil, geopandas as gpd
rss = lambda: psutil.Process(os.getpid()).memory_info().rss / 1e6
before = rss()
gdf = gpd.read_file("parcels.gpkg")
print(f"{len(gdf):,} rows, {rss() - before:,.0f} MB, {gdf.memory_usage(deep=True).sum() / 1e6:,.0f} MB reported")
"""
print(subprocess.run([sys.executable, "-c", measure], capture_output=True, text=True, check=True).stdout, end="")

1,000,000 rows, 1,068 MB, 97 MB reported

Run fresh, because the same lines inside this page's session, after the reads above, reported 52 MB: the allocator reused memory the process already held. The gap between 97 MB reported and 1,068 MB actual is the point: memory_usage(deep=True) counts the geometry column as one 8-byte pointer per row, and the GEOS geometries behind them are not counted. Real memory was eleven times the reported figure, and peak during parsing is higher still.

Every operation that returns a new frame adds to it, up to a full second copy:

gdf = gpd.read_file("parcels.gpkg")
gdf2 = gdf[gdf["class"] == "residential"]        # a copy of the matching rows
gdf3 = gdf.copy()                                # a full duplicate
gdf4 = gdf.to_crs(27700)                         # new geometries — a full duplicate

The fix is not holding 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).strip_dirs().sort_stats("cumulative").print_stats(top)
    print(s.getvalue())
    return result

def assign_zones(points, zones):
    return [zones.index[zones.contains(pt)] for pt in points.geometry]


def my_pipeline(points, zones):
    parcels = gpd.read_file("parcels.gpkg")
    return parcels, assign_zones(points, zones)


profile(lambda: my_pipeline(points.head(20_000), polys))

         7085360 function calls (6925090 primitive calls) in 5.008 seconds

   Ordered by: cumulative time
   List reduced from 450 to 12 due to restriction <12>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.000    0.000    5.008    5.008 cell_279.py:20(<lambda>)
        1    0.000    0.000    5.007    5.007 cell_279.py:15(my_pipeline)
        1    0.013    0.013    2.767    2.767 cell_279.py:11(assign_zones)
        1    0.000    0.000    2.240    2.240 file.py:197(_read_file)
        1    0.058    0.058    2.240    2.240 file.py:482(_read_file_pyogrio)
        1    0.000    0.000    2.182    2.182 geopandas.py:138(read_dataframe)
    20000    0.006    0.000    1.995    0.000 base.py:2583(contains)
    20000    0.025    0.000    1.989    0.000 base.py:83(_binary_op)
        1    1.726    1.726    1.726    1.726 raw.py:39(read)
    20000    0.034    0.000    1.518    0.000 base.py:27(_delegate_binary_method)
    20000    0.006    0.000    0.885    0.000 array.py:803(contains)
    20000    0.011    0.000    0.879    0.000 array.py:782(_binary_method)

The time splits almost evenly: 2.2 s reading the GeoPackage and 2.8 s in assign_zones, most of it 20,000 calls to contains. They need different fixes. GeoParquet would cut the read (step 3); the loop is bottlenecks 1 and 2 together, sjoin would remove nearly all of it, and it is the part that grows — 49 s on all 184,204 points in Example 2. Changing the file format alone would have saved under 2 s. Full technique in how to profile a slow Python GIS script.

A benchmark that names your bottleneck

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):
    """Run it in a fresh interpreter: memory deltas mean little in a process that has already freed memory."""
    import psutil
    proc = psutil.Process(os.getpid())
    rss = lambda: proc.memory_info().rss / 1e6

    findings = []

    # ── I/O ──────────────────────────────────────────────────────────────
    gpd.read_file(path, rows=10)          # load the drivers first, so they are not counted as data
    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)")
    import pyogrio
    total_rows = pyogrio.read_info(path)["features"]
    if total_rows > len(gdf):
        print(f"          full file has {total_rows:,} rows — expect roughly "
              f"{read_mb * total_rows / len(gdf):,.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     0.14 s       87 MB  (2.8 ms per 1k rows)
area      vectorised 0.003 s   apply(axis=1) would be 0.3 s (83x)
sindex    built in 0.011 s for 50,000 rows
memory    5 MB reported, 87 MB actual (18.0x)
          full file has 1,000,000 rows — expect roughly 1,744 MB in memory
vertices  mean 25  p95 25  max 25

  → row-wise apply would cost 83x the vectorised call — never use apply(axis=1)

That was run in a fresh interpreter on the test parcels. The projection from a sample is the useful part: a fraction of a second of reading predicts about 1.7 GB for the whole file, which is a decision you want before starting rather than after twenty minutes. Treat it as an upper bound — the full read took 1,068 MB, because part of what any read allocates is fixed cost that a small sample over-weights.

The vertex check catches a bottleneck people rarely look for. GEOS predicate cost scales with vertex count, so one geometry with hundreds of thousands of vertices can dominate a join where every other feature has a few dozen. The test parcels are uniform — 25 coordinates each — so it found nothing here.

Example 2: the same task, four ways

import time
import geopandas as gpd
import numpy as np
from shapely.geometry import Point

points = make_points(184_204)                # from step 2: random points in a 100 km square
zones = make_polygons(215).rename(columns={"name": "ward_name"})   # 215 Voronoi "wards" tiling it

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             56.02 s   184,204 matched
apply + contains        49.21 s   184,204 matched
manual sindex            0.67 s   184,204 matched
sjoin                    0.11 s   184,204 matched

Identical results, and 509× between the extremes. The two jumps are the two bottlenecks separately: adding the index takes 49 s to 0.67 s (that is bottleneck 2), and moving from a Python loop to the vectorised sjoin takes 0.67 s to 0.11 s (bottleneck 1).

Note that approach_3 is already fast enough for most work. The index is the change that matters — 73× — and vectorising the remaining loop is a further 6×.

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: 33 → 33 MB
  · simplified: 25,000,000 → 16,152,800 vertices (35% fewer)
  · spatial index built for 1,000,000 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.

A 1 m tolerance removed 35% of the test parcels' vertices, whose edges wobble by up to 15% of their radius; how much it removes from real boundaries depends on how finely they were digitised, so measure it on yours. Every vertex removed makes subsequent predicates cheaper. It is lossy, so it belongs in display and screening work rather than anywhere the exact boundary matters.

Two stacks: pandas for attributes, GEOS for shapes

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 was 53× faster than a list comprehension and 160× 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 an allocation and a copy of every column's value, for every row. On pandas 3.0 the fixed part dominates — a 41-column frame took only 24% longer than a 2-column one — so the cost is mostly the Series itself. itertuples avoids it and was 6.8× 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. Geometry lives outside what pandas counts, so a million-row frame reported as 97 MB took 1,068 MB, 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.

Per-row Series, cached indexes and vertex counts

• apply(axis=1) builds a Series per row. The cost scales with column count, so dropping unused columns helps even before vectorising.
• itertuples was 6.8× faster than iterrows and is the right escape hatch when a loop is genuinely needed.
• sindex is 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 badly: it counts one 8-byte pointer per geometry. On the test parcels, real memory was eleven times the reported figure.
• pyogrio was twelve times faster than fiona for a bulk read and is the GeoPandas 1.0 default.
• bbox= on read_file uses 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.values is a GeometryArray that 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 53× slower than the vectorised call but 3× faster than apply(axis=1). What is catastrophic is iterating rows — apply(axis=1) and iterrows.

How much memory will my file need?

Far more than memory_usage(deep=True) reports — eleven times more on the test parcels — because it counts a pointer per geometry and none of the geometry itself. Read a 50,000-row sample and extrapolate before loading the whole thing.

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.