⚙️ How-To APIs & tile services Performance & scale Serving Spatial Data · step 3 of 20 Updated

How to Add Bounding Box and Attribute Filters to a Spatial API

Problem statement

An endpoint that returns everything is a download with extra latency. The measured difference on a 4,596-feature dataset:

GET /features                      54.06 MB   (18.05 MB gzipped)   2.3 req/s
GET /features?bbox=…               58.8 kB    (20.8 kB gzipped)  105.1 req/s

Nine hundred times less data and forty-five times the throughput, for one query parameter.

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:

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.

Bar chart of throughput for unfiltered, paged and bounding-box-filtered requests.
Filters are not an optimisation to add later — they are why the service exists.

Step-by-step solution

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 — maxx before minx selects 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 nothing measurable: on 4,596 polygons, an indexed bounding-box selection took 0.43 ms and a brute-force intersects against a box took 0.37 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 time

3. 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))]

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"}

for parameter, column in FILTERABLE.items():
    value = request.query_params.get(parameter)
    if value is not None:
        subset = subset[subset[column] == value]

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: a page of 1,000 large polygons was 12.4 MB in measurement, and the server sustained 2.3 requests per second. A page of 50 was 857 kB at 32.5 requests per second.

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 a measured payload from 54.06 MB to 29.56 MB.

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.

Triage table of four invalid bounding boxes and their silent consequences.
Reject with a 400 and a message that says what was expected.

Code examples

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 subset

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 out

Measured on the full 4,596-feature layer, the three reductions individually:

as-is                             54.06 MB   (18.05 MB gzipped)
1 m coordinate precision          29.56 MB   ( 8.85 MB gzipped)
simplified to 0.01°               16.42 MB   ( 5.31 MB gzipped)

Precision alone removed 45% and changes nothing a consumer can detect. 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, rows

Parameter placeholders rather than string formatting, for the same reason as in any other SQL: a bounding box is client input.

Explanation

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.06 MB to 58.8 kB for a city-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 linear scan is 0.37 ms. Descending a tree, collecting candidates and testing them measured 0.43 ms — marginally slower.

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 unfiltered endpoint served 2.3 requests per second while returning 12.4 MB each time; ten such clients saturate the service.

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.

Two panels showing an envelope filter and an exact intersects test on the survivors.
That two-stage structure is what a spatial index does internally.

Edge cases or notes

  • bbox is 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 sindex at start-up, not per request.
  • X-Total-Count or numberMatched lets clients avoid pointless paging.
  • Round coordinates. Six decimal places is 11 cm and halves the payload.

FAQ

How much does a bounding box filter save?

Measured on a 4,596-feature dataset: the unfiltered response was 54.06 MB at 2.3 requests per second; a city-sized bounding box returned 58.8 kB at 105 requests per second.

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 a gain at small sizes: on 4,596 polygons an indexed selection took 0.43 ms against 0.37 ms for a brute-force test.

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 a 12 MB response at 2.3 requests per second, which any client can trigger.