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 every building in the city
• a spatial index, or every request scans the whole layer
• EPSG:4326 output, or nothing renders

Measured on Leiden's 58,879 OpenStreetMap buildings, the whole layer as GeoJSON is 20.80 MB — 2.98 MB gzipped (Example 2). 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,
    })

from fastapi.testclient import TestClient

client = TestClient(app)
response = client.get("/features", params={"bbox": "4.48,52.155,4.50,52.165"})
body = response.json()
print({k: body[k] for k in ("numberMatched", "numberReturned", "truncated")},
      f"{len(response.content) / 1e6:.2f} MB")
print(client.get("/features").status_code, client.get("/features").json()["detail"][0]["msg"])
print(client.get("/features", params={"bbox": "4.50,52.16,4.48,52.17"}).json())

{'numberMatched': 6278, 'numberReturned': 1000, 'truncated': True} 0.42 MB
422 Field required
{'detail': 'bbox is empty or inverted'}

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

Require a bbox, use the index, cap and report

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

import timeit

bounds = (4.48, 52.155, 4.50, 52.165)
via_index = timeit.timeit(lambda: DATA.iloc[list(SINDEX.intersection(bounds))], number=50) / 50
via_cx = timeit.timeit(lambda: DATA.cx[4.48:4.50, 52.155:52.165], number=50) / 50
print(f"sindex {via_index * 1000:.2f} ms, .cx {via_cx * 1000:.2f} ms per query")

sindex 0.92 ms, .cx 2.65 ms per query

gdf.cx[...] rebuilds its mask each call, which is linear in the feature count. sindex.intersection is an R-tree lookup. Over 58,879 features the difference is the one measured above — 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

import time

DATA.to_file("buildings.geojson", driver="GeoJSON")
for path, read in (("buildings.parquet", gpd.read_parquet), ("buildings.geojson", gpd.read_file)):
    started = time.perf_counter()
    read(path)
    print(f"{path:18} {time.perf_counter() - started:.2f} s")

buildings.parquet  0.08 s
buildings.geojson  0.40 s

Neither belongs inside a request handler. Load at startup and build the spatial index there too.

A service with the essentials and a lighter payload

The data used on this page

The layer is every building in Leiden in OpenStreetMap — about 59,000 footprints with their building type and, where mapped, a name, height and number of levels — fetched with OSMnx through an Overpass mirror and saved as GeoParquet in EPSG:4326. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

import io
import os
import requests
import zipfile

import geopandas as gpd
import osmnx as ox

# OpenStreetMap as Overpass returned it on 10 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-geojson-api-fastapi-python-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"   # overpass-api.de mirror
ox.settings.overpass_rate_limit = False
if not os.path.exists("buildings.parquet"):
    found = ox.features_from_place("Leiden, Netherlands", {"building": True}).reset_index()
    found = found[found.geom_type.isin(["Polygon", "MultiPolygon"])]
    keep = [c for c in ("osmid", "building", "name", "height", "building:levels") if c in found.columns]
    found[keep + ["geometry"]].rename(columns={"building:levels": "levels"}).to_parquet("buildings.parquet")
buildings = gpd.read_parquet("buildings.parquet")
print(f"{len(buildings):,} OSM buildings in Leiden, EPSG:{buildings.crs.to_epsg()}; "
      f"columns {[c for c in buildings.columns if c != 'geometry']}")

58,879 OSM buildings in Leiden, EPSG:4326; columns ['building', 'name', 'height', 'levels']

Example 1 — a service with the essentials

from contextlib import asynccontextmanager

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

STATE = {}
MAX_FEATURES = 5000
MAX_BBOX_DEGREES = 0.1


@asynccontextmanager
async def lifespan(app):
    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']}")
    yield


app = FastAPI(title="Features API", lifespan=lifespan)


def parse_bbox_unlimited(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")
    return west, south, east, north


def parse_bbox(bbox: str):
    west, south, east, north = parse_bbox_unlimited(bbox)
    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"})

with TestClient(app) as client:                     # runs the lifespan: loads the data once
    ok = client.get("/collections/features/items",
                    params={"bbox": "4.48,52.155,4.50,52.165", "where": "building == 'house'",
                            "properties": "building,name", "limit": 100}).json()
    print(ok["numberMatched"], "houses matched,", ok["numberReturned"], "returned")
    print(client.get("/collections/features/items", params={"bbox": "4,52,5,53"}).json())

  loaded 58,879 features, fields ['building', 'name', 'height', 'levels']
2721 houses matched, 100 returned
{'detail': 'bbox exceeds 0.1 degrees; request a smaller area'}

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

The whole layer, at full precision and then rounded:

everything = STATE["data"]
for label, text in (("to_json()", everything.to_json()),
                    ("__geo_interface__", json.dumps(everything.__geo_interface__))):
    print(f"{label:18} {len(text) / 1e6:6.2f} MB, {len(gzip.compress(text.encode(), 6)) / 1e6:5.2f} MB gzipped")
for decimals in (6, 5, 4):
    text = json.dumps(compact_geojson(everything, decimals=decimals))
    emptied = everything.geometry.set_precision(10 ** -decimals).is_empty.sum()
    print(f"{decimals} decimals: {len(text) / 1e6:6.2f} MB, "
          f"{len(gzip.compress(text.encode(), 6)) / 1e6:5.2f} MB gzipped, {emptied:,} buildings collapsed to nothing")

to_json()           20.80 MB,  2.98 MB gzipped
__geo_interface__   24.07 MB,  3.71 MB gzipped
6 decimals:  22.96 MB,  3.15 MB gzipped, 0 buildings collapsed to nothing
5 decimals:  21.37 MB,  2.44 MB gzipped, 46 buildings collapsed to nothing
4 decimals:  13.16 MB,  0.92 MB gzipped, 32,384 buildings collapsed to nothing

The two unrounded lines differ only in serialiser: __geo_interface__ adds a bounding box to every feature, and to_json() leaves it out. The rounded lines use the __geo_interface__ route, so compare them with the second line. OpenStreetMap stores coordinates to seven decimals — about a centimetre — so rounding to seven would change nothing. Six decimal places is about 11 cm and cost nothing. Five collapsed 46 small buildings. Four — about 11 m — is not a size optimisation at all: the file shrank because 32,384 of the 58,879 buildings snapped to a grid coarser than themselves and vanished. Check for empty geometries whenever you round.

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())}

with TestClient(app) as client:
    grid = client.get("/summary", params={"bbox": "4.43,52.12,4.54,52.19", "cells": 20})
    body = grid.json()
    print(f"{body['total']:,} buildings in {len(body['features'])} occupied cells, "
          f"{len(grid.content) / 1e3:.1f} kB; busiest cell {max(f['properties']['count'] for f in body['features']):,}")

  loaded 58,879 features, fields ['building', 'name', 'height', 'levels']
58,879 buildings in 174 occupied cells, 35.2 kB; busiest cell 1,396

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.

Why the bounding box must be required

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 the whole layer — tens of megabytes for one city's buildings — and blocks a worker while it does.

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 58,879 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.

Caps, startup indexes and WGS 84

• 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 — and never to four for buildings; it emptied 55% of them.
• Gzip the response; the Leiden buildings compressed from 20.80 MB to 2.98 MB.
• Validate filter expressions rather than passing them to query unchecked.
• Offer a summary endpoint so large-area questions have a cheap answer.

Internal links

• How to prepare a GeoDataFrame for the web — the output format rules
• My GeoJSON is too big for the browser — why the cap matters
• How to use the spatial index directly in GeoPandas (sindex) — the index used here
• Spatial web services explained: WFS, WMS and OGC API Features — the standards this imitates
• How to build vector tiles from a GeoDataFrame in Python — the rendering alternative
• How to write partitioned GeoParquet and query it with filters — a larger-than-memory backend
• How to add interactivity and popups to a web map — the client that consumes this
• How to handle credentials and secrets in an automated GIS job — securing the service

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 all 58,879 Leiden buildings is 20.80 MB of GeoJSON, built while a worker waits.

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 — and rounding much harder deletes features: at four decimals, 32,384 of Leiden's 58,879 buildings collapsed to nothing.

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.