How to Profile a Slow Python GIS Script and Find the Real Bottleneck

Problem statement

The script takes forty minutes. You believe you know why, and you are usually wrong.

The common experience: you spend an afternoon converting everything to GeoParquet, and the run goes from 41 minutes to 39 — because reading was 2 minutes and a nested loop was 37. Or you parallelise across eight cores and it gets slower, because the job was memory-bound and eight processes each loaded a copy.

Intuition is unreliable here for a specific reason. GIS code mixes four cost centres — Python-level iteration, GEOS predicate work, file I/O, and memory pressure — and they produce similar-looking symptoms while needing opposite fixes. Chunking a CPU-bound loop makes it slower. Vectorising an I/O-bound script changes nothing.

Profiling takes about ten minutes and settles it.

Quick answer

Start coarse, then narrow:

import cProfile, pstats, io

pr = cProfile.Profile()
pr.enable()
main()
pr.disable()

s = io.StringIO()
pstats.Stats(pr, stream=s).sort_stats("cumulative").print_stats(20)
print(s.getvalue())

Each rung is more detailed and more expensive. Most problems are found on the first two.

Tool	Answers	Overhead
`time.perf_counter` around phases	which phase	none
`cProfile`	which function	2–3×
`line_profiler` (`@profile`)	which line	10–100×
`memray` / `tracemalloc`	which allocation	2–10×
`py-spy`	a running process, no code change	~0

pip install line_profiler memray py-spy snakeviz

Do the coarse timing first. It costs nothing, needs no tooling, and usually points straight at the phase to profile properly.

Step-by-step solution

1. Time the phases before profiling anything

import time
from contextlib import contextmanager

TIMINGS = {}

@contextmanager
def phase(name):
    t0 = time.perf_counter()
    try:
        yield
    finally:
        TIMINGS[name] = TIMINGS.get(name, 0.0) + time.perf_counter() - t0

with phase("read"):
    parcels = gpd.read_file("parcels.gpkg")
with phase("clean"):
    parcels = clean(parcels)
with phase("join"):
    parcels = assign_zones(parcels, zones)
with phase("aggregate"):
    summary = summarise(parcels)
with phase("write"):
    summary.to_file("out.gpkg", driver="GPKG")

total = sum(TIMINGS.values())
for name, secs in sorted(TIMINGS.items(), key=lambda kv: -kv[1]):
    print(f"{name:<12} {secs:>8.1f} s  {100 * secs / total:>5.1f}%")

join           2214.8 s   88.4%
read             28.1 s    1.1%
clean           184.2 s    7.4%
write            62.4 s    2.5%
aggregate        15.9 s    0.6%

Eighteen lines, and the answer is already clear: the join is 88% of the run and the file format is 1%. Everything after this is about the join.

Amdahl's law is the reason this matters. Making read ten times faster improves the total run by 1%. Making join twice as fast improves it by 44%.

2. Profile the phase that dominates

import cProfile, pstats

def profile_call(fn, *args, sort="cumulative", top=20, dump=None, **kwargs):
    pr = cProfile.Profile()
    pr.enable()
    result = fn(*args, **kwargs)
    pr.disable()
    stats = pstats.Stats(pr)
    if dump:
        stats.dump_stats(dump)
    stats.sort_stats(sort).print_stats(top)
    return result

profile_call(assign_zones, parcels, zones, dump="join.prof")

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1    0.412    0.412 2214.802 2214.802 pipeline.py:88(assign_zones)
  4012884  118.204    0.000 2189.441    0.001 pipeline.py:102(<listcomp>)
  4012884 1841.220    0.000 1841.220    0.000 {method 'contains' of 'BaseGeometry'}
  4012884  204.118    0.000  204.118    0.000 geopandas/geoseries.py:...(__getitem__)

Two columns matter and they answer different questions:

tottime — time in the function itself, excluding calls it makes. The 1,841 seconds in contains is real work.
cumtime — time including everything it called. Useful for finding the caller responsible.

Sort by cumulative to find where the time goes; sort by tottime to find what is slow.

Four million calls to contains says the loop is testing every parcel against something without an index. That is the diagnosis; the rest is applying a spatial join.

Visualise a large profile rather than reading columns:

snakeviz join.prof            # opens an interactive flame graph in a browser

3. Narrow to the line

cProfile reports per function, which is not enough when one function is fifty lines:

# pip install line_profiler
# decorate the function, then run: kernprof -lv script.py

@profile                       # line_profiler injects this name
def clean(gdf):
    gdf = gdf[gdf.geometry.notna()]
    gdf = gdf[~gdf.geometry.is_empty]
    gdf["valid"] = gdf.geometry.is_valid
    gdf.loc[~gdf["valid"], "geometry"] = gdf.loc[~gdf["valid"], "geometry"].make_valid()
    gdf["area"] = gdf.geometry.area
    gdf = gdf[gdf["area"] > 1.0]
    gdf["centroid"] = gdf.geometry.representative_point()
    return gdf

Line #  Hits    Time   Per Hit   % Time  Line Contents
    12     1    412.0     412.0     0.2  gdf = gdf[gdf.geometry.notna()]
    13     1    388.0     388.0     0.2  gdf = gdf[~gdf.geometry.is_empty]
    14     1  18204.0   18204.0    10.0  gdf["valid"] = gdf.geometry.is_valid
    15     1 141882.0  141882.0    77.6  gdf.loc[~gdf["valid"], "geometry"] = ...
    16     1   9204.0    9204.0     5.0  gdf["area"] = gdf.geometry.area
    17     1    204.0     204.0     0.1  gdf = gdf[gdf["area"] > 1.0]
    18     1  12881.0   12881.0     7.0  gdf["centroid"] = ...

Line 15 is 78% of the function. make_valid on 8,412 invalid geometries out of 4 million is inherently expensive — but the .loc assignment on a boolean mask is doing more work than it needs to. Rewriting it to operate on the subset and reassign is a straightforward win.

line_profiler adds 10–100× overhead, so use it on one function with a small input, never on a whole run.

4. Profile memory separately — it does not show in time profiles

Panels contrasting the CPU and memory signatures of a CPU-bound job and a memory-bound one. — A memory-bound job looks slow with idle CPU. A time profile attributes the slowness to whatever line happened to touch memory.

# pip install memray
memray run -o profile.bin script.py
memray flamegraph profile.bin        # writes an HTML flame graph
memray summary profile.bin

📦 Total allocated: 41.2 GiB
📈 Peak memory: 18.4 GiB

  Location                                              Peak      Allocations
  pipeline.py:44 gpd.read_file(...)                    14.2 GiB          8,412
  pipeline.py:102 parcels.to_crs(27700)                 6.1 GiB      4,012,884
  pipeline.py:118 gpd.sjoin(...)                        2.8 GiB        184,204

to_crs allocating 6.1 GiB across 4 million allocations is the finding: it constructs a new Shapely object per geometry. Called once that is a cost; called inside a loop it is the bottleneck.

A lighter option with no dependency:

import tracemalloc

tracemalloc.start()
result = expensive_step(gdf)
current, peak = tracemalloc.get_traced_memory()
print(f"current {current/1e6:,.0f} MB, peak {peak/1e6:,.0f} MB")
for stat in tracemalloc.take_snapshot().statistics("lineno")[:5]:
    print(f"  {stat}")
tracemalloc.stop()

tracemalloc sees only Python-level allocations, so it under-reports GEOS and GDAL memory. memray intercepts the allocator and sees everything, which for GIS work is the difference between a useful number and a misleading one.

5. Profile a running process when you cannot restart it

A job that has been running for six hours should not be restarted to profile it:

# pip install py-spy
py-spy dump --pid 48122                 # a one-off stack snapshot
py-spy top --pid 48122                  # a live, top-like view
py-spy record -o profile.svg --pid 48122 --duration 60

Thread 48122 (active): "MainThread"
    contains (shapely/geometry/base.py:748)
    <listcomp> (pipeline.py:102)
    assign_zones (pipeline.py:88)
    main (pipeline.py:21)

One command, no code changes, no restart, and the stack names the line. py-spy reads the process's memory from outside, so overhead is effectively zero and it works on a production job.

This is the first thing to reach for when a scheduled job is mysteriously slow — see why a GIS pipeline fails silently for the related diagnostic problem.

6. Confirm the fix with a measurement, not a feeling

import time, statistics

def benchmark(fn, *args, runs=5, warmup=1, **kwargs):
    for _ in range(warmup):
        fn(*args, **kwargs)
    times = []
    for _ in range(runs):
        t0 = time.perf_counter()
        result = fn(*args, **kwargs)
        times.append(time.perf_counter() - t0)
    print(f"{fn.__name__:<24} min {min(times):.3f}s  "
          f"median {statistics.median(times):.3f}s  "
          f"max {max(times):.3f}s")
    return result

before = benchmark(assign_zones_loop, parcels, zones, runs=3)
after = benchmark(assign_zones_sjoin, parcels, zones, runs=3)
assert (before["zone"].fillna("") == after["zone"].fillna("")).all(), \
    "the optimised version returns different results"

The assertion is the part people skip. A rewrite that is faster and returns different answers is not an optimisation, and spatial rewrites change answers easily — a predicate swapped from within to intersects, a boundary case handled differently, a join that now duplicates rows.

Report min rather than mean: the minimum is the closest to the true cost, with the least contamination from other processes.

Code examples

Example 1: a profiling harness for a pipeline

import time, io, cProfile, pstats, os, gc
from contextlib import contextmanager
from dataclasses import dataclass, field

@dataclass
class Profiler:
    phases: dict = field(default_factory=dict)
    rows: dict = field(default_factory=dict)
    peaks: dict = field(default_factory=dict)

    def _rss(self):
        try:
            import psutil
            return psutil.Process(os.getpid()).memory_info().rss / 1e6
        except ImportError:
            return float("nan")

    @contextmanager
    def phase(self, name, rows=None):
        gc.collect()
        before = self._rss()
        t0 = time.perf_counter()
        try:
            yield
        finally:
            elapsed = time.perf_counter() - t0
            after = self._rss()
            self.phases[name] = self.phases.get(name, 0.0) + elapsed
            self.peaks[name] = max(self.peaks.get(name, 0.0), after - before)
            if rows is not None:
                self.rows[name] = rows

    def report(self):
        total = sum(self.phases.values()) or 1.0
        print(f"{'phase':<18}{'seconds':>10}{'share':>8}{'ΔRSS':>10}{'per 1k rows':>14}")
        for name, secs in sorted(self.phases.items(), key=lambda kv: -kv[1]):
            n = self.rows.get(name)
            per = f"{secs / n * 1000 * 1000:>11.1f} ms" if n else f"{'':>14}"
            print(f"{name:<18}{secs:>9.1f}s{100*secs/total:>7.1f}%"
                  f"{self.peaks[name]:>9.0f}M{per}")
        print(f"{'TOTAL':<18}{total:>9.1f}s")

        biggest = max(self.phases, key=self.phases.get)
        share = 100 * self.phases[biggest] / total
        print(f"\n  → '{biggest}' is {share:.0f}% of the run. "
              f"Halving it saves {share/2:.0f}% overall; "
              f"halving everything else saves {(100-share)/2:.0f}%.")
        return biggest

prof = Profiler()

with prof.phase("read", rows=None):
    parcels = gpd.read_file("parcels.gpkg")
with prof.phase("clean", rows=len(parcels)):
    parcels = clean(parcels)
with prof.phase("reproject", rows=len(parcels)):
    parcels = parcels.to_crs(27700)
with prof.phase("join", rows=len(parcels)):
    parcels = assign_zones(parcels, zones)
with prof.phase("write", rows=len(parcels)):
    parcels.to_parquet("out.parquet")

worst = prof.report()

phase                seconds   share      ΔRSS   per 1k rows
join                 2214.8s   88.4%     2841M      551.9 ms
clean                 184.2s    7.4%     6104M       45.9 ms
write                  62.4s    2.5%      412M       15.6 ms
read                   28.1s    1.1%    14204M
reproject              15.9s    0.6%     6188M        4.0 ms
TOTAL                2505.4s

  → 'join' is 88% of the run. Halving it saves 44% overall; halving everything
    else saves 6%.

Two columns beyond the obvious earn their place. ΔRSS shows read allocating 14 GB while taking only 1% of the time — a memory finding a time profile would never surface, and the reason the job needs a large machine. And per 1k rows normalises the cost so phases with different inputs are comparable, and makes it easy to project what happens when the data doubles.

The Amdahl sentence at the end is there because it is the calculation people skip. It converts "the join is slow" into "optimising anything else is capped at 6%".

Example 2: finding which geometries are expensive

Sometimes the bottleneck is not a line of code but a handful of features:

import time
import numpy as np
import geopandas as gpd
import shapely

def profile_by_geometry(gdf, operation, *, sample=2_000, seed=0):
    """Time an operation per feature and correlate cost with complexity."""
    rng = np.random.default_rng(seed)
    idx = rng.choice(len(gdf), size=min(sample, len(gdf)), replace=False)
    subset = gdf.iloc[idx]

    verts = shapely.get_num_coordinates(subset.geometry.values)
    times = np.empty(len(subset))
    for i, geom in enumerate(subset.geometry):
        t0 = time.perf_counter()
        operation(geom)
        times[i] = time.perf_counter() - t0

    order = np.argsort(-times)
    print(f"sampled {len(subset):,} of {len(gdf):,} features")
    print(f"  total sampled time {times.sum():.2f}s")
    print(f"  slowest 1% account for {100 * times[order[:len(order)//100 or 1]].sum() / times.sum():.0f}%")
    print(f"  correlation(cost, vertices) = {np.corrcoef(times, verts)[0,1]:.2f}")
    print("\n  slowest features:")
    for i in order[:5]:
        print(f"    index {subset.index[i]}  {verts[i]:>8,} vertices  "
              f"{times[i]*1000:>8.1f} ms")

    all_verts = shapely.get_num_coordinates(gdf.geometry.values)
    projected = times.mean() * len(gdf) * (all_verts.mean() / verts.mean())
    print(f"\n  projected full-layer cost: {projected:.0f}s")
    return times, verts

profile_by_geometry(parcels, lambda g: g.buffer(50))

sampled 2,000 of 4,012,884 features
  total sampled time 1.84s
  slowest 1% account for 61%
  correlation(cost, vertices) = 0.94

  slowest features:
    index 412881   402,113 vertices    884.2 ms
    index  88104   288,044 vertices    602.1 ms
    index   4102     1,204 vertices      3.1 ms

  projected full-layer cost: 3691s

A correlation of 0.94 confirms cost scales with vertex count, and "the slowest 1% account for 61%" is the actionable finding: simplifying or subdividing a few hundred features would cut the operation roughly in half.

This kind of skew is common and invisible to cProfile, which reports one aggregate number for four million calls to buffer. Only per-feature timing reveals that most of it is a handful of coastline polygons.

Example 3: automated regression detection

Profiling once fixes today's problem. Recording the numbers catches tomorrow's:

import json, time, platform, sys
from pathlib import Path

BASELINE = Path("perf_baseline.json")

def timed_run(name, fn, *args, **kwargs):
    t0 = time.perf_counter()
    result = fn(*args, **kwargs)
    return name, time.perf_counter() - t0, result

def check_performance(results, *, tolerance=0.25, update=False):
    """Compare against a committed baseline and fail on a regression."""
    baseline = json.loads(BASELINE.read_text()) if BASELINE.exists() else {}
    current = {name: round(secs, 3) for name, secs, _ in results}

    regressions, improvements = [], []
    for name, secs in current.items():
        was = baseline.get(name)
        if was is None:
            print(f"  new   {name:<24} {secs:>8.2f}s")
            continue
        change = (secs - was) / was
        flag = "SLOWER" if change > tolerance else "faster" if change < -tolerance else "ok"
        print(f"  {flag:<6}{name:<24} {was:>8.2f}s → {secs:>8.2f}s "
              f"({change:+.0%})")
        if change > tolerance:
            regressions.append((name, was, secs))
        elif change < -tolerance:
            improvements.append((name, was, secs))

    if update or not BASELINE.exists():
        BASELINE.write_text(json.dumps(
            {**current, "_recorded": {"python": sys.version.split()[0],
                                      "platform": platform.platform()}}, indent=2))
        print(f"\n  baseline written to {BASELINE}")
    elif regressions:
        raise SystemExit(
            f"\n{len(regressions)} performance regression(s) beyond {tolerance:.0%}:\n" +
            "\n".join(f"  {n}: {w:.2f}s → {c:.2f}s" for n, w, c in regressions))
    return regressions

results = [
    timed_run("read_parcels", gpd.read_file, "tests/data/parcels_sample.gpkg"),
    timed_run("clean", clean, parcels),
    timed_run("join_zones", assign_zones, parcels, zones),
]
check_performance(results)

  ok    read_parcels               1.42s →     1.38s (-3%)
  ok    clean                      4.18s →     4.02s (-4%)
  SLOWER join_zones                 1.92s →     8.41s (+338%)

1 performance regression(s) beyond 25%:
  join_zones: 1.92s → 8.41s

A 338% regression in the join is the kind of change that would otherwise be noticed weeks later, after several more commits, when nobody can say which one caused it. Here it fails on the pull request that introduced it.

The 25% tolerance is deliberately loose, because CI runners are noisy and a tight threshold produces false alarms that get ignored. Recording the platform in the baseline is what stops a comparison between a laptop and a runner being treated as a regression.

Explanation

Halving the 88% phase saves 44% of the run. Halving everything else saves 6%.

Profiling is worth doing because intuition about performance is systematically wrong, and in GIS it is wrong in a particular way: the four cost centres produce overlapping symptoms.

A CPU-bound loop and a memory-bound job both look like "the script is slow". They need opposite fixes — one wants vectorisation, the other wants chunking, and each fix makes the other problem worse. A time profile cannot distinguish them, because a memory-bound job spends its time in whatever line happens to touch memory when the swapping starts, which attributes the cost to an innocent line.

Amdahl's law is the reason to start coarse. If one phase is 88% of the run, no amount of work on the other 12% can improve the total by more than 12%. Phase timing takes eighteen lines and tells you where optimisation can pay, before any tooling is installed. Skipping it is how people spend an afternoon on file formats worth 1%.

tottime and cumtime answer different questions, and confusing them is the commonest misreading of a profile. cumtime includes everything a function called, so the top of a cumulative sort is always main — useful for finding the caller responsible for a cost. tottime excludes callees, so it finds the function actually burning cycles. In GIS profiles the tottime leader is usually a GEOS method, and the useful question is which caller invokes it four million times.

Sampling profilers and tracing profilers make different trade-offs. cProfile traces every call, which is accurate on call counts and adds 2–3× overhead — enough to distort the ratio between cheap and expensive calls. py-spy samples the stack at intervals, which has almost no overhead and gives statistically accurate proportions without exact counts. For a running production job, sampling is the only option; for finding a hot function during development, tracing is more informative.

Memory needs its own tool because Python's does not see most of it. tracemalloc hooks the Python allocator, so it reports Python objects. But a GeoDataFrame's real weight is GEOS structures and GDAL buffers allocated by C code, which tracemalloc never sees. memray intercepts at the allocator level and sees everything — which is why its numbers for GIS code are often several times larger, and closer to the truth.

Finally, the discipline that makes profiling pay off twice: record the numbers. A profile answers "why is this slow now". A committed baseline answers "which commit made it slow", which is a much harder question to answer after the fact. It costs a JSON file and a CI step, and it turns performance from something you investigate into something you notice.

Edge cases or notes

cProfile adds 2–3× overhead and inflates the apparent cost of many small calls relative to a few large ones.
line_profiler adds 10–100×. Use it on one function with a small input.
tracemalloc misses C allocations, which is most of a GeoDataFrame. Use memray for real numbers.
py-spy needs no code change and no restart — the right tool for a job already running.
Report min of several runs, not the mean. The minimum is least contaminated by other processes.
Warm up before timing. The first call pays for imports, lazy index builds and disk cache misses.
Assert the results match after an optimisation. Spatial rewrites change answers easily.
gc.collect() before a memory measurement, or you measure garbage that has not been collected yet.
snakeviz profile.prof turns a dump into an interactive flame graph, which is far easier than reading columns.
Profiling in a container may report host CPU counts, which distorts any parallelism you are measuring.

Internal links

Why GeoPandas is slow: the four real bottlenecks — what you are looking for
How to replace a row loop with vectorised GeoPandas operations — the usual fix
How to use the spatial index directly in GeoPandas (sindex) — the other usual fix
GeoPandas .apply() takes hours: how to fix it — the commonest profile result
Fixing memory errors in GeoPandas when working with large files — when the profile says memory
How to speed up GeoPandas: tips for large datasets — the checklist afterwards
GIS tests run out of memory or time out in CI — profiling in a constrained environment
How to add a progress bar and ETA to a long GIS batch job — knowing it is slow before profiling why

FAQ

Where should I start?

Time the phases with perf_counter. It costs nothing and usually identifies the one phase worth profiling properly. Anything else is premature.

What is the difference between `tottime` and `cumtime`?

tottime is time in the function itself; cumtime includes everything it called. Sort by cumtime to find the responsible caller, by tottime to find the function burning cycles.

Why does my time profile not explain the slowness?

Probably memory. Once the process swaps, everything slows while CPU sits idle, and the time is attributed to whichever line happened to touch memory. Profile with memray instead.

How do I profile a job that is already running?

py-spy dump --pid <pid> gives a stack snapshot with almost no overhead and no restart. It is the right tool for a production job.

Does profiling change the timings?

Yes. cProfile adds 2–3× and biases against many small calls; line_profiler adds far more. Use profiles for relative cost, and untraced timing for absolute numbers.

How do I know my optimisation actually worked?

Benchmark before and after with several runs, report the minimum, and assert that the results are identical. A faster function returning different answers is not an optimisation.

Can I catch performance regressions automatically?

Record phase timings to a committed JSON baseline and fail CI when one regresses beyond a tolerance. Use a loose threshold — CI runners are noisy, and false alarms get ignored.

How to Profile a Slow Python GIS Script and Find the Real Bottleneck #

Problem statement #

Quick answer #

Step-by-step solution #

1. Time the phases before profiling anything #

2. Profile the phase that dominates #

3. Narrow to the line #

4. Profile memory separately — it does not show in time profiles #

5. Profile a running process when you cannot restart it #

6. Confirm the fix with a measurement, not a feeling #

Code examples #

Example 1: a profiling harness for a pipeline #

Example 2: finding which geometries are expensive #

Example 3: automated regression detection #

Explanation #

Edge cases or notes #

Internal links #

FAQ #

Where should I start? #

What is the difference between tottime and cumtime? #

Why does my time profile not explain the slowness? #

How do I profile a job that is already running? #

Does profiling change the timings? #

How do I know my optimisation actually worked? #

Can I catch performance regressions automatically? #