This guide is part of Mastering Geospatial Data Ingestion in Python.

Ingesting WFS and OGC API — Features Services in Python

The Problem: A Standards-Compliant Endpoint Is Not a Predictable One

Web Feature Service and its successor, OGC API — Features, are the two interfaces through which most public authorities publish vector data. Both are standardised, both are well documented, and neither behaves the same way twice across publishers. A pipeline that reads one national cadastre successfully will fail against the next one for reasons that are entirely legal under the specification: a different default coordinate order, a server-side feature cap the capabilities document does not mention, an output format list that includes GeoJSON but returns GML anyway.

Naive approaches fail in three specific ways. Reading the whole layer in one request works in development against a test extent and times out against the real one. Trusting the declared CRS produces silently transposed coordinates for half the services in Europe. And paging with a numeric startIndex against a service whose ordering is unstable returns duplicated and missing features in the same run — the pattern described in parsing GeoJSON and shapefile APIs, and worse here because the specification does not require a stable sort.

The ingestion loop this page builds negotiates capabilities first, pins the coordinate order explicitly, pages with a bounded request size, and records enough state to resume.

Prerequisites and Environment

pip install "owslib>=0.31" "geopandas>=1.0" "pyogrio>=0.9" "requests>=2.32" "shapely>=2.0"

OWSLib handles the WFS capabilities negotiation; pyogrio reads the GeoJSON and GML payloads through GDAL without a temporary file for the JSON case. Verify the GDAL build has the GML driver before relying on it, since minimal wheels sometimes omit it:

import pyogrio

drivers = pyogrio.list_drivers()
assert "GML" in drivers, "GDAL build lacks the GML driver — install libgdal with the full driver set"
assert "GeoJSON" in drivers

Version and Compatibility Matrix

Step 1 — Negotiate Capabilities Before Requesting Anything

The capabilities document tells you which layers exist, which output formats they support, which coordinate reference systems they can serve, and — critically — whether the service caps the number of features per request.

import logging
from owslib.wfs import WebFeatureService

logger = logging.getLogger(__name__)


def describe_wfs_layer(url: str, layer: str, version: str = "2.0.0") -> dict:
    """Return the facts an ingestion loop needs before its first GetFeature."""
    wfs = WebFeatureService(url=url, version=version, timeout=60)
    if layer not in wfs.contents:
        raise KeyError(f"layer {layer!r} not in service; available: {sorted(wfs.contents)[:10]}")

    meta = wfs.contents[layer]
    operation = next(op for op in wfs.operations if op.name == "GetFeature")
    formats = operation.parameters.get("outputFormat", {}).get("values", [])

    facts = {
        "title": meta.title,
        "crs_options": [str(c) for c in meta.crsOptions],
        "bbox_wgs84": meta.boundingBoxWGS84,
        "output_formats": formats,
        "prefers_geojson": any("json" in f.lower() for f in formats),
    }
    logger.info("layer %s: %d CRS options, geojson=%s", layer, len(facts["crs_options"]),
                facts["prefers_geojson"])
    return facts

The prefers_geojson flag decides the rest of the pipeline. A service that offers GeoJSON removes the GML parsing step entirely; one that does not requires the approach in converting GML responses to GeoDataFrames.

Step 2 — Count Before You Download

resultType=hits costs one cheap request and changes what the ingestion does. Under a few thousand features, page once and be done; over a few hundred thousand, split by bounding box before paging at all.

import requests
from xml.etree import ElementTree

WFS_NS = {"wfs": "http://www.opengis.net/wfs/2.0"}


def count_features(url: str, layer: str, bbox: tuple[float, float, float, float] | None = None) -> int:
    """Ask the service how many features match, without transferring any of them."""
    params = {
        "service": "WFS", "version": "2.0.0", "request": "GetFeature",
        "typeNames": layer, "resultType": "hits",
    }
    if bbox:
        params["bbox"] = ",".join(str(v) for v in bbox) + ",urn:ogc:def:crs:EPSG::4326"

    response = requests.get(url, params=params, timeout=60)
    response.raise_for_status()
    root = ElementTree.fromstring(response.content)
    matched = root.attrib.get("numberMatched", "unknown")
    return -1 if matched == "unknown" else int(matched)

A service that returns unknown is telling you it cannot count without scanning. Treat that as a signal to partition spatially rather than as a reason to skip the check.

Step 3 — Pin the Coordinate Order Explicitly

This is the single most common source of silently wrong output from a standards-compliant service. WFS 1.1.0 and 2.0.0 follow the axis order defined by the EPSG authority, and for EPSG:4326 that order is latitude first. A client that assumes longitude first plots Edinburgh in the Indian Ocean.

from pyproj import CRS
import geopandas as gpd


def normalize_axis_order(gdf: gpd.GeoDataFrame, requested_crs: str) -> gpd.GeoDataFrame:
    """Return a frame in longitude-latitude order, whatever the service sent."""
    crs = CRS.from_user_input(requested_crs)
    axis_names = [ax.abbrev.lower() for ax in crs.axis_info]

    if axis_names[:2] == ["lat", "lon"]:
        # Authority order — swap into the x=lon, y=lat convention the rest of the pipeline assumes.
        flipped = gdf.geometry.map(lambda geom: shapely.ops.transform(lambda x, y: (y, x), geom))
        gdf = gdf.set_geometry(gpd.GeoSeries(flipped, crs=4326))
    return gdf.set_crs(4326, allow_override=True)

The reliable alternative is to avoid the ambiguity: request urn:ogc:def:crs:OGC:1.3:CRS84, which is WGS 84 with longitude first by definition. Every OGC API — Features service supports it, and most modern WFS deployments do too. Where it is unavailable, verify the result against a known point before trusting a full ingest — a check that costs one assertion and prevents an entire dataset landing transposed, in the same spirit as the CRS invariants set out in CRS normalization across mixed datasets.

Step 4 — Page With Bounded Requests and Persisted State

from pathlib import Path
import json

import geopandas as gpd
import pandas as pd
import requests


def page_ogc_items(items_url: str, page_size: int = 1000,
                   state_path: Path | None = None) -> gpd.GeoDataFrame:
    """Walk an OGC API — Features collection through its next links, resumably."""
    session = requests.Session()
    session.headers["User-Agent"] = "spatial-etl/1.0 (data@example.org)"

    next_url = items_url
    if state_path and state_path.exists():
        next_url = json.loads(state_path.read_text())["next_url"]

    frames: list[gpd.GeoDataFrame] = []
    while next_url:
        response = session.get(next_url, params={"limit": page_size} if next_url == items_url else None,
                               timeout=120)
        response.raise_for_status()
        payload = response.json()

        page = gpd.GeoDataFrame.from_features(payload["features"], crs="OGC:CRS84")
        frames.append(page)

        links = {link["rel"]: link["href"] for link in payload.get("links", [])}
        next_url = links.get("next")
        if state_path:
            state_path.write_text(json.dumps({"next_url": next_url, "pages": len(frames)}))

    if not frames:
        return gpd.GeoDataFrame(geometry=[], crs="OGC:CRS84")
    return gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs="OGC:CRS84")

Persisting the next link after each page is what makes the loop resumable: a worker that dies on page 240 restarts on page 240 rather than page 1. The equivalent for WFS is the startIndex offset, which is resumable only when the service applies a stable sortBy — see paging through OGC API — Features with Python for the failure modes when it does not.

Step 5 — Validate What Arrived Before It Leaves the Stage

def assert_ingest_contract(gdf: gpd.GeoDataFrame, expected_count: int,
                           layer: str, tolerance: float = 0.001) -> None:
    """Fail loudly rather than pass a partial or transposed layer downstream."""
    if gdf.crs is None:
        raise ValueError(f"{layer}: no CRS on the ingested frame")

    if expected_count > 0:
        shortfall = abs(len(gdf) - expected_count) / expected_count
        if shortfall > tolerance:
            raise ValueError(f"{layer}: got {len(gdf)} features, service reported {expected_count}")

    minx, miny, maxx, maxy = gdf.total_bounds
    if not (-180 <= minx <= 180 and -90 <= miny <= 90):
        raise ValueError(f"{layer}: bounds {gdf.total_bounds} outside WGS84 — axis order suspect")

    invalid = (~gdf.geometry.is_valid).sum()
    if invalid:
        logging.getLogger(__name__).warning("%s: %d invalid geometries on arrival", layer, invalid)

The bounds check is the axis-order guard: latitude values above 90 cannot appear in a longitude column, so a transposed layer fails here rather than three stages later. Geometry validity is logged rather than raised, because repair belongs in geometry repair with Shapely and GeoPandas rather than in the ingestion boundary.

Advanced Patterns and Edge Cases

Filtering Server-Side With CQL2 Rather Than in Pandas

An OGC API — Features service advertising the filter conformance class accepts CQL2 expressions, which moves the predicate to the database that holds the data. Filtering land_use = 'residential' server-side over a national layer transfers thousands of features instead of millions, and the difference is entirely in transfer time rather than in cleverness. Check /conformance first, because a filter parameter an endpoint does not implement is usually ignored rather than rejected — and an ignored filter returns everything, which looks like success.

Layers That Change Underneath a Long Walk

A paged read of half a million features can take an hour, and the source may be edited during it. Cursor-based paging protects the walk’s consistency but not its completeness: features added behind the cursor are missed until the next run. Where the service exposes a timestamp property, record the run’s start time and re-query for anything modified after it on the following run, using the overlap-window logic from incremental spatial loading with watermark timestamps.

Services That Cap Silently

Many deployments enforce a maximum feature count per request that the capabilities document does not mention. The symptom is a page that returns exactly 1 000 or 5 000 features regardless of the count you asked for. Detect it by comparing the returned page length against the requested one on the first page and logging the effective cap — then size the paging loop to the real limit rather than the hoped-for one.

Performance Notes

Three measurements dominate WFS ingestion cost, and only one of them is in your control.

Transfer is usually the largest, and GeoJSON over gzip is materially smaller than GML for the same features — request compression explicitly with an Accept-Encoding header, since some servers only compress when asked. Parsing comes second: pyogrio reads GeoJSON several times faster than a pure-Python loop over json.loads, and the gap widens with attribute count. Server-side query time is third and outside your influence, which is the argument for asking for less rather than for asking faster.

Where a layer is read repeatedly, cache the raw response bodies keyed by request URL and validator, exactly as in detecting dataset changes with ETag and Last-Modified. Most public services set validators correctly, and a cached ingest costs one conditional request per page.

Integration Into ETL Pipelines

The ingestion function returns a GeoDataFrame and nothing else — no writes, no side effects — which keeps it testable against recorded fixtures. The orchestration layer decides the schedule, the retry policy and the concurrency budget, along the boundary described in orchestrating spatial ETL pipelines. Failed pages go to the dead-letter path with their request URL attached, so a replay re-issues exactly the request that failed rather than restarting the layer.

Failure-Mode Reference

Handling Authentication Where It Exists

Most public WFS endpoints are open, but the interesting ones increasingly are not. Three schemes cover nearly all of them, and each interacts differently with paging.

HTTP basic or bearer tokens are the simplest: set the header once on a requests.Session and every page inherits it. The only trap is a token that expires mid-walk, which for a long paged read is a real risk — refresh proactively rather than on the first 401, as set out in automating government portal downloads.

API keys in the query string survive paging badly, because the next link a service returns may or may not carry the key forward. Check the first next URL for the key and re-append it when the service has stripped it, rather than discovering the problem as a 401 on page two.

Session cookies issued by a login endpoint work but tie the walk to one client. Persist the cookie jar alongside the resume token if the ingestion must survive a process restart, and treat an unexpected redirect to a login page as an error rather than following it — an HTML login form parsed as GeoJSON produces a confusing exception several layers away from its cause.

Choosing Between Bulk Download and Service Reads

A WFS endpoint is not always the right way to obtain a layer, even when it is the documented one. Where the publisher also offers a full extract — a nightly GeoPackage, a regional shapefile archive — that path is usually faster, cheaper for both parties, and immune to paging bugs entirely.

The service read wins in three situations: when you need a small spatial or attribute subset of a large layer, when you need data fresher than the extract cadence, and when the extract is not published for the layer you want. Otherwise the extract wins, and the ingestion becomes a conditional download plus a local read — the pattern in detecting dataset changes with ETag and Last-Modified.

It is worth writing the comparison down per source, because the answer changes as a project grows: a pipeline that started needing one city’s parcels and now needs the country’s has usually crossed the line without anyone revisiting the decision.

Recording the Service Contract Alongside the Data

Everything discovered during negotiation — the effective page cap, the coordinate order actually returned, whether the filter class is implemented, the observed numberMatched — is knowledge that cost real requests to acquire and is lost at the end of the run unless it is written down.

Store it as a small JSON document next to the layer’s configuration and refresh it on a slow schedule rather than on every run. Two things then become possible. A change in the service — a new cap, a dropped conformance class, a different default CRS — shows up as a diff on that document rather than as a mysterious failure. And a new pipeline against the same endpoint starts from the observed contract instead of rediscovering it from scratch.

Related

• Parsing GeoJSON & Shapefile APIs
• Paging Through OGC API — Features with Python
• Handling WFS GetFeature Requests with OWSLib
• Converting GML Responses to GeoDataFrames
• Automating Government Portal Downloads