How to Add Bounding Box and Attribute Filters to a Spatial API
On this page
Problem statement
An endpoint that returns everything is a download with extra latency. Measured on Natural Earth's 4,596 states and provinces, one client making requests back to back (Example 4):
GET /features?limit=200 3,380.2 kB (1,174.3 kB gzipped) 14.2 req/s
GET /features?limit=50 859.3 kB ( 300.9 kB gzipped) 47.6 req/s
GET /features?bbox=<London> 116.8 kB ( 36.0 kB gzipped) 118.7 req/s
The whole layer — what an endpoint without a cap would send — is 54.29 MB (Example 2). A London-sized bounding box is 116.8 kB of it, and served eight times as many requests per second as the largest page the cap allows.
Filters are therefore not a feature to add when the API gets slow — they are the reason it is an API. The work is in doing them safely: parsing untrusted input, choosing the right predicate, deciding when an index is worth having, and refusing requests that would return too much.
Quick answer
Parse, validate, filter, and cap:
%%writefile bbox_app.py
from fastapi import FastAPI, HTTPException, Query, Response
import geopandas as gpd
app = FastAPI()
LAYER = gpd.read_file("provinces.gpkg")
MAX_LIMIT = 200
def parse_bbox(raw: str) -> tuple[float, float, float, float]:
try:
x0, y0, x1, y1 = (float(v) for v in raw.split(","))
except ValueError:
raise HTTPException(400, "bbox must be minx,miny,maxx,maxy")
if x0 >= x1 or y0 >= y1:
raise HTTPException(400, "bbox min must be less than max")
if not (-180 <= x0 <= 180 and -180 <= x1 <= 180
and -90 <= y0 <= 90 and -90 <= y1 <= 90):
raise HTTPException(400, "bbox is outside the valid CRS84 range")
return x0, y0, x1, y1
@app.get("/features")
def features(bbox: str | None = None, admin: str | None = None,
limit: int = Query(50, ge=1, le=MAX_LIMIT), offset: int = 0):
subset = LAYER
if bbox:
x0, y0, x1, y1 = parse_bbox(bbox)
subset = subset.cx[x0:x1, y0:y1] # uses the spatial index
if admin:
subset = subset[subset["admin"] == admin]
page = subset.iloc[offset:offset + limit]
return Response(page.to_json(), media_type="application/geo+json",
headers={"X-Total-Count": str(len(subset))})bbox is minx,miny,maxx,maxy in longitude-first order, matching GeoJSON and OGC API Features. The service is a file so that it can run under uvicorn and be measured; Example 4 does that.
What you get — from 2. Use the spatial index — and know when it stops mattering
4,596 polygons: .cx 0.39 ms, brute-force intersects 0.38 ms
Validate the box, use the index, choose the predicate
1. Validate the bounding box before using it
Four failures worth rejecting explicitly, because each one produces a confusing result rather than an error:
- Wrong number of values — three or five numbers.
- Reversed order —
maxxbeforeminxselects nothing, silently. - Out of range — a latitude of 400, or coordinates in metres sent to a degrees API.
- Enormous — a bounding box covering the world, which is the unfiltered request wearing a filter.
Reject with a 400 and a message that says what was expected. A silently empty response is the least useful outcome available.
2. Use the spatial index — and know when it stops mattering
GeoPandas' .cx[] indexer uses the layer's spatial index. At small sizes that buys little:
import timeit
import geopandas as gpd
from shapely.geometry import box
LAYER = gpd.read_file("provinces.gpkg")
LONDON = (-0.51, 51.28, 0.33, 51.69)
x0, y0, x1, y1 = LONDON
indexed = timeit.timeit(lambda: LAYER.cx[x0:x1, y0:y1], number=200) / 200
brute = timeit.timeit(lambda: LAYER[LAYER.intersects(box(*LONDON))], number=200) / 200
print(f"{len(LAYER):,} polygons: .cx {indexed * 1000:.2f} ms, brute-force intersects {brute * 1000:.2f} ms")4,596 polygons: .cx 0.39 ms, brute-force intersects 0.38 ms
The index earns its place as the layer grows, and the crossover is well above a few thousand features. Build it once at start-up rather than per request:
LAYER.sindex # touching it builds the index; do this at import time3. Choose the predicate deliberately
.cx[] selects by bounding box intersection, which is a coarse filter: a feature whose envelope overlaps the box is returned even if the geometry does not.
For a viewport that is exactly right and much cheaper. For "features actually inside this polygon" it is wrong, and the exact test has to follow:
from shapely.geometry import box
candidates = LAYER.cx[x0:x1, y0:y1] # cheap envelope filter
exact = candidates[candidates.intersects(box(x0, y0, x1, y1))]
print(f"{len(candidates)} candidates, {len(exact)} actually intersect the box")39 candidates, 39 actually intersect the box
Here every candidate passed: provinces are compact, so their envelopes hug their shapes. The difference appears with long or diagonal features, as in the explanation below.
Two stages, cheap first. That is the same structure a spatial index uses internally.
4. Add attribute filters with an allowlist
Never interpolate a client-supplied column name into a query. Map the parameters you accept to the columns they filter:
FILTERABLE = {"admin": "admin", "type": "type", "region": "region"}
query_params = {"admin": "United Kingdom", "type": "London Borough", "column": "anything"}
subset = LAYER
for parameter, column in FILTERABLE.items():
value = query_params.get(parameter)
if value is not None:
subset = subset[subset[column] == value]
print(len(subset), "London boroughs; the unlisted 'column' parameter was ignored")28 London boroughs; the unlisted 'column' parameter was ignored
The dictionary is the API's documented surface. It also stops a filter on an unindexed column becoming a full scan by accident.
5. Cap the response, not just the page
A page limit bounds the row count and not the payload: provinces vary so much in detail that a page of 200 and a page of 50 differ by more than the factor of four in their counts would suggest — both are in the table at the top of the page.
Two additional caps are worth having:
- A maximum matched count without a bounding box — above it, require one.
- A coordinate precision limit. Rounding to six decimal places took the whole layer from 54.29 MB to 32.38 MB (Example 2).
6. Tell the client how much matched
headers={"X-Total-Count": str(len(subset))}or, for an OGC API Features service, numberMatched in the body. Either way the client can decide to narrow the filter rather than paging through four hundred pages.
A validated filter pipeline over world provinces
The data used on this page
The layer is Natural Earth's 1:10m admin-1 states and provinces — 4,596 polygons with their name, country (admin), type and region — saved as provinces.gpkg, and loaded into DuckDB for Example 3.
import io
import os
import zipfile
import duckdb
import geopandas as gpd
import requests
if not os.path.exists("provinces.gpkg"):
archive = requests.get("https://naciscdn.org/naturalearth/10m/cultural/"
"ne_10m_admin_1_states_provinces.zip", timeout=300)
with zipfile.ZipFile(io.BytesIO(archive.content)) as z:
for ext in (".shp", ".shx", ".dbf", ".prj", ".cpg"):
open(f"ne_provinces{ext}", "wb").write(z.read(f"ne_10m_admin_1_states_provinces{ext}"))
gpd.read_file("ne_provinces.shp")[["name", "admin", "type", "region", "geometry"]].to_file("provinces.gpkg")
if not os.path.exists("provinces.duckdb"):
build = duckdb.connect("provinces.duckdb")
build.execute("install spatial; load spatial;")
build.execute("create table provinces as select name, admin, geom from st_read('provinces.gpkg')")
build.close()
print(len(gpd.read_file("provinces.gpkg")), "provinces")4596 provinces
Example 1 — a filter pipeline with validation and limits
from dataclasses import dataclass
from fastapi import HTTPException
from shapely.geometry import box
import geopandas as gpd
@dataclass
class FilterLimits:
max_bbox_area_deg2: float = 100.0
max_unfiltered_features: int = 1000
max_limit: int = 200
coordinate_precision: int = 6
FILTERABLE = {"admin": "admin", "type": "type"}
def apply_filters(layer: gpd.GeoDataFrame, params: dict, limits: FilterLimits):
subset = layer
if (raw := params.get("bbox")):
x0, y0, x1, y1 = parse_bbox(raw)
area = (x1 - x0) * (y1 - y0)
if area > limits.max_bbox_area_deg2:
raise HTTPException(400, f"bbox covers {area:.1f} deg², maximum is "
f"{limits.max_bbox_area_deg2}")
subset = subset.cx[x0:x1, y0:y1]
if params.get("exact") == "true":
subset = subset[subset.intersects(box(x0, y0, x1, y1))]
for parameter, column in FILTERABLE.items():
if (value := params.get(parameter)) is not None:
if column not in subset.columns:
raise HTTPException(500, f"filterable column {column!r} is missing")
subset = subset[subset[column] == value]
if not params.get("bbox") and len(subset) > limits.max_unfiltered_features:
raise HTTPException(
400, f"{len(subset):,} features match; add a bbox or a filter "
f"(maximum without one is {limits.max_unfiltered_features:,})")
return subsetfrom bbox_app import parse_bbox
limits = FilterLimits()
print(len(apply_filters(LAYER, {"bbox": "-0.51,51.28,0.33,51.69", "exact": "true"}, limits)), "provinces in the London box")
for params in ({"admin": "Brazil"}, {"bbox": "-180,-60,180,80"}, {}):
try:
print(len(apply_filters(LAYER, params, limits)), "for", params)
except HTTPException as exc:
print(exc.status_code, exc.detail)39 provinces in the London box
27 for {'admin': 'Brazil'}
400 bbox covers 50400.0 deg², maximum is 100.0
400 4,596 features match; add a bbox or a filter (maximum without one is 1,000)
Example 2 — reducing the payload before serialising
import shapely
def prepare_response(gdf, precision=6, keep_columns=None, simplify_deg=None):
"""Three cheap reductions, applied in the order that compounds best."""
out = gdf
if keep_columns:
out = out[[*keep_columns, out.geometry.name]]
if simplify_deg:
out = out.copy()
out["geometry"] = out.geometry.simplify(simplify_deg, preserve_topology=True)
if precision is not None:
out = out.copy()
out["geometry"] = shapely.set_precision(out.geometry.values,
10 ** -precision)
return outMeasured on the full 4,596-feature layer, the reductions individually:
import gzip
for label, frame in (("as-is", LAYER),
("6 decimal places", prepare_response(LAYER, precision=6)),
("simplified to 0.01°", prepare_response(LAYER, precision=None, simplify_deg=0.01))):
text = frame.to_json().encode()
print(f"{label:22} {len(text) / 1e6:6.2f} MB ({len(gzip.compress(text, 6)) / 1e6:5.2f} MB gzipped)")as-is 54.29 MB (18.10 MB gzipped)
6 decimal places 32.38 MB (10.15 MB gzipped)
simplified to 0.01° 16.62 MB ( 5.34 MB gzipped)
Precision alone removed 40% and changes nothing a consumer can detect at map scales. Simplification removes more and does change the geometry, so it belongs on a drawing endpoint rather than on a query endpoint.
Example 3 — pushing the filter into the database instead
import duckdb
def query_features(con, bbox=None, admin=None, limit=50, offset=0):
"""For datasets too large to hold in memory, filter in the engine."""
clauses, params = [], []
if bbox:
x0, y0, x1, y1 = bbox
clauses.append("st_intersects(geom, st_makeenvelope(?, ?, ?, ?))")
params += [x0, y0, x1, y1]
if admin:
clauses.append("admin = ?")
params.append(admin)
where = ("where " + " and ".join(clauses)) if clauses else ""
matched = con.execute(
f"select count(*) from provinces {where}", params).fetchone()[0]
rows = con.execute(f"""
select name, admin, st_asgeojson(geom) as geometry
from provinces {where}
limit ? offset ?
""", params + [limit, offset]).fetchall()
return matched, rowscon = duckdb.connect("provinces.duckdb", read_only=True)
con.execute("load spatial")
matched, rows = query_features(con, bbox=LONDON, admin="United Kingdom", limit=3)
print(matched, "matched;", [row[0] for row in rows])39 matched; ['Essex', 'Thurrock', 'Kent']
Parameter placeholders rather than string formatting, for the same reason as in any other SQL: a bounding box is client input.
Example 4 — measuring the service
import socket
import subprocess
import sys
import time
import httpx
with socket.socket() as probe:
probe.bind(("127.0.0.1", 0))
port = probe.getsockname()[1]
server = subprocess.Popen([sys.executable, "-m", "uvicorn", "bbox_app:app", "--port", str(port),
"--log-level", "warning"], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL)
base = f"http://127.0.0.1:{port}"
for _ in range(200):
try:
httpx.get(base + "/docs", timeout=1); break
except httpx.TransportError:
time.sleep(0.1)
def throughput(path, seconds=5.0):
"""Sequential requests from one client for a few seconds."""
n, start = 0, time.perf_counter()
while time.perf_counter() - start < seconds:
body = httpx.get(base + path, timeout=None).content
n += 1
return body, n / (time.perf_counter() - start)
for path in ("/features?limit=200", "/features?limit=50", "/features?bbox=-0.51,51.28,0.33,51.69"):
body, rate = throughput(path)
print(f"GET {path:42} {len(body) / 1e3:9.1f} kB ({len(gzip.compress(body, 6)) / 1e3:7.1f} kB gzipped) "
f"{rate:6.1f} req/s")
server.terminate()
server.wait()GET /features?limit=200 3380.2 kB ( 1174.3 kB gzipped) 14.2 req/s
GET /features?limit=50 859.3 kB ( 300.9 kB gzipped) 47.6 req/s
GET /features?bbox=-0.51,51.28,0.33,51.69 116.8 kB ( 36.0 kB gzipped) 118.7 req/s
Why the bounding box beats every other optimisation
Why the bounding box is worth more than any other optimisation
Every other technique on this page reduces the response by a factor of two or three. The bounding box reduces it by whatever fraction of the world the client is looking at — measured, from 54.29 MB for the whole layer to 116.8 kB for a London-sized window.
That is not an optimisation, it is a different question being asked. The API returns what the client needs instead of what the dataset contains, and everything downstream — bandwidth, serialisation, browser memory — scales with it.
Why the spatial index does not always help
An index turns a linear scan into a tree descent, and at 4,596 features the two measured the same: 0.39 ms and 0.38 ms.
The crossover depends on the geometry and the machine, and it is well above a few thousand features. The practical rule: build the index once at start-up, because it costs nothing to have and matters a great deal at a hundred thousand features.
Why the envelope filter is not the exact answer
.cx[] compares bounding boxes. A long diagonal river's envelope covers a large rectangle, most of which the river does not touch, so an envelope query returns it for boxes it does not intersect at all.
For a map viewport that is harmless and desirable — better to draw a feature slightly outside the view than to miss one. For "which features are in this administrative area?", the exact predicate is required, and running it only on the envelope-filtered candidates keeps it cheap.
Why a cap belongs in the first version
Without one, the API's worst case is the whole dataset, and the worst case is what a misbehaving client will find. Measured, the largest page the service allows — 3.4 MB — served 14.2 requests per second to one client; a service without the cap would serve the whole layer to anyone who asked.
A cap converts that into a 400 with an explanation, which is a much better failure. It also documents the API's intended use: an endpoint that requires a bounding box above a thousand features is telling clients what it is for.
Axis order, the antimeridian and attribute allowlists
bboxis longitude-first —minx,miny,maxx,maxy— in CRS84.- A bounding box crossing the antimeridian has
minx > maxx; either support it explicitly or reject it. - A bbox in the wrong CRS silently selects nothing; validate the range.
.cx[]is envelope-based. Follow with an exact predicate where correctness requires it.- Attribute filters need an allowlist, never a client-supplied column name.
- Build
sindexat start-up, not per request. X-Total-CountornumberMatchedlets clients avoid pointless paging.- Round coordinates. Six decimal places is 11 cm and cut the whole layer by 40%.
Internal links
- Serving spatial data explained: files, features and tiles — where filters sit in the design
- OGC API Features explained — the standard parameter names
- How to paginate a large feature API in Python — the other half of bounding a response
- How to build a GeoJSON API with FastAPI — the base implementation
- Fixing a spatial API that is slow under load — when the filters are not enough
- How to use a spatial index in GeoPandas — what
.cx[]uses - How to select features by location — the predicates
- Fixing an API that returns empty or truncated GeoJSON — when a filter silently matches nothing
FAQ
How much does a bounding box filter save?
Measured on a 4,596-feature dataset: the whole layer is 54.29 MB; a London-sized bounding box returned 116.8 kB at 118.7 requests per second, against 14.2 for a 200-feature page.
What order are bounding box coordinates in?
minx,miny,maxx,maxy — longitude first — following GeoJSON and OGC API Features, which use CRS84 rather than EPSG:4326's official latitude-first order.
Do I need a spatial index?
Build it at start-up because it is free to have, but do not expect much at small sizes: on 4,596 polygons, .cx took 0.39 ms and a brute-force intersects 0.38 ms.
Is .cx[] an exact spatial filter?
No, it compares bounding boxes. For a viewport that is what you want; for exact containment, follow it with an intersects or within test on the candidates.
How do I filter on attributes safely?
With an allowlist mapping accepted query parameters to columns. Never interpolate a client-supplied column name.
What should happen when a request would return too much?
A 400 with a message asking for a bounding box or a narrower filter. The alternative is the whole layer, which any client can trigger.