Profiling GEOS Predicate Hot Paths

9 min read

Profiling a spatial suite is different from profiling ordinary Python in one respect that changes everything: most of the time is spent inside a C extension, where a Python profiler sees a single opaque call. This guide sits beneath performance benchmarking spatial suites and covers reading a profile whose hot path is a GEOS call: which profiler to use, what the numbers mean when the work is invisible, and the two changes — prepared geometry and vertex reduction — that account for most of the available improvement.

The finding that recurs is worth stating first. Predicate cost is dominated by vertex count, not by feature count, and the two are uncorrelated in real data. A layer of a thousand highly-detailed polygons is far more expensive than one of ten thousand simple ones, and no amount of index tuning changes that.

Root cause: the profiler cannot see inside the predicate

A deterministic Python profiler records every call and return, which gives exact call counts and attributes time precisely — to the boundary of the C extension. Everything GEOS does inside intersects appears as time spent in one function, so the profile says which predicate is expensive and nothing about why.

That is less limiting than it sounds, because the actionable question is usually which call site dominates rather than which GEOS routine does. But it changes the reading: a profile showing 80 per cent of time in a single predicate call is not a dead end, it is an instruction to look at the inputs to that call.

The profiler's horizon in a spatial call stack A call stack drawn vertically. The upper three frames — the test function, a helper looping over features, and the predicate call itself — are shown as distinct with exact call counts and attributed durations, and are marked as visible to the profiler. A horizontal boundary marks the profiler's horizon at the C extension edge. Beneath it, the GEOS entry point, the noding step and the relate computation are drawn as a single undifferentiated block labelled opaque, because the profiler records only the total time spent in the predicate call. A closing note records that the actionable information is therefore the call count and the properties of the inputs, rather than any internal breakdown. test_no_overlaps · 1 call · 4.0 s check_pairs · 1 call · 3.9 s intersects · 41 203 calls · 3.7 s the profiler’s horizon — the C extension boundary GEOS: noding, relate, prepared cache one opaque block — no internal breakdown the actionable number The call count and the inputs are what the profile gives you. That is enough, because both available remedies act on exactly those.

Profiler reference

Tool Sees Overhead Right for
Deterministic call profiler Exact call counts, Python frames High — distorts short calls Finding the dominant call site
Sampling profiler Approximate distribution, native frames Low Confirming a hot path under realistic load
Native profiler Inside the C extension Moderate, needs symbols Only when the engine itself is suspect
Manual counters Exactly what you instrument Negligible A number you want in every run

For spatial work the first two cover almost everything. The deterministic profiler answers “how many times is this predicate called”, which is usually the actionable question; the sampling profiler answers “where is the time actually going” without distorting the measurement, which matters when the calls are short and numerous.

Step-by-step implementation

Step 1 — Count the calls before timing them

The call count is the number that leads somewhere, and it is cheap enough to record on every run rather than only when profiling.

import cProfile, pstats, io

def profile_call_counts(func, *args, top=10):
    pr = cProfile.Profile()
    pr.enable(); func(*args); pr.disable()
    s = io.StringIO()
    pstats.Stats(pr, stream=s).sort_stats("tottime").print_stats(top)
    return s.getvalue()

A predicate called forty thousand times when the dataset has a thousand features means the filter stage is not filtering — the same finding the candidate-count benchmark produces, arrived at from the other direction.

Step 2 — Use prepared geometry for repeated comparisons

When one geometry is compared against many, GEOS can build an internal index of it once and reuse it. This is the single largest available improvement in most spatial suites, and it is one line.

import shapely

def count_intersections(one, many):
    shapely.prepare(one)                  # builds the internal index once
    return sum(1 for g in many if shapely.intersects(one, g))

Without preparation the geometry is re-analysed on every comparison. With it, the analysis happens once and each comparison is a lookup — a difference that grows with the vertex count of the prepared geometry, which is exactly where the cost lives.

Step 3 — Measure the vertex count, not the feature count

Feature count is the number everybody records and vertex count is the one that predicts cost.

import shapely
import numpy as np

def vertex_profile(gdf) -> dict:
    counts = shapely.get_num_coordinates(gdf.geometry.values)
    return {
        "features": int(len(counts)),
        "vertices_total": int(counts.sum()),
        "vertices_p95": int(np.percentile(counts, 95)),
        "vertices_max": int(counts.max()),
    }

The 95th percentile is the useful one. A layer whose median feature has 40 vertices and whose 95th percentile has 12,000 will be dominated by a handful of features, and simplifying only those changes the profile more than any other single action.

Step 4 — Confirm with a sampling profiler under load

A deterministic profiler distorts short calls by adding per-call overhead, which for a predicate called tens of thousands of times can invert the ranking. Confirm the finding without that distortion.

# Attach a sampling profiler to the real run rather than a profiled one.
py-spy record -o profile.svg -- pytest -q tests/bench/test_join.py
Ten thousand simple features against one thousand detailed ones Two layers compared on three quantities. The first layer has ten thousand features, a ninety-fifth percentile vertex count of forty, a modest total vertex count, and a predicate pass that completes quickly. The second layer has only one thousand features but a ninety-fifth percentile vertex count of twelve thousand, a far larger total vertex count, and a predicate pass an order of magnitude slower despite having a tenth as many features. A closing note identifies the ninety-fifth percentile vertex count as the number worth recording, because the cost is dominated by a small number of highly detailed features rather than by the feature count. LAYER FEATURES P95 VERTICES PREDICATE PASS simple polygons administrative grid cells 10 000 40 0.4 s detailed polygons coastline, river network 1 000 12 000 6.1 s A tenth as many features, fifteen times the cost. Feature count predicts almost nothing about predicate expense. Record the 95th percentile vertex count — the cost is dominated by a small number of detailed features. Simplifying only those changes the profile more than any other single action available.

Verify the fix

Compare the profile before and after preparing the geometry:

pytest -q tests/bench/test_join.py --benchmark-only --benchmark-columns=median

Preparation should reduce the median substantially where one geometry is compared against many, and change nothing where each comparison involves a different pair — which is itself a useful confirmation that the call pattern is what you assumed.

The two changes that account for most of the improvement

Profiling a spatial suite tends to produce the same two findings, and knowing them in advance shortens the exercise considerably.

Prepared geometry, wherever one shape meets many. Any containment test against a boundary, any clip against a study area, any check of many features against one reference geometry. The improvement scales with the vertex count of the prepared shape, so it is largest exactly where the cost is worst.

Vertex reduction on the heavy tail. Not simplification of the whole layer — that changes the data — but a decision about whether the pipeline needs full detail at this stage. A validity check does; a bounding-box prefilter does not. Running cheap stages against a simplified copy and expensive ones against the original is often a large win with no loss of correctness, provided the simplified copy is used only where its use is sound.

A third finding appears less often and is worth naming because it is invisible otherwise: repeated construction of a transformer, an index, or a prepared geometry inside a loop. The profile shows many short calls whose total dominates, and the fix is to hoist the construction. This is the same class as the repeated-fixture finding, one level down.

Preparation moves the analysis out of the loop Two arrangements of the same comparison loop. Without preparation, every comparison re-analyses the reference geometry, so each iteration costs an amount proportional to that geometry's vertex count and the total is the vertex count multiplied by the number of comparisons. With preparation, the analysis is performed once before the loop and each comparison becomes an index lookup at small constant cost, so the total is one analysis plus a constant per comparison. A closing note observes that the saving grows with both the vertex count of the reference geometry and the number of comparisons, which means it is largest exactly where the cost was worst. Without prepare() each block re-analyses the reference geometry total ≈ vertices × comparisons With prepare() analyse once each comparison is an index lookup total ≈ one analysis + k × comparisons The saving grows with both the reference geometry’s vertex count and the number of comparisons — largest exactly where the cost was worst. It is one line, and it is the single largest improvement available in most spatial suites.

When to stop profiling

Profiling has diminishing returns and a recognisable stopping point. Once the two standard remedies have been applied — preparation where one shape meets many, and construction hoisted out of loops — the remaining profile is usually flat: no single call site dominates, and further work means algorithmic change rather than tuning.

That flatness is the signal to stop. A flat profile means the cost is distributed across genuinely necessary work, and the next available improvement is a different approach rather than a faster version of this one: pushing the operation into the database, reducing the data the stage sees, or accepting the cost and adjusting the budget.

Recognising it early matters because profiling is absorbing. It is always possible to find another five per cent, and the effort is rarely worth it once the profile has flattened — particularly against the alternative of moving an expensive check to the scheduled tier, which is a five-minute change with a larger effect than any micro-optimisation available.

Failure modes and edge cases

  1. Profiling with a deterministic profiler and trusting the ranking. Per-call overhead can invert the order when calls are short and numerous. Confirm with a sampling profiler before acting.
  2. Preparing a geometry used once. Preparation costs an index build; if the shape is compared once, that is pure overhead. Prepare only for repeated use.
  3. Simplifying the data rather than a working copy. A simplified geometry that reaches the output is a correctness change wearing a performance justification. Simplify a copy, and use it only where the loss is sound.
  4. Profiling a fixture-dominated run. If the fixture build is most of the time, the predicate profile is a small slice of a small slice. Profile the operation, not the test.
  5. Optimising from a profile taken at fixture scale. The ranking at a thousand features can differ from the ranking at a million, particularly where an index is involved. Profile at a size that resembles production.
  6. Ignoring memory. A predicate pass that fits in cache and one that does not differ by more than any algorithmic change available, and the profile shows only the symptom.

Conclusion

A spatial profile has a horizon at the C extension boundary, and the information above that line is enough: call counts identify the dominant site, and the inputs to that site identify the remedy. Recording vertex counts rather than feature counts, preparing geometry wherever one shape meets many, hoisting construction out of loops, and confirming with a sampling profiler covers most of the improvement available — which is what makes the measurements in performance benchmarking spatial suites actionable rather than merely descriptive.