Spatial Web Services Explained: WFS, WMS and OGC API Features

Problem statement

A government data portal offers you the same dataset five ways: a shapefile ZIP, a WFS endpoint, a WMS endpoint, an "OGC API" link and a tile URL. Nothing on the page says which one you want, and picking wrong is expensive in a specific way:

• Choose WMS and you get a picture. There is no way to count anything in it, ever.
• Choose WFS and you get features — but the server may hand you 14,515 of them, or the first 1,000 with no indication that there were more.
• Choose the bulk file and you get everything, correct today and stale in a month.

The words are not helping. "Service", "API" and "endpoint" all describe all five. What actually distinguishes them is one question: do you receive geometry, or pixels?

Quick answer

In GeoPandas, the two that matter read almost identically:

import geopandas as gpd

# WFS — XML-era protocol, parameters in the query string
wfs = ("https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0"
       "?service=WFS&version=2.0.0&request=GetFeature"
       "&typeNames=wijkenbuurten:gemeenten&outputFormat=application/json")
gdf = gpd.read_file(wfs)
print(len(gdf), gdf.crs)

424 EPSG:28992

Ask what the service has, count, then page

1. Ask the service what it has

Every OGC service self-describes. For WFS the request is GetCapabilities, and it returns an XML document listing every layer:

import re
import requests

base = "https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0"
caps = requests.get(base, params={
    "service": "WFS", "request": "GetCapabilities", "version": "2.0.0",
}, timeout=60)

names = re.findall(r"<(?:wfs:)?Name>([^<]+)</(?:wfs:)?Name>", caps.text)
print(names)

['wijkenbuurten:buurten', 'wijkenbuurten:wijken', 'wijkenbuurten:gemeenten']

Those strings are what you pass as typeNames. They are namespaced — wijkenbuurten: is the workspace — and the prefix is part of the name, not decoration.

For OGC API Features the equivalent is a JSON document at /collections:

r = requests.get("https://demo.pygeoapi.io/master/collections", params={"f": "json"}, timeout=60)
print([c["id"] for c in r.json()["collections"]][:5])

['obs', 'lakes', 'dutch_windmills', 'dutch_castles', 'dutch_georef_stations']

Same information, one is XML and one is JSON. That difference is most of what separates the two standards in practice.

2. Find out how many features exist before you ask for them

This is the step people skip, and it is where truncation comes from. WFS has a dedicated cheap request:

r = requests.get(base, params={
    "service": "WFS", "version": "2.0.0", "request": "GetFeature",
    "typeNames": "wijkenbuurten:buurten", "resultType": "hits",
}, timeout=60)
print(re.search(r'numberMatched="(\d+)"', r.text).group(1))

14515

Fourteen thousand features. Whatever comes back from your actual request, compare it against this number.

Note the trap: resultType=hits returns XML even when you asked for outputFormat=application/json. The count lives in an attribute of the root element, not in a JSON body.

3. Page through the result

WFS 2.0 uses count and startIndex:

def page(type_name, start, size=3):
    r = requests.get(base, params={
        "service": "WFS", "version": "2.0.0", "request": "GetFeature",
        "typeNames": type_name, "outputFormat": "application/json",
        "count": size, "startIndex": start,
    }, timeout=60)
    return [f["properties"]["buurtcode"] for f in r.json()["features"]]


for start in (0, 3, 6):
    print(start, page("wijkenbuurten:buurten", start))

0 ['BU09989999', 'BU00349997', 'BU00509997']
3 ['BU00609998', 'BU00729998', 'BU00889998']
6 ['BU00939998', 'BU00969998', 'BU01669997']

No overlap between pages, which is what you need — but only because this server applies a stable ordering. A server without one can return the same feature on two pages and skip another entirely. Add sortBy on a unique attribute when the total matters.

OGC API Features uses limit and offset, plus a next link you can follow instead of computing offsets:

r = requests.get("https://demo.pygeoapi.io/master/collections/lakes/items",
                 params={"f": "json"}, timeout=60).json()
print(r["numberMatched"], r["numberReturned"], [l["rel"] for l in r["links"]])

25 10 ['self', 'alternate', 'alternate', 'alternate', 'next', 'collection']

numberMatched and numberReturned in the same response is the single best feature of the newer standard. WFS makes you ask twice.

4. Filter on the server, not in pandas

Downloading 14,515 features to keep 40 is slow for you and rude to the service. Both standards filter server-side.

Bounding box — supported everywhere:

params = {..., "bbox": "120000,480000,125000,485000,urn:ogc:def:crs:EPSG::28992"}

Always send the trailing CRS URI. Without it the server assumes its own default CRS — the Dutch grid, on this server, so metres happen to work — and a bbox in any other CRS silently returns nothing. With the EPSG:4326 URN, degrees go latitude first.

Attribute filter — the standard WFS 2.0 way is an OGC Filter Encoding expression in FILTER:

params = {..., "FILTER": '<fes:Filter xmlns:fes="http://www.opengis.net/fes/2.0"><fes:PropertyIsEqualTo>'
                         '<fes:ValueReference>gemeentenaam</fes:ValueReference><fes:Literal>Utrecht</fes:Literal>'
                         '</fes:PropertyIsEqualTo></fes:Filter>'}

The shorter CQL_FILTER that many tutorials show is a GeoServer extension. PDOK runs MapServer, which ignores it and returns the whole layer with HTTP 200 — the WFS guide measures exactly that.

OGC API Features puts attributes straight in the query string:

params = {"f": "json", "bbox": "4.3,51.9,4.5,52.1", "datetime": "2026-01-01/.."}

5. Check the CRS you were handed

print(gdf.crs, gdf.total_bounds.round(1))          # the municipalities from the quick answer

EPSG:28992 [ 10425.2 306846.2 278026.1 621876.3]

Those are metres in the Dutch national grid, not degrees. GeoJSON from a WFS often carries a non-WGS84 CRS declaration, which technically contradicts the GeoJSON specification but is extremely common and is what you want. GeoPandas reads it correctly; code that assumes "GeoJSON means lat/lon" does not.

A service inventory before committing

Example 1 — a service inventory before you commit to one

import re
import requests

HEADERS = {"User-Agent": "spatialworkflow-example/1.0"}


def wfs_inventory(base):
    caps = requests.get(base, params={
        "service": "WFS", "request": "GetCapabilities", "version": "2.0.0",
    }, headers=HEADERS, timeout=60)
    caps.raise_for_status()

    layers = re.findall(r"<(?:wfs:)?Name>([^<]+)</(?:wfs:)?Name>", caps.text)
    formats = sorted(set(re.findall(r"<(?:ows:)?Value>(application/[^<]+)</", caps.text)))

    rows = []
    for layer in layers:
        hits = requests.get(base, params={
            "service": "WFS", "version": "2.0.0", "request": "GetFeature",
            "typeNames": layer, "resultType": "hits",
        }, headers=HEADERS, timeout=120)
        match = re.search(r'numberMatched="(\d+)"', hits.text)
        rows.append((layer, int(match.group(1)) if match else -1))

    print(f"formats: {formats}")
    for layer, n in rows:
        print(f"  {layer:34} {n:>8,} features")
    return rows


wfs_inventory("https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0")

formats: ['application/gml+xml; version=3.2', 'application/json', 'application/json; subtype=geojson']
  wijkenbuurten:buurten                14,515 features
  wijkenbuurten:wijken                  3,434 features
  wijkenbuurten:gemeenten                 424 features

Three numbers that decide the whole approach. 424 municipalities is one request. 14,515 neighbourhoods needs paging or a bbox. And application/json being in the format list is what makes gpd.read_file() on the URL work at all — plenty of servers offer only GML.

Example 2 — the same layer from a WFS and an OGC API service

import geopandas as gpd
import requests

def from_wfs(base, layer, bbox=None, crs_urn="urn:ogc:def:crs:EPSG::28992"):
    params = {
        "service": "WFS", "version": "2.0.0", "request": "GetFeature",
        "typeNames": layer, "outputFormat": "application/json",
    }
    if bbox:
        params["bbox"] = ",".join(str(v) for v in bbox) + "," + crs_urn
    r = requests.get(base, params=params, headers=HEADERS, timeout=180)
    r.raise_for_status()
    return gpd.read_file(r.content)


def from_ogcapi(base, collection, bbox=None, limit=1000):
    params = {"f": "json", "limit": limit}
    if bbox:
        params["bbox"] = ",".join(str(v) for v in bbox)     # always WGS84 here
    r = requests.get(f"{base}/collections/{collection}/items",
                     params=params, headers=HEADERS, timeout=180)
    r.raise_for_status()
    payload = r.json()
    print(f"  returned {payload.get('numberReturned')} of {payload.get('numberMatched')}")
    return gpd.GeoDataFrame.from_features(payload["features"], crs="EPSG:4326")


gemeenten = from_wfs("https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0",
                     "wijkenbuurten:gemeenten")
lakes = from_ogcapi("https://demo.pygeoapi.io/master", "lakes")
print(len(gemeenten), gemeenten.crs, "|", len(lakes), lakes.crs)

  returned 10 of 25
424 EPSG:28992 | 10 EPSG:4326

Two things the output tells you immediately. The WFS layer arrived in the national grid, in metres — usable for area and distance without reprojection. The OGC API request asked for 1,000 and got 10, because the server caps limit and says so in numberReturned. That is the check WFS makes you do by hand.

Example 3 — deciding between a service and a bulk download

def access_advice(n_features, area_fraction, refetch_per_month):
    """Rough guidance: services suit small slices of changing data."""
    if n_features > 500_000:
        return "bulk file — too many features for any service"
    if area_fraction < 0.1 and refetch_per_month <= 4:
        return "service with a bbox — small slice, low volume"
    if refetch_per_month > 20:
        return "bulk file + local cache — you are hammering the service"
    return "service, paged, cached locally"


cases = [
    ("all Dutch neighbourhoods, monthly", 14_515, 1.0, 1),
    ("one city's neighbourhoods, daily", 14_515, 0.02, 30),
    ("national buildings, once", 10_000_000, 1.0, 1),
]
for label, n, frac, freq in cases:
    print(f"{label:38} -> {access_advice(n, frac, freq)}")

all Dutch neighbourhoods, monthly      -> service, paged, cached locally
one city's neighbourhoods, daily       -> bulk file + local cache — you are hammering the service
national buildings, once               -> bulk file — too many features for any service

The middle case is the counter-intuitive one. A small area queried often is worse for a service than a large area queried rarely, because every request costs the provider the same overhead. Fetch once, cache locally, and refresh on a schedule.

WMS cannot answer questions; WFS truncates silently

Why WMS can never answer a question

A WMS GetMap request runs the query, styles the result and rasterises it on the server. What crosses the network is an image. The features that produced it — their geometry, their attributes, their count — do not exist in the response.

GetFeatureInfo is the partial exception: click a pixel, get the attributes of whatever rendered there. It is one feature at a time, at the server's discretion, and it is not a substitute for WFS.

The practical rule from choosing a data source: if your deliverable contains a number, you needed WFS. If it contains only a backdrop, WMS is cheaper and faster.

Why WFS truncation is silent by design

WFS 2.0 servers set a maximum features limit — PDOK's is 1,000 — and applying it is a normal, spec-compliant response. The document you get back is valid and complete-looking. numberMatched and numberReturned exist as attributes on the root element of the GML response, but the GeoJSON output of many servers drops them.

So the defence is procedural, not technical: always issue a resultType=hits request first and compare. That is the entire content of GeoJSON downloaded from an API is empty or truncated, and it applies to every service in this article.

Why OGC API Features is worth preferring

It is not that REST is better than XML. It is three concrete things:

• numberMatched and numberReturned in every response, so truncation is visible without a second request.
• A next link, so paging does not require you to compute offsets correctly.
• JSON throughout, so the capabilities document is json.loads-able rather than regex-able.

WFS is not going away — the installed base is enormous — but where a provider offers both, the newer one has fewer ways to be silently wrong.

Why the CRS parameter is fiddly

WFS 1.0 interpreted EPSG:4326 as longitude-latitude. WFS 1.1 and 2.0 follow the authority's own axis order, which for EPSG:4326 is latitude first. This is the single most common source of "my features are in the Indian Ocean" from a WFS.

The fix is to be explicit everywhere: pass the full CRS URN in the bbox parameter, request srsName explicitly, and check total_bounds against something you know before trusting the result. See points plotting in the ocean for the general form of this failure.

Output formats and WFS version differences

• outputFormat=application/json is not guaranteed. Many servers offer only GML. gpd.read_file() can read GML through GDAL, but attribute types and namespaces come through less cleanly.
• typeNames is typeName (singular) in WFS 1.1. Sending the wrong one produces an exception report, not an error status.
• Paging without sortBy is not stable on all servers. If the total matters, sort on a unique attribute.
• A bbox without a CRS URN uses the layer's declared CRS. Coordinates in the wrong units usually return zero features rather than an error.
• Exception reports come back with HTTP 200 on many WFS servers. The body is an XML <ows:ExceptionReport>. Check for it before parsing.
• WCS is the raster sibling and returns actual cell values, unlike WMS. If you need elevation or reflectance numbers rather than a picture, WCS or a STAC catalogue is what you want.
• Rate limits are rarely documented. Treat any public service as fragile: cache, back off, and do not run an unthrottled loop against one.

Internal links

• How to download data from a WFS service in Python — the practical, paged implementation
• GIS data sources explained — where services sit among bulk files and APIs
• GeoJSON downloaded from an API is empty or truncated — the truncation failure in detail
• STAC explained: how satellite imagery catalogues work — the equivalent standard for imagery
• How to cache downloaded GIS data so you fetch it once — being a good citizen of a shared service
• How to add a basemap to a GeoPandas plot — the legitimate use of rendered tiles
• Points plot in the ocean off Africa — the axis-order failure this article warns about
• WKT, WKB and GeoJSON: how geometry is actually stored — the formats these services speak

FAQ

What is the difference between WFS and WMS?

WFS returns features — geometry and attributes you can analyse. WMS returns a rendered image. If you need to count, measure or join anything, you need WFS.

Why did my WFS request return exactly 1,000 features?

That is the server's maximum features limit, applied silently. Issue a resultType=hits request to learn the real total, then page with count and startIndex.

Is OGC API Features a replacement for WFS?

It is the successor standard. Where a provider offers both, prefer it: every response reports how many features matched versus how many were returned, which removes the commonest source of silent error.

Why are my WFS coordinates swapped?

WFS 2.0 follows EPSG axis order, which is latitude-first for EPSG:4326. Request a projected CRS explicitly with srsName, or reproject and verify against known bounds.

Can I use gpd.read_file() directly on a WFS URL?

Yes, if the server supports outputFormat=application/json and the result fits in one response. For anything paged, build the requests yourself so you can check the counts.

How do I filter on an attribute?

The standard OGC Filter Encoding FILTER parameter works on any WFS 2.0 server. CQL_FILTER works only where the server supports it (GeoServer) — elsewhere, as on PDOK's MapServer, it is silently ignored and you get the whole layer. OGC API Features takes attribute parameters directly in the query string.

What is WCS for?

Raster values rather than rendered pixels — elevation, temperature, reflectance. It is the WFS of the raster world, and much less widely deployed than WMS.