How to Implement OGC API Features in Python

On this page

Problem statement

You need a feature service and you would rather not invent one. OGC API Features gives you the URL structure, the parameter names, the paging convention and the response shape — and a conformant service is consumable by QGIS, GDAL and several JavaScript clients without a line of custom code.

Two routes:

Write it. A conformant core service is about 120 lines of FastAPI — Example 1 is 122. You control the backend, the filters and the performance.
Deploy pygeoapi. A configuration file over an existing backend. Faster to stand up, and less flexible when the data does not fit its providers.

The measurements that matter are the same either way. On Natural Earth's 4,596 provinces, measured on the filtering page, a page of 50 features was 859 kB and one client got 47.6 requests per second; a page of 200 was 3.4 MB at 14.2; the whole layer is 54.29 MB. Whichever route you take, the limits decide whether the service works.

Quick answer

The five endpoints, minimally — an excerpt of the complete service in Example 1, which is what actually runs:

from fastapi import FastAPI, Query, Request, Response
import geopandas as gpd, json, datetime

app = FastAPI()
COLLECTIONS = {"provinces": gpd.read_file("provinces.gpkg")}
CRS84 = "http://www.opengis.net/def/crs/OGC/1.3/CRS84"
MAX_LIMIT = 200


@app.get("/")
def landing(request: Request):
    base = str(request.base_url).rstrip("/")
    return {"title": "Features service", "links": [
        {"href": f"{base}/", "rel": "self", "type": "application/json"},
        {"href": f"{base}/conformance", "rel": "conformance", "type": "application/json"},
        {"href": f"{base}/collections", "rel": "data", "type": "application/json"}]}


@app.get("/conformance")
def conformance():
    return {"conformsTo": [
        "http://www.opengis.net/spec/ogcapi-features-1/1.0/conf/core",
        "http://www.opengis.net/spec/ogcapi-features-1/1.0/conf/geojson"]}

The /collections, /collections/{id} and /collections/{id}/items endpoints follow the same pattern, and /items is where all the real work is.

max_items is not a tuning parameter; it is the difference between a service and an outage.

What you get — from The data used on this page

4596 provinces

Derived metadata, standard parameters, exact paging links

1. Build the collection description from the data

Everything in a collection's metadata can be derived, which means it cannot go stale:

def describe(name, gdf, base):
    minx, miny, maxx, maxy = gdf.total_bounds
    return {
        "id": name,
        "title": name.replace("_", " ").title(),
        "extent": {"spatial": {"bbox": [[float(minx), float(miny),
                                         float(maxx), float(maxy)]],
                               "crs": CRS84}},
        "itemType": "feature",
        "crs": [CRS84],
        "links": [{"href": f"{base}/collections/{name}/items",
                   "rel": "items", "type": "application/geo+json"}],
    }

The extent is the field clients use most: it lets a map zoom to the data before requesting any features.

2. Implement `/items` with the standard parameters

bbox, datetime, limit, offset and simple property equality. The bounding box is the one that matters: measured, it takes a response from the whole 54.29 MB layer to 116.8 kB for a London-sized window on a 4,596-feature dataset.

@app.get("/collections/{name}/items")
def items(request: Request, name: str,
          bbox: str | None = None, datetime_: str | None = Query(None, alias="datetime"),
          limit: int = Query(50, ge=1), offset: int = Query(0, ge=0)):
    limit = min(limit, MAX_LIMIT)          # clamp rather than refuse
    gdf = COLLECTIONS[name]

    if bbox:
        x0, y0, x1, y1 = (float(v) for v in bbox.split(","))
        gdf = gdf.cx[x0:x1, y0:y1]
    if datetime_:
        gdf = filter_datetime(gdf, datetime_)
    for key, value in request.query_params.items():
        if key in FILTERABLE:
            gdf = gdf[gdf[FILTERABLE[key]] == value]

    matched = len(gdf)
    page = gdf.iloc[offset:offset + limit]
    body = json.loads(page.to_json())
    body |= {"numberMatched": matched, "numberReturned": len(page),
             "timeStamp": datetime.datetime.now(datetime.UTC).isoformat(),
             "links": paging_links(request, offset, limit, matched)}
    return Response(json.dumps(body), media_type="application/geo+json")

3. Get the paging links right

Clients follow rel="next" rather than constructing URLs, so the links are the paging contract:

def paging_links(request, offset, limit, matched):
    self_url = str(request.url)
    base = str(request.url.remove_query_params(["offset"]))
    joiner = "&" if "?" in base else "?"
    links = [{"href": self_url, "rel": "self", "type": "application/geo+json"}]
    if offset + limit < matched:
        links.append({"href": f"{base}{joiner}offset={offset + limit}",
                      "rel": "next", "type": "application/geo+json"})
    if offset > 0:
        previous = max(0, offset - limit)
        links.append({"href": f"{base}{joiner}offset={previous}",
                      "rel": "prev", "type": "application/geo+json"})
    return links

4. Handle `datetime` in its three forms

The standard allows an instant, a closed interval and an open-ended one:

?datetime=2025-06-01T00:00:00Z
?datetime=2025-01-01/2025-12-31
?datetime=2025-01-01/..

def filter_datetime(gdf, spec, column="observed_at"):
    if "/" not in spec:
        return gdf[gdf[column] == spec]
    start, end = spec.split("/", 1)
    if start not in ("", ".."):
        gdf = gdf[gdf[column] >= start]
    if end not in ("", ".."):
        gdf = gdf[gdf[column] <= end]
    return gdf

5. Serve the single-feature endpoint

/collections/{name}/items/{featureId} returns a bare GeoJSON Feature — not a FeatureCollection. Clients rely on that difference.

@app.get("/collections/{name}/items/{feature_id}")
def item(name: str, feature_id: str):
    gdf = COLLECTIONS[name]
    match = gdf[gdf["id"].astype(str) == feature_id]
    if match.empty:
        raise HTTPException(404, f"no feature {feature_id!r} in {name!r}")
    feature = json.loads(match.to_json())["features"][0]
    return Response(json.dumps(feature), media_type="application/geo+json")

6. Decide between writing it and deploying pygeoapi

Write it when the backend is unusual, when the filters are domain-specific, or when performance matters enough that you want control of the query.

Deploy pygeoapi when the data sits in something it already supports — PostGIS, a file, an OGC service — and the requirement is a standards-compliant endpoint rather than a bespoke one. It is a YAML file and a container.

Either way, the limits are yours to set, and the measurements above are the reason to set them low.

The implementer already knows how big the dataset is; the client does not.

A complete core service and its conformance tests

The data used on this page

The collection is Natural Earth's 1:10m admin-1 states and provinces — 4,596 polygons with their name, country (admin) and type — saved as provinces.gpkg with an id column.

import io
import os
import zipfile

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}"))
    provinces = gpd.read_file("ne_provinces.shp")[["name", "admin", "type", "geometry"]]
    provinces.insert(0, "id", range(1, len(provinces) + 1))
    provinces.to_file("provinces.gpkg")
print(len(gpd.read_file("provinces.gpkg")), "provinces")

4596 provinces

Example 1 — the complete core service

Everything above, in one file:

%%writefile service.py
import datetime
import json

import geopandas as gpd
from fastapi import FastAPI, HTTPException, Query, Request, Response

app = FastAPI(title="Features service")
CRS84 = "http://www.opengis.net/def/crs/OGC/1.3/CRS84"
MAX_LIMIT = 200
FILTERABLE = {"admin": "admin", "type": "type"}

COLLECTIONS = {"provinces": gpd.read_file("provinces.gpkg").set_index("id", drop=False)}
for _gdf in COLLECTIONS.values():
    _gdf.sindex                      # build the spatial index at start-up


def base_of(request):
    return str(request.base_url).rstrip("/")


def describe(name, gdf, base):
    minx, miny, maxx, maxy = gdf.total_bounds
    return {"id": name, "title": name.replace("_", " ").title(),
            "extent": {"spatial": {"bbox": [[float(minx), float(miny), float(maxx), float(maxy)]],
                                   "crs": CRS84}},
            "itemType": "feature", "crs": [CRS84],
            "links": [{"href": f"{base}/collections/{name}/items", "rel": "items",
                       "type": "application/geo+json"},
                      {"href": f"{base}/collections/{name}", "rel": "self", "type": "application/json"}]}


def paging_links(request, offset, limit, matched):
    self_url = str(request.url)
    base = str(request.url.remove_query_params(["offset"]))
    joiner = "&" if "?" in base else "?"
    links = [{"href": self_url, "rel": "self", "type": "application/geo+json"}]
    if offset + limit < matched:
        links.append({"href": f"{base}{joiner}offset={offset + limit}",
                      "rel": "next", "type": "application/geo+json"})
    if offset > 0:
        links.append({"href": f"{base}{joiner}offset={max(0, offset - limit)}",
                      "rel": "prev", "type": "application/geo+json"})
    return links


def filter_datetime(gdf, spec, column="observed_at"):
    if "/" not in spec:
        return gdf[gdf[column] == spec]
    start, end = spec.split("/", 1)
    if start not in ("", ".."):
        gdf = gdf[gdf[column] >= start]
    if end not in ("", ".."):
        gdf = gdf[gdf[column] <= end]
    return gdf


@app.get("/")
def landing(request: Request):
    base = base_of(request)
    return {"title": "Features service", "links": [
        {"href": f"{base}/", "rel": "self", "type": "application/json"},
        {"href": f"{base}/conformance", "rel": "conformance", "type": "application/json"},
        {"href": f"{base}/collections", "rel": "data", "type": "application/json"}]}


@app.get("/conformance")
def conformance():
    return {"conformsTo": [
        "http://www.opengis.net/spec/ogcapi-features-1/1.0/conf/core",
        "http://www.opengis.net/spec/ogcapi-features-1/1.0/conf/geojson"]}


@app.get("/collections")
def collections(request: Request):
    base = base_of(request)
    return {"collections": [describe(name, gdf, base) for name, gdf in COLLECTIONS.items()],
            "links": [{"href": f"{base}/collections", "rel": "self", "type": "application/json"}]}


@app.get("/collections/{name}")
def collection(request: Request, name: str):
    if name not in COLLECTIONS:
        raise HTTPException(404, f"no collection {name!r}")
    return describe(name, COLLECTIONS[name], base_of(request))


@app.get("/collections/{name}/items")
def items(request: Request, name: str,
          bbox: str | None = None, datetime_: str | None = Query(None, alias="datetime"),
          limit: int = Query(50, ge=1), offset: int = Query(0, ge=0)):
    if name not in COLLECTIONS:
        raise HTTPException(404, f"no collection {name!r}")
    limit = min(limit, MAX_LIMIT)    # clamp, don't refuse: GDAL asks for 1,000
    gdf = COLLECTIONS[name]
    if bbox:
        x0, y0, x1, y1 = (float(v) for v in bbox.split(","))
        gdf = gdf.cx[x0:x1, y0:y1]
    if datetime_:
        gdf = filter_datetime(gdf, datetime_)
    for key, value in request.query_params.items():
        if key in FILTERABLE:
            gdf = gdf[gdf[FILTERABLE[key]] == value]

    matched = len(gdf)
    page = gdf.iloc[offset:offset + limit]
    body = json.loads(page.to_json())
    body |= {"numberMatched": matched, "numberReturned": len(page),
             "timeStamp": datetime.datetime.now(datetime.UTC).isoformat(),
             "links": paging_links(request, offset, limit, matched)}
    return Response(json.dumps(body), media_type="application/geo+json")


@app.get("/collections/{name}/items/{feature_id}")
def item(name: str, feature_id: str):
    if name not in COLLECTIONS:
        raise HTTPException(404, f"no collection {name!r}")
    gdf = COLLECTIONS[name]
    match = gdf[gdf["id"].astype(str) == feature_id]
    if match.empty:
        raise HTTPException(404, f"no feature {feature_id!r} in {name!r}")
    feature = json.loads(match.to_json())["features"][0]
    return Response(json.dumps(feature), media_type="application/geo+json")

A conformant service is one a generic client can use with nothing but its URL. GDAL's OAPIF driver is such a client, so the real test is to start the service and let GDAL find the collection, page through a filtered query, and count what it got:

import socket
import subprocess
import sys
import time

import httpx
from osgeo import ogr

gdal_ok = ogr.UseExceptions()
with socket.socket() as probe:
    probe.bind(("127.0.0.1", 0))
    port = probe.getsockname()[1]
server = subprocess.Popen([sys.executable, "-m", "uvicorn", "service: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 + "/conformance", timeout=1); break
    except httpx.TransportError:
        time.sleep(0.1)

try:
    source = ogr.Open(f"OAPIF:{BASE}")
    layer = source.GetLayer(0)
    print("GDAL sees", source.GetLayerCount(), "collection:", layer.GetName())
    layer.SetSpatialFilterRect(-0.51, 51.28, 0.33, 51.69)          # London
    names = sorted(f.GetField("name") for f in layer)
    print(len(names), "features in the London box, first three:", names[:3])
    layer.SetSpatialFilter(None)
    layer.SetAttributeFilter("admin = 'Brazil'")
    print(sum(1 for _ in layer), "features for admin = 'Brazil'")
    source = None
finally:
    server.terminate(); server.wait()

GDAL sees 1 collection: provinces
39 features in the London box, first three: ['Barking and Dagenham', 'Barnet', 'Bexley']
27 features for admin = 'Brazil'

The first version of this service answered limit above 200 with a 422, and GDAL stopped at its first request: its OAPIF driver asks for limit=1000. Clamping the limit to the maximum and saying how many were returned is what lets a generic client page through — numberReturned and the next link tell it the rest.

Example 2 — testing conformance in CI

%%writefile test_oapif.py
import pytest
from fastapi.testclient import TestClient


CORE = "http://www.opengis.net/spec/ogcapi-features-1/1.0/conf/core"


@pytest.fixture(scope="module")
def client():
    from service import app
    return TestClient(app)


def test_landing_page_links(client):
    body = client.get("/").json()
    rels = {link["rel"] for link in body["links"]}
    assert {"self", "conformance", "data"} <= rels


def test_conformance_declares_core(client):
    assert CORE in client.get("/conformance").json()["conformsTo"]


def test_items_shape(client):
    body = client.get("/collections/provinces/items?limit=5").json()
    assert body["type"] == "FeatureCollection"
    assert body["numberReturned"] == len(body["features"]) == 5
    assert "numberMatched" in body


def test_paging_follows_links(client):
    body = client.get("/collections/provinces/items?limit=5").json()
    next_link = next(l for l in body["links"] if l["rel"] == "next")
    second = client.get(next_link["href"]).json()
    first_ids = {f["id"] for f in body["features"]}
    second_ids = {f["id"] for f in second["features"]}
    assert not (first_ids & second_ids), "pages overlap"


def test_limit_is_capped(client):
    body = client.get("/collections/provinces/items?limit=100000").json()
    assert body["numberReturned"] == 200 < body["numberMatched"]


def test_single_feature_is_a_feature(client):
    items = client.get("/collections/provinces/items?limit=1").json()
    feature_id = items["features"][0]["id"]
    body = client.get(f"/collections/provinces/items/{feature_id}").json()
    assert body["type"] == "Feature"

!{sys.executable} -m pytest -q -p no:cacheprovider test_oapif.py 2>&1 | tail -2

......                                                                   [100%]
6 passed in 0.58s

test_paging_follows_links is the one that catches the most: it follows the link the service produced rather than a URL the test constructed, which is exactly what a real client does.

Example 3 — a pygeoapi configuration for the same data

server:
  bind: {host: 0.0.0.0, port: 5000}
  url: https://example.org/oapif
  limits: {default_items: 50, max_items: 200}

resources:
  provinces:
    type: collection
    title: Administrative provinces
    extents:
      spatial: {bbox: [-180, -90, 180, 90], crs: "…/CRS84"}
    providers:
      - type: feature
        name: PostgreSQL
        data:
          host: db
          dbname: gis
          user: reader
        id_field: id
        table: provinces
        geom_field: geom

docker run -p 5000:80 -v $(pwd)/config.yml:/pygeoapi/local.config.yml \
  geopython/pygeoapi:latest

Note max_items: 200. Set it explicitly rather than trusting whatever the deployment ships with; the measurements at the top are the reason. This configuration was not run here — it needs a PostgreSQL database — but it describes the same collection the Python service serves.

What a conformant service gets for free

Why a conformant service pays for itself immediately

QGIS, GDAL's OAPIF driver and several JavaScript clients consume conformant services with no configuration beyond a URL. That removes the wrapper each of those would otherwise need.

It also removes design decisions. The bounding-box parameter is bbox, paging is limit and offset, the count is numberMatched, and nobody has to have an opinion.

Why the paging links are the part to get exactly right

A client that follows links is coupled only to the standard; one that builds URLs is coupled to your URL structure. The link-driven design lets a service move to cursor paging, change its base path, or sit behind a CDN with different URLs without breaking a single client.

The bug to avoid is a next link that drops the other query parameters. A next link without the original bbox starts paging the whole dataset — and it looks like a paging bug rather than a filter bug.

Why `numberMatched` matters more than it looks

It lets a client decide not to page. A query matching 4,596 features at a page size of 50 is ninety-two requests; a client that knows the number can narrow the filter instead.

It is also what makes progress reporting possible, and it is the field most often omitted from home-made services, because the implementer already knows how large the dataset is.

Why the limits are the most important configuration

Measured on Natural Earth's provinces: 50 features per page was 859 kB at 47.6 requests per second, and 200 was 3.4 MB at 14.2.

A service that allows limit=10000 has a worst case of the whole 54.29 MB layer — sixteen times the 200-feature page — reachable by any client with a URL bar. max_items is not a tuning parameter; it is the difference between a service and an outage.

It presents as a paging problem and is a URL-construction problem.

Single features, carried parameters and clamped limits

/items/{id} returns a Feature, not a FeatureCollection.
The next link must carry every other parameter, especially bbox.
CRS84 is longitude-first. Do not label it EPSG:4326 and hope.
limit should clamp, never be unbounded — and not 422: GDAL's OAPIF driver asks for 1,000 and gives up on the refusal.
datetime has three forms including open-ended intervals with ...
f=json as a query parameter is far friendlier for browser testing than an Accept header.
Build the spatial index at start-up, not per request.
Set pygeoapi's max_items explicitly before deploying.

Internal links

OGC API Features explained — the standard itself
Serving spatial data explained: files, features and tiles — whether a feature API is the right shape
How to paginate a large feature API in Python — the paging details
How to add bounding box and attribute filters to a spatial API — the filters
How to test a spatial API with pytest and httpx — the test suite
How to build a GeoJSON API with FastAPI — a simpler service
Spatial web services explained — WFS and the older standards
How to set cache headers on a spatial API and tile service — making it fast

FAQ

Should I write an OGC API Features service or deploy pygeoapi?

Deploy pygeoapi when the data sits in a backend it supports and you want a standard endpoint. Write it when the backend is unusual, the filters are domain-specific, or you need control of the queries.

How much code is a conformant service?

About 120 lines of FastAPI for the core plus GeoJSON conformance classes — the one in Example 1 is 122 — five endpoints, the paging links and the standard parameters.

What is the most important limit to set?

max_items. A page of 50 features was 859 kB at 47.6 requests per second; a page of 200 was 3.4 MB at 14.2, and the whole layer is 54.29 MB.

What is the commonest implementation bug?

A next link that drops the other query parameters, so paging silently starts returning the whole dataset instead of the filtered result.

Does the single-feature endpoint return a FeatureCollection?

No — a bare Feature. Clients depend on the difference.

How do I test conformance?

Follow the links the service produces rather than constructing URLs, and assert the response shapes: landing links, conformance classes, numberMatched, non-overlapping pages and a capped limit.