⚙️ How-To Web maps & tiles Pipelines & automation Publishing Web Maps · step 15 of 15 Updated

How to Publish a GeoJSON API from a GeoDataFrame with FastAPI

Problem statement

Sometimes a static file is not enough: the data changes, the client needs a filtered subset, or the dataset is too large to send whole. A small API answers "give me the features in this box matching this filter".

Four things decide whether it survives contact with real use:

  • a bounding-box requirement, or one client asks for everything
  • a result cap, or one query returns 59,391 features
  • a spatial index, or every request scans the whole layer
  • EPSG:4326 output, or nothing renders

Measured, 59,391 building polygons as GeoJSON is 27.04 MB — 2.65 MB gzipped. An unbounded endpoint over that dataset is a denial-of-service button.

Quick answer 

import geopandas as gpd
from fastapi import FastAPI, HTTPException, Query
from fastapi.responses import JSONResponse

app = FastAPI()
DATA = gpd.read_parquet("buildings.parquet").to_crs(4326)
SINDEX = DATA.sindex
MAX_FEATURES = 5000


@app.get("/features")
def features(
    bbox: str = Query(..., description="west,south,east,north"),
    limit: int = Query(1000, le=MAX_FEATURES),
):
    try:
        west, south, east, north = (float(v) for v in bbox.split(","))
    except ValueError:
        raise HTTPException(400, "bbox must be west,south,east,north")
    if east <= west or north <= south:
        raise HTTPException(400, "bbox is empty or inverted")

    idx = list(SINDEX.intersection((west, south, east, north)))
    subset = DATA.iloc[idx]
    truncated = len(subset) > limit

    return JSONResponse({
        **subset.head(limit).__geo_interface__,
        "numberMatched": len(subset),
        "numberReturned": min(len(subset), limit),
        "truncated": truncated,
    })

Requiring the bounding box rather than defaulting it is the decision that keeps the service up.

A request validated for bbox and limit, filtered through a spatial index, capped, and returned as GeoJSON with match counts.
Validate, index, cap, report. Every step exists because of a failure mode.

Step-by-step solution

1. Require a bounding box

An endpoint that returns everything when no filter is given will be called that way, by a crawler if not by a user. Make bbox required and validate it — inverted or zero-area boxes are common client bugs and should be a 400, not an empty result.

2. Use the spatial index 

idx = list(SINDEX.intersection(bounds))
subset = DATA.iloc[idx]

gdf.cx[...] rebuilds its mask each call, which is linear in the feature count. sindex.intersection is an R-tree lookup. Over 59,391 features that is the difference between milliseconds and tens of milliseconds per request — and it compounds under load.

Note that the R-tree query is on bounding boxes, so it returns candidates. If exact containment matters, follow it with a precise test on the candidates only.

3. Cap the result and say that you did 

"numberMatched": len(subset),
"numberReturned": min(len(subset), limit),

A silently truncated response is worse than an error: the client draws a partial map and believes it is complete. Reporting both counts — the OGC API Features convention — lets the client detect it.

4. Return EPSG:4326

GeoJSON is defined in WGS 84 longitude and latitude, and clients assume it. Reproject once at load, not per request.

5. Load the data once, at startup

Reading a GeoParquet file takes 0.06 seconds; reading GeoJSON takes 0.49. Neither belongs inside a request handler. Load at startup and build the spatial index there too.

An endpoint without a required bbox returning 59,391 features and 27 MB, against a bounded and capped endpoint returning a few hundred.
An unbounded endpoint over a large layer is a denial-of-service button with a friendly URL.

Code examples

Example 1 — a service with the essentials 

import geopandas as gpd
from fastapi import FastAPI, HTTPException, Query
from fastapi.responses import JSONResponse

app = FastAPI(title="Features API")

STATE = {}
MAX_FEATURES = 5000
MAX_BBOX_DEGREES = 1.0


@app.on_event("startup")
def load():
    gdf = gpd.read_parquet("buildings.parquet").to_crs(4326)
    gdf = gdf[gdf.geometry.notna() & ~gdf.geometry.is_empty]
    STATE["data"] = gdf
    STATE["sindex"] = gdf.sindex
    STATE["fields"] = [c for c in gdf.columns if c != "geometry"]
    print(f"  loaded {len(gdf):,} features, fields {STATE['fields']}")


def parse_bbox(bbox: str):
    try:
        west, south, east, north = (float(v) for v in bbox.split(","))
    except ValueError:
        raise HTTPException(400, "bbox must be four numbers: "
                                 "west,south,east,north")
    if east <= west or north <= south:
        raise HTTPException(400, "bbox is empty or inverted")
    if (east - west) > MAX_BBOX_DEGREES or (north - south) > MAX_BBOX_DEGREES:
        raise HTTPException(400, f"bbox exceeds {MAX_BBOX_DEGREES} degrees; "
                                 "request a smaller area")
    return west, south, east, north


@app.get("/collections/features/items")
def items(bbox: str = Query(...), limit: int = Query(1000, ge=1, le=MAX_FEATURES),
          properties: str | None = None, where: str | None = None):
    bounds = parse_bbox(bbox)
    gdf = STATE["data"]

    idx = list(STATE["sindex"].intersection(bounds))
    subset = gdf.iloc[idx]

    if where:
        try:
            subset = subset.query(where)
        except Exception as exc:
            raise HTTPException(400, f"invalid filter: {exc}")

    matched = len(subset)
    subset = subset.head(limit)

    if properties:
        wanted = [p for p in properties.split(",") if p in gdf.columns]
        subset = subset[wanted + ["geometry"]]

    payload = subset.__geo_interface__
    payload["numberMatched"] = matched
    payload["numberReturned"] = len(subset)
    payload["links"] = [{"rel": "self", "href": f"/collections/features/items"
                                                f"?bbox={bbox}&limit={limit}"}]
    return JSONResponse(payload, headers={"Cache-Control": "public, max-age=60"})

Capping the bounding box size as well as the result count is the second line of defence. A client asking for the whole world with limit=1000 still forces an index query over every feature.

Example 2 — cutting the payload 

from fastapi import Response
import gzip
import json


def compact_geojson(gdf, decimals=6, drop_columns=()):
    """Round coordinates and drop attributes before serialising."""
    out = gdf.drop(columns=[c for c in drop_columns if c in gdf.columns])
    out = out.copy()
    out["geometry"] = out.geometry.set_precision(10 ** -decimals)
    return out.__geo_interface__


@app.get("/features.geojson")
def features_compact(bbox: str, limit: int = 1000, decimals: int = 6):
    bounds = parse_bbox(bbox)
    subset = STATE["data"].iloc[
        list(STATE["sindex"].intersection(bounds))].head(limit)

    payload = json.dumps(compact_geojson(subset, decimals=decimals))
    compressed = gzip.compress(payload.encode(), 6)
    print(f"  {len(subset):,} features, {len(payload) / 1e6:.2f} MB, "
          f"{len(compressed) / 1e6:.2f} MB gzipped")
    return Response(compressed, media_type="application/geo+json",
                    headers={"Content-Encoding": "gzip"})

Six decimal places is about 11 cm. Measured on a full dataset, going from seven decimals to four took a GeoJSON from 2.83 MB to 1.11 MB gzipped — the compressed saving is larger than the raw one, because shorter numbers are more repetitive.

Example 3 — a summary endpoint for large areas 

from fastapi import Query
import numpy as np


@app.get("/summary")
def summary(bbox: str = Query(...), cells: int = Query(20, ge=2, le=100)):
    """Counts on a grid, so a large area returns something small."""
    west, south, east, north = parse_bbox_unlimited(bbox)
    gdf = STATE["data"]
    idx = list(STATE["sindex"].intersection((west, south, east, north)))
    subset = gdf.iloc[idx]

    if subset.empty:
        return {"cells": [], "total": 0}

    points = subset.geometry.representative_point()
    col = np.clip(((points.x - west) / (east - west) * cells).astype(int),
                  0, cells - 1)
    row = np.clip(((north - points.y) / (north - south) * cells).astype(int),
                  0, cells - 1)
    counts = np.zeros(cells * cells, int)
    np.add.at(counts, row * cells + col, 1)

    step_x = (east - west) / cells
    step_y = (north - south) / cells
    features = []
    for i, count in enumerate(counts):
        if not count:
            continue
        r, c = divmod(i, cells)
        features.append({
            "type": "Feature",
            "properties": {"count": int(count)},
            "geometry": {"type": "Polygon", "coordinates": [[
                [west + c * step_x, north - (r + 1) * step_y],
                [west + (c + 1) * step_x, north - (r + 1) * step_y],
                [west + (c + 1) * step_x, north - r * step_y],
                [west + c * step_x, north - r * step_y],
                [west + c * step_x, north - (r + 1) * step_y]]]},
        })

    return {"type": "FeatureCollection", "features": features,
            "total": int(counts.sum())}

Giving clients a cheap way to ask about a large area is what stops them asking for it feature by feature. A grid of counts over a whole city is a few kilobytes and answers most "where is the data dense" questions.

Explanation

Why the bounding box must be required

Optional filters get omitted. A crawler, a misconfigured client or a developer testing in a browser will call the bare endpoint, and if that returns everything, the service serialises 27 MB and blocks a worker for seconds.

Making it required turns that into a 422 from FastAPI's validation, before any work happens. Capping its size closes the remaining hole, where a client asks for the whole world.

Why the index matters under load

A single gdf.cx[...] over 59,391 features is fast enough not to notice. At a hundred requests a second it is the whole CPU budget.

The R-tree turns each query into a tree descent proportional to the number of results rather than the size of the dataset. It is built once at startup and shared, since it is read-only.

That also means the service is effectively stateless per request, which is what allows it to scale by running more processes.

Why to report numberMatched

A truncated GeoJSON is indistinguishable from a complete one. The client draws it, the user reads it as the full picture, and nothing anywhere indicates otherwise.

The OGC API Features convention of returning both numberMatched and numberReturned makes truncation detectable. A client that cares can page; one that does not at least has the information.

Why an API instead of tiles

Tiles are better for rendering: they are pre-generated, cacheable, and the client only loads the current view.

An API is better when the client needs the data rather than a picture — attribute queries, exact geometry, arbitrary filters, or a subset defined by something other than the viewport. It is also better when the data changes faster than a pyramid can be regenerated.

Many systems need both: tiles for the map, an API for the query panel and the export button.

A features endpoint requiring a small bounding box against a summary endpoint returning grid counts for a large one.
Most "where is the data dense" questions do not need geometry at all.

Edge cases or notes

  • Require the bounding box and cap its size.
  • Cap the result count and report numberMatched.
  • Build the spatial index at startup, not per request.
  • Return EPSG:4326; reproject once at load.
  • Round coordinates to six decimals — about 11 cm.
  • Gzip the response; GeoJSON compresses about tenfold.
  • Validate filter expressions rather than passing them to query unchecked.
  • Offer a summary endpoint so large-area questions have a cheap answer.

FAQ

How do I serve a GeoDataFrame as an API?

Load it once at startup, reproject to EPSG:4326, build a spatial index, and expose an endpoint that requires a bounding box, caps the result count and reports how many matched.

Why require a bounding box?

Because an optional filter will be omitted. Serialising 59,391 features is 27 MB and blocks a worker for seconds.

How do I stop a client asking for too much?

Cap the result count, cap the bounding-box size, and offer a summary endpoint that answers large-area questions cheaply.

Why report numberMatched?

A truncated response looks identical to a complete one. Reporting both counts is the OGC API Features convention and the only way a client can detect truncation.

Should I use .cx or the spatial index?

The index. .cx rebuilds its mask each call, which is linear in the dataset size and becomes the whole CPU budget under load.

How much precision should the API return?

Six decimal places, about 11 cm. Beyond that is invisible on any map and costs real bytes.

API or vector tiles?

Tiles for rendering; an API when the client needs the data itself — attribute queries, exact geometry, or arbitrary filters. Large systems usually have both.