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

Problem statement

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

The pipeline on this page is an ordinary one: read HM Land Registry's 293,729 title polygons for Leeds, clean them, find the census area (LSOA) each one falls in, count parcels per area, write the result. It has a Python loop in the join, and anyone reading the code would bet the join is the slow part. On this machine it takes about three minutes, and the join is under a quarter of that.

Intuition is unreliable here for a specific reason. GIS code mixes several cost centres — Python-level iteration, GEOS predicate work, file I/O, the network, memory pressure — and they produce similar-looking symptoms while needing opposite fixes. Chunking a CPU-bound loop makes it slower. Vectorising a script that is waiting on something 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).strip_dirs().sort_stats("cumulative").print_stats(20)
print(s.getvalue())

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.

read            139.9 s   76.4%
join             42.6 s   23.3%
clean             0.6 s    0.3%
write             0.0 s    0.0%
aggregate         0.0 s    0.0%

Time the phases, then profile, narrow and measure memory

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, columns=["INSPIREID"])
with phase("clean"):
    parcels = clean(parcels)
with phase("join"):
    parcels = assign_zones(parcels, zones)
with phase("aggregate"):
    summary = summarise(parcels, zones)
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}%")

read            139.9 s   76.4%
join             42.6 s   23.3%
clean             0.6 s    0.3%
write             0.0 s    0.0%
aggregate         0.0 s    0.0%

Eighteen lines, and the bet was wrong: reading a 46 MB download is most of the run, and the loop everyone would have rewritten first is a fraction of it.

Amdahl's law is the reason this matters. Making the join infinitely fast would save only its share of the run. Fixing whatever the read is doing is worth more than every other phase put together.

2. Profile the phase that dominates

import cProfile
import pstats


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

Profile the read — and, because a Python profiler cannot see inside GDAL, ask GDAL to log what it does at the same time:

import pyogrio

pyogrio.set_gdal_config_options({"CPL_DEBUG": "ON", "CPL_LOG": os.path.abspath("gdal.log")})
parcels = profile_call(gpd.read_file, PARCELS, columns=["INSPIREID"], sort="tottime", top=4)
pyogrio.set_gdal_config_options({"CPL_DEBUG": None, "CPL_LOG": None})

print("\n".join(line.strip()[:150] for line in open("gdal.log") if "HTTP" in line))

         4252 function calls (3985 primitive calls) in 140.816 seconds

   Ordered by: internal time
   List reduced from 388 to 4 due to restriction <4>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        1  140.573  140.573  140.574  140.574 raw.py:39(read)
        1    0.213    0.213    0.213    0.213 io.py:308(from_wkb)
        1    0.021    0.021  140.816  140.816 file.py:482(_read_file_pyogrio)
        9    0.003    0.000    0.003    0.000 {method 'copy' of 'numpy.ndarray' objects}


HTTP: Fetch(http://192.168.4.19:5080/geoserver/LR/wfs?service=WFS&version=2.0.0&request=DescribeFeatureType&typeName=LR%3APREDEFINED)

Two columns matter and they answer different questions:

• tottime — time in the function itself, excluding calls it makes. Almost all of it is in pyogrio's read, which is a single call into GDAL's C code.
• cumtime — time including everything it called: it finds the caller responsible.

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

Here cProfile runs out of road: the time is inside one C call it cannot see into. GDAL's debug log names the culprit. The GML file declares its schema at http://192.168.4.19:5080/… — a private-network address, presumably a server inside HM Land Registry — and GDAL tries to download it before reading a single feature. Where the network drops the request silently, GDAL waits until it times out; where it is refused at once, there is no stall — a re-run on 8 October 2026, after this machine's router began refusing it, read the file in 6.7 s, not 140.8 s. No guess would have blamed the network.

The fix is one setting, telling GDAL not to fetch schemas:

pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"})

t0 = time.perf_counter()
parcels = gpd.read_file(PARCELS, columns=["INSPIREID"])
print(f"read {time.perf_counter() - t0:.1f} s")

read 4.9 s

Now the join really is the biggest phase, so profile it — on a sample, because cProfile's per-call overhead is large on a loop of small calls:

parcels = clean(parcels)
sample = parcels.sample(20_000, random_state=0)
profile_call(assign_zones, sample, zones, sort="tottime", top=6, dump="join.prof")

         14859362 function calls (14439436 primitive calls) in 5.842 seconds

   Ordered by: internal time
   List reduced from 369 to 6 due to restriction <6>

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
    20000    0.599    0.000    0.599    0.000 predicates.py:551(contains)
    20000    0.312    0.000    0.318    0.000 compute.py:467(take)
3319956/3279954    0.309    0.000    0.481    0.000 {built-in method builtins.isinstance}
    59974    0.162    0.000    0.216    0.000 array.py:778(__getitem__)
919796/579874    0.132    0.000    0.199    0.000 {built-in method builtins.len}
    60002    0.121    0.000    0.169    0.000 generic.py:6119(__finalize__)

Twenty thousand calls to contains — one per parcel — each testing the point against all 488 zones, plus pandas indexing on every iteration. That is the diagnosis; the fix is a spatial join, measured in step 6.

Visualise a large profile rather than reading columns:

snakeviz join.prof            # opens an interactive view in a browser

3. Narrow to the line

cProfile reports per function, which is not enough when one function is many lines. line_profiler times each line; it can wrap a function directly, with no decorator and no kernprof:

from line_profiler import LineProfiler

raw = gpd.read_file(PARCELS, columns=["INSPIREID"])
lp = LineProfiler(clean)
lp.runcall(clean, raw)
lp.print_stats(output_unit=1e-3, stripzeros=True)

Timer unit: 0.001 s

Total time: 0.557815 s
File: cell_297.py
Function: clean at line 37

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    37                                           def clean(gdf):
    38         1         21.0     21.0      3.8      gdf = gdf[gdf.geometry.notna()]
    39         1         16.4     16.4      2.9      gdf = gdf[~gdf.geometry.is_empty]
    40         1        393.8    393.8     70.6      gdf["valid"] = gdf.geometry.is_valid
    41         1          5.2      5.2      0.9      gdf.loc[~gdf["valid"], "geometry"] = gdf.loc[~gdf["valid"], "geometry"].make_valid()
    42         1         27.0     27.0      4.8      gdf["area"] = gdf.geometry.area
    43         1         12.2     12.2      2.2      gdf = gdf[gdf["area"] > 1.0]
    44         1         82.2     82.2     14.7      gdf["centroid"] = gdf.geometry.representative_point()
    45         1          0.0      0.0      0.0      return gdf

The line to look at is the one with the largest % Time, not the one that looks expensive. On these parcels make_valid is free, because the Land Registry polygons are all valid, and the cost is in testing validity at all and in the representative points. On a dirtier layer the same function would profile completely differently — which is the argument for profiling on your real data rather than a tidy sample.

line_profiler adds a lot of overhead per line, so use it on one function with a modest input, never on a whole run.

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

memray records every allocation, including those made in C:

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

The standard library's tracemalloc needs no install, but it only sees allocations made through Python's own allocator. Compare it with the process's resident memory on the same read:

import gc
import tracemalloc

import psutil

gc.collect()
rss_before = psutil.Process().memory_info().rss
tracemalloc.start()
parcels = gpd.read_file(PARCELS, columns=["INSPIREID"])
current, peak = tracemalloc.get_traced_memory()
tracemalloc.stop()
rss_after = psutil.Process().memory_info().rss
print(f"tracemalloc: peak {peak / 1e6:,.0f} MB, still held {current / 1e6:,.0f} MB")
print(f"process RSS grew by {(rss_after - rss_before) / 1e6:,.0f} MB")

tracemalloc: peak 147 MB, still held 20 MB
process RSS grew by 298 MB

The gap is GEOS: the 293,729 Shapely geometries live in C memory that tracemalloc never sees. For GIS work that is most of the memory, so tracemalloc alone gives a misleading number and memray (or RSS, crudely) a useful 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

py-spy reads the process's memory from outside, so there is no code change, no restart and almost no overhead, and the stack it prints names the line. On many Linux systems attaching to a process you did not start needs sudo — Ubuntu's default kernel.yama.ptrace_scope = 1 only lets a process trace its own children — so either run it with sudo, or start the job under it: py-spy record -o profile.svg -- python pipeline.py.

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. It would also have found the stall in step 2 without any debug logging: a stack stuck in a network call is hard to misread.

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

import statistics


def assign_zones_sjoin(parcels, zones):
    """The same result as assign_zones, with an indexed spatial join."""
    points = parcels.set_geometry("centroid")[["centroid"]]
    joined = gpd.sjoin(points, zones[["LSOA21CD", "geometry"]], how="left", predicate="within")
    joined = joined[~joined.index.duplicated(keep="first")]
    return parcels.assign(zone=joined["LSOA21CD"])


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__:<20} min {min(times):.3f}s  "
          f"median {statistics.median(times):.3f}s  max {max(times):.3f}s")
    return result


before = benchmark(assign_zones, sample, zones, runs=3)
after = benchmark(assign_zones_sjoin, sample, zones, runs=3)
same = (before["zone"].fillna("") == after["zone"].fillna("")).all()
print(f"same zones: {same}")
assert same, "the optimised version returns different results"

assign_zones         min 2.922s  median 2.927s  max 2.995s
assign_zones_sjoin   min 0.021s  median 0.023s  max 0.028s
same zones: True

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.

A profiling harness and the expensive geometries

The data used on this page

HM Land Registry's INSPIRE Index Polygons for Leeds — the extent of every registered freehold title, 293,729 polygons — and the 488 Lower Layer Super Output Areas of Leeds from the ONS, with their Census 2021 population from Nomis. The block also defines the pipeline's functions, so every step above runs against them. The parcels contain HM Land Registry data © Crown copyright and database right 2026, licensed under the Open Government Licence v3.0.

import io
import os
import time

import geopandas as gpd
import numpy as np
import pandas as pd
import requests

# HM Land Registry INSPIRE Index Polygons for Leeds. The download sets a cookie
# and redirects to itself, so it needs a session to keep the cookie.
reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/Leeds.zip",
                               timeout=600)
reply.raise_for_status()
with open("leeds_inspire.zip", "wb") as f:
    f.write(reply.content)
PARCELS = "zip://leeds_inspire.zip!Land_Registry_Cadastral_Parcels.gml"

# zones: the 488 Lower Layer Super Output Areas of Leeds, with Census 2021 population
ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "Lower_layer_Super_Output_Areas_December_2021_Boundaries_EW_BGC_V5/FeatureServer/0/query")
features, offset = [], 0
while True:                                   # the service pages at 1,000 features
    page = requests.get(ONS, params={"where": "LSOA21NM LIKE 'Leeds %'", "outFields": "LSOA21CD",
                                     "outSR": 27700, "f": "geojson", "resultOffset": offset}, timeout=120).json()
    features += page["features"]
    if not page.get("properties", {}).get("exceededTransferLimit"):
        break
    offset += len(page["features"])
zones = gpd.GeoDataFrame.from_features(features, crs=27700).sort_values("LSOA21CD").reset_index(drop=True)
population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(zones['LSOA21CD'])}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
zones["population"] = zones["LSOA21CD"].map(population)


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


def assign_zones(parcels, zones):
    """The zone of each parcel — one Python iteration per parcel, no index."""
    zone = []
    for point in parcels["centroid"]:
        hit = zones.loc[zones.contains(point), "LSOA21CD"]
        zone.append(hit.iloc[0] if len(hit) else None)
    return parcels.assign(zone=zone)


def summarise(parcels, zones):
    counts = parcels.groupby("zone").size()
    out = zones.set_index("LSOA21CD")
    out["parcels"] = counts.reindex(out.index, fill_value=0)
    out["parcels_per_1000"] = 1000 * out["parcels"] / out["population"]
    return out.reset_index()


print(f"{os.path.getsize('leeds_inspire.zip') / 1e6:.0f} MB of parcels downloaded; "
      f"{len(zones)} zones, {zones['population'].sum():,} residents")

46 MB of parcels downloaded; 488 zones, 811,957 residents

Example 1: a profiling harness for a pipeline

import gc
import os
import time
from contextlib import contextmanager
from dataclasses import dataclass, field


@dataclass
class Profiler:
    phases: dict = field(default_factory=dict)
    rows: dict = field(default_factory=dict)
    growth: 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
            self.phases[name] = self.phases.get(name, 0.0) + elapsed
            self.growth[name] = max(self.growth.get(name, 0.0), self._rss() - before)
            if rows is not None:
                self.rows[name] = rows

    def report(self):
        total = sum(self.phases.values()) or 1.0
        print(f"{'phase':<12}{'seconds':>10}{'share':>8}{'ΔRSS':>9}{'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:<12}{secs:>9.1f}s{100 * secs / total:>7.1f}%"
                  f"{self.growth[name]:>8.0f}M{per}")
        print(f"{'TOTAL':<12}{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"):
    parcels = gpd.read_file(PARCELS, columns=["INSPIREID"])      # the schema fetch is off now
with prof.phase("clean", rows=len(parcels)):
    parcels = clean(parcels)
with prof.phase("reproject", rows=len(parcels)):
    parcels_wgs84 = parcels.to_crs(4326)
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             42.7s   80.3%      17M      145.5 ms
read              6.1s   11.5%       0M              
reproject         2.5s    4.6%       9M        8.4 ms
write             1.0s    1.9%     185M        3.4 ms
clean             0.9s    1.6%       1M        2.9 ms
TOTAL            53.2s

  → 'join' is 80% of the run. Halving it saves 40% overall; halving everything else saves 10%.

Two columns beyond the obvious earn their place. ΔRSS shows memory a time profile would never surface — which phase needs the big 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 how much optimising anything else could possibly buy.

Example 2: finding which geometries are expensive

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

import numpy as np
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[:3]:
        print(f"    {subset.index[i]:>7}  {verts[i]:>8,} vertices  {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


times, verts = profile_by_geometry(parcels, lambda g: g.buffer(50))

sampled 2,000 of 293,485 features
  total sampled time 0.07s
  slowest 1% account for 18%
  correlation(cost, vertices) = 0.85

  slowest features:
     148658       337 vertices       1.3 ms
      34958       455 vertices       1.3 ms
     226787       552 vertices       1.2 ms

  projected full-layer cost: 11s

The correlation says how far cost follows vertex count, and the share taken by the slowest 1% says whether a handful of features dominate. When they do, simplifying or subdividing those few is the fix — and the per-feature timing is the only way to see it, because cProfile reports one aggregate number for every call to buffer.

Example 3: automated regression detection

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

import json
import platform
import 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 = []
    for name, secs in current.items():
        was = baseline.get(name)
        if was is None:
            print(f"  new    {name:<14} {secs:>8.2f}s")
            continue
        change = (secs - was) / was
        flag = "SLOWER" if change > tolerance else "faster" if change < -tolerance else "ok"
        print(f"  {flag:<7}{name:<14} {was:>8.2f}s → {secs:>8.2f}s ({change:+.0%})")
        if change > tolerance:
            regressions.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"  baseline written to {BASELINE}")
    elif regressions:
        raise SystemExit(
            f"{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


check_performance([timed_run("join_zones", assign_zones_sjoin, sample, zones)])     # the baseline

# later, a refactor quietly brings the loop back
try:
    check_performance([timed_run("join_zones", assign_zones, sample, zones)])
except SystemExit as failure:
    print(failure)

  new    join_zones         0.02s
  baseline written to perf_baseline.json
  SLOWER join_zones         0.02s →     2.88s (+13000%)
1 performance regression(s) beyond 25%:
  join_zones: 0.02s → 2.88s

A regression like this 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.

Intuition about performance is systematically wrong

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 76% of the run, as the read was here, no amount of work on the other 24% can improve the total by more than 24%. Phase timing takes eighteen lines and tells you where optimisation can pay, before any tooling is installed. Skipping it is how the join on this page would have been rewritten first — a real improvement, on less than a quarter of the problem.

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 often a GEOS method, and the useful question is which caller invokes it once per row. And sometimes the leader is a single call into C — pyogrio's read, here — where the profile's honest answer is "not in Python", and the next tool has to come from the library itself.

Sampling profilers and tracing profilers make different trade-offs. cProfile traces every call, which is accurate on call counts and adds overhead per call — on this page's join, 5.8 s profiled against 2.9 s unprofiled on the same 20,000 parcels, 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. On the Leeds read, tracemalloc saw a small fraction of what the process's resident memory grew by.

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.

Profiler overhead, C calls and tracemalloc's blind spot

• cProfile doubled the cost of the loop it profiled here and inflates the apparent cost of many small calls relative to a few large ones.
• line_profiler's overhead is per line executed — negligible on vectorised code like clean, heavy inside a Python loop. Use it on one function with a modest input.
• When the time is inside one C call, ask the library. GDAL's CPL_DEBUG log named the schema URL that cProfile could only call "read".
• 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 doubled the time of the loop profiled on this page and biases against many small calls; line_profiler costs more again inside loops. 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.