Spatial Dashboard Streamlit · Panel · Production
Maps & Components

Filtering a Layer by a User-Drawn Polygon

Run the index pass, then the exact predicate on its candidates, then the attribute masks — and put the cache boundary between the second and third. A traced corridor returns ten times the candidates a rectangle does for the same answer, which is worth knowing before somebody reports it as a bug.

Why this matters

Filtering by a drawn shape is mechanically the same as filtering by a viewport, and it behaves differently enough in practice that reusing the viewport code without thinking produces a feature that is fine in testing and slow in use.

The difference is selectivity. A viewport is an axis-aligned rectangle, so the bounding box the spatial index works with is the query shape, and the cheap first pass is very nearly the whole answer. A shape somebody traced around a catchment, a corridor, or a district with a notch in it has a bounding box that can be several times its own area — and the index can only ever prune by bounding box. Everything inside that box and outside the shape comes back as a candidate and has to be rejected by the exact predicate, one feature at a time.

The two-pass filter a drawn shape needsA drawn polygon runs through the same two passes as a viewport, and for the same reason. The first pass asks the spatial index which features have bounding boxes overlapping the shape's bounding box, which is cheap and returns a superset. The second runs the exact predicate against only those candidates. The difference from a viewport is the selectivity: a viewport is a rectangle, so its bounding box is itself, and the first pass is nearly exact. A hand-drawn shape is usually long, curved or L-shaped, so its bounding box can be several times its area and the first pass returns many features the second must reject. That is not a fault in the index; it is the reason the second pass exists, and it is why treating the candidate list as the answer is wrong for exactly the shapes users draw.DRAWN SHAPE → CANDIDATES → MATCHESfiltersindex.query(shape)exact intersectsrender subsetA rectangle's bounding box is the rectangle. A drawn shape's can be several times its area — which is precisely whythe candidate list is never the answer.

Prerequisites

  • The validated, simplified shape from storing drawn polygons in session state. A raw freehand trace with four thousand vertices makes every exact predicate call expensive in proportion.
  • A layer with a built .sindex, held by whatever caches the layer, per speeding up GeoPandas spatial joins.
  • Both the shape and the layer in the same CRS. This is not a detail — a mismatch here returns an empty frame with no error at all.

Step-by-step solution

Step 1 — Align, then reduce, then test

python
import geopandas as gpd
from shapely.geometry.base import BaseGeometry

def filter_by_shape(gdf: gpd.GeoDataFrame, shape: BaseGeometry,
                    predicate: str = "intersects") -> gpd.GeoDataFrame:
    """Two passes: bounding-box candidates from the index, then the exact test."""
    assert gdf.crs is not None, "layer has no CRS — the predicate would return nothing"
    candidates = gdf.sindex.query(shape, predicate=None)      # bbox pass only
    subset = gdf.iloc[candidates]
    exact = getattr(subset.geometry, predicate)(shape)        # intersects / within
    return subset[exact].reset_index(drop=True)

Passing predicate=None to sindex.query asks for the bounding-box pass alone, which is what makes the two stages visible and measurable. Passing the predicate straight to query lets GeoPandas do both in one call and is what you would ship — the split here is so you can count the candidates, which is the number that explains the performance.

Step 2 — Know what the shape costs before the user draws it

Selectivity by shape, for the same feature countFour query shapes of identical area are run against the same layer of two hundred thousand points and the number of candidates returned by the bounding-box pass is compared with the number that survive the exact test. A rectangle returns 4,100 candidates of which 4,100 match, because its bounding box is itself. A compact drawn blob returns 5,200 of which 4,050 match, a modest overshoot. An L-shaped area returns 9,800 of which 4,100 match, because the bounding box covers the missing corner as well. A long river corridor returns 41,000 of which 4,000 match — a ten-to-one overshoot, since an axis-aligned box around a diagonal strip is mostly empty. The exact pass therefore does ten times more work for the corridor than for the rectangle, which is the cost worth knowing about before somebody traces a river and reports the dashboard as slow.SAME AREA, FOUR SHAPES — CANDIDATES vs MATCHESrectangle4,100 candidates → 4,100 matchcompact blob5,200 → 4,050L-shaped area9,800 → 4,100river corridor41,000 → 4,000a ten-to-one overshootThe exact pass does ten times the work for the corridor, for the same answer size. Worth knowing before a usertraces a river and files a performance ticket.

There is no fix for a corridor’s selectivity, because the index has only bounding boxes to work with. What there is instead is a mitigation: for shapes whose bounding box is much larger than their area, splitting the query into several smaller boxes along the shape and unioning the candidate sets prunes far more. It is worth doing only when the ratio is extreme, and the ratio is one line to compute:

python
overshoot = shape.envelope.area / shape.area          # 1.0 for a rectangle
if overshoot > 6:
    parts = split_along_shape(shape)                  # a handful of sub-boxes
    candidates = np.unique(np.concatenate([gdf.sindex.query(p) for p in parts]))

Step 3 — Order the predicates, and put the cache boundary in the right place

Combining a drawn shape with the other filtersA drawn area is one predicate among several, and the order they run in decides both speed and what can be cached. The spatial index pass runs first because it is the only one with an index behind it and the only one that reduces the candidate set logarithmically. The exact geometric predicate runs second, against candidates rather than against the layer. The attribute predicates — category, date, status — run last, as boolean masks over what is by then a small frame, costing microseconds. The cache boundary sits between the second and third: the spatial result depends only on the drawn shape, so it is worth caching under that shape's hash, while the attribute masks change every time a dropdown moves and are cheap enough to recompute. Reverse the order and both properties are lost at once.PREDICATE ORDER, AND WHERE THE CACHE BOUNDARY SITS1 · sindex.query(drawn shape)the only pass with an index behind it — logarithmic, and it does the real reduction2 · exact intersects on candidatesexpensive per feature, cheap in total because there are few candidates left— cache boundary —everything above depends only on the drawn shape, so it caches under that shape's hash3 · category, date, status masksboolean masks over a small frame; microseconds, and they change on every dropdownRun the attribute masks first and you scan the whole layer with no index, and produce a cache entry keyed to adropdown value that is evicted the moment anybody changes their mind.
python
import hashlib
import streamlit as st

@st.cache_data(ttl=900, max_entries=16, show_spinner=False)
def spatial_pass(shape_wkt: str) -> gpd.GeoDataFrame:
    """Cached on the shape alone — attribute filters stay outside."""
    from shapely import wkt
    return filter_by_shape(load_layer(), wkt.loads(shape_wkt))

subset = spatial_pass(st.session_state.aoi_wkt)
if category != "All":
    subset = subset[subset["category"] == category]     # microseconds, uncached

Keying the cache on the WKT string means an analyst who draws a shape, changes their mind about the category three times, and changes it back pays for the spatial work once. Keying it on the whole filter set would produce a separate entry per combination, each holding its own copy of the geometry, all of them discarded the moment the shape changes.

Step 4 — Say which predicate the map is using

intersects and within give visibly different answers for features that straddle the drawn boundary, and users have a firm intuition about which one they meant that does not survive contact with a shape they drew slightly too small. Put it in the interface:

python
mode = st.radio("Include features that…", ["touch the area", "are fully inside"],
                horizontal=True)
predicate = "intersects" if mode.startswith("touch") else "within"

Two words of interface removes an entire category of “the filter is wrong” reports, because the disagreement was never about the code.

Verification

python
from shapely.geometry import box

# 1. A rectangle's index pass is exact, so both passes agree.
rect = box(-0.2, 51.45, -0.05, 51.55)
cands = layer.sindex.query(rect, predicate=None)
assert len(layer.iloc[cands][layer.iloc[cands].intersects(rect)]) == len(cands)

# 2. The two-pass result equals the naive full scan.
fast = filter_by_shape(layer, drawn)
slow = layer[layer.intersects(drawn)]
assert len(fast) == len(slow), "the index pass dropped a match"

# 3. within is a strict subset of intersects.
assert len(filter_by_shape(layer, drawn, "within")) <= len(fast)

The second assertion is the one that matters, and it is the one that catches a stale index built before a reprojection — the failure described in speeding up spatial joins, where the query is fast and wrong and nothing raises.

Edge cases and gotchas

  • A shape drawn entirely outside the layer. Returns an empty frame, which is correct and looks identical to a CRS mismatch. Report “no features in the drawn area” explicitly so the two cases are distinguishable to the user.
  • A multi-part shape. If the draw control allows several polygons, decide whether they union or intersect. Union is almost always what users mean by drawing two areas, and the index pass takes the union’s parts individually rather than the union’s bounding box, which is much more selective.
  • Very large layers. Above a few million features the candidate pass itself becomes noticeable, and the answer is to push the predicate into PostGIS with a GiST index rather than filtering a resident frame at all.
  • Holes. A polygon with an interior ring behaves correctly with intersects, but its bounding box ignores the hole entirely, so the overshoot ratio is high and the exact pass rejects everything in the middle. Expected, and worth knowing when a doughnut-shaped catchment is slow.

FAQ

Should the exact predicate run in Python or in the database?

In the database when the layer is large enough that holding it resident is a problem, and in Python when it is already cached. The crossover is roughly where the layer stops fitting comfortably in the container’s memory budget — below it, a resident frame with an R-tree beats a round trip; above it, PostGIS with a GiST index wins on every axis.

Why does within return fewer features than users expect?

Because a feature has to be entirely inside the drawn shape, and a hand-drawn boundary almost always clips something the user thought they had enclosed. intersects matches the intuition of “I drew around these” far more often, which is why it is the better default — but say which one is active rather than relying on the default being intuitive.

Can I reuse the viewport filter code?

Yes — the two-pass structure is identical. What does not carry over is the assumption that the first pass is nearly exact, which holds for a rectangle and not for a traced shape. If the viewport code logs candidate counts, the drawn filter will show you the difference immediately.

Back to Dynamic Spatial Filtering.