⚙️ How-To Data sources & downloads Files & formats Finding & Fetching Data · step 7 Notebook By · Updated

How to Download Data from a WFS Service in Python

On this page

Problem statement

You have a WFS URL and you want the data as a GeoDataFrame. The one-liner works, right up until it does not:

import geopandas as gpd

url = ("https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0"
       "?service=WFS&version=2.0.0&request=GetFeature"
       "&typeNames=wijkenbuurten:buurten&outputFormat=application/json")
gdf = gpd.read_file(url)
print(len(gdf))
1000

There are 14,515 neighbourhoods in that layer. You received 1,000, the server considers that a complete and correct response, and nothing in the GeoDataFrame indicates otherwise.

The fix is not a bigger count parameter — most servers cap it and ignore anything larger. The fix is to find out the real total first, then page until you have all of it.

Quick answer

Ask for the count, then loop:

import re
import geopandas as gpd
import pandas as pd
import requests

BASE = "https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0"
HEADERS = {"User-Agent": "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"}


def wfs_count(base, layer):
    r = requests.get(base, params={
        "service": "WFS", "version": "2.0.0", "request": "GetFeature",
        "typeNames": layer, "resultType": "hits",
    }, headers=HEADERS, timeout=120)
    r.raise_for_status()
    return int(re.search(r'numberMatched="(\d+)"', r.text).group(1))


print(wfs_count(BASE, "wijkenbuurten:buurten"))
14515

That number is the contract. Whatever you download must match it, or you know the download is incomplete.

A paged WFS download that verifies itself against the server hit count. the final assertion is the point of the loop1GetCapabilities — layer names and formatsonce2resultType=hits — the real totalXML, not JSON3GetFeature with count, startIndex, sortByrepeat4concat, then assert against the totalraise on mismatchsortBy is what makes startIndex mean anything spatialworkflow.io
The count request costs one round trip and turns silent truncation into an assertion.

Layer names, JSON support, hit counts, paging

1. Discover the layer names

typeNames must match exactly, prefix included. Get them from the capabilities document rather than the portal page:

caps = requests.get(BASE, params={
    "service": "WFS", "request": "GetCapabilities", "version": "2.0.0",
}, headers=HEADERS, timeout=60)
print(re.findall(r"<(?:wfs:)?Name>([^<]+)</(?:wfs:)?Name>", caps.text))
['wijkenbuurten:buurten', 'wijkenbuurten:wijken', 'wijkenbuurten:gemeenten']

2. Confirm the server can give you JSON

print(sorted(set(re.findall(r"<(?:ows:)?Value>(application/[^<]+)</", caps.text))))
['application/gml+xml; version=3.2', 'application/json', 'application/vnd.google-earth.kml+xml']

If application/json is absent, drop outputFormat entirely and let the server return GML — GDAL reads it, though attribute types arrive less cleanly and namespaced column names need tidying.

3. Get the hit count

resultType=hits returns metadata with no features. Two things to know: it is cheap, and it comes back as XML even when you asked for JSON, so parse the attribute rather than calling .json().

4. Page with count and startIndex

def wfs_page(base, layer, start, size, sort_by, extra=None):
    params = {
        "service": "WFS", "version": "2.0.0", "request": "GetFeature",
        "typeNames": layer, "outputFormat": "application/json",
        "count": size, "startIndex": start,
        "sortBy": sort_by,                 # stable ordering — see below
    }
    params.update(extra or {})
    r = requests.get(base, params=params, headers=HEADERS, timeout=180)
    r.raise_for_status()
    return gpd.read_file(r.content)

sortBy is not optional if the total matters. Without a defined order, a server is free to return rows in whatever order the query planner produced, and two pages can overlap while a third feature is never returned at all. Sort on something unique — and on an attribute: the GeoJSON from this server carries an id column, but that is the feature identifier, and sortBy=id is refused with HTTP 400. The layer's own key is buurtcode (or gemeentecode for municipalities).

5. Stop on the right condition

Two conditions end the loop, and you need both:

  • a page comes back with fewer features than you asked for — the last page
  • the accumulated total reaches numberMatched — the assertion

If the first fires before the second, the download is short and you should raise rather than continue.

6. Check the CRS you were given

gdf = wfs_page(BASE, "wijkenbuurten:gemeenten", 0, 5, sort_by="gemeentecode")
print(gdf.crs, gdf.total_bounds.round(0))
EPSG:28992 [132850. 479382. 269942. 592678.]

Metres in the Dutch national grid. A GeoJSON response carrying a non-WGS84 CRS is common from WFS servers and is what you want here — but any code that assumes "GeoJSON is lat/lon" will now be wrong by 100 km.

Requests for progressively more features all returning the server cap. asking for more returns the cap, not an errorcount=100100 returnedcount=10001,000 returnedcount=1000001,000 returnedactually exist14,515014.5kEvery row above is HTTP 200 with a valid document spatialworkflow.io
Asking for more than the cap does not raise. It returns the cap.

A paged download that verifies itself

Example 1 — a complete paged download that verifies itself

import re
import time

import geopandas as gpd
import pandas as pd
import requests

HEADERS = {"User-Agent": "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"}


def wfs_download(base, layer, *, sort_by, page_size=1000, bbox=None,
                 crs_urn=None, fes_filter=None, pause=0.5):
    """Download every feature of a WFS layer, verified against the server's own count."""
    common = {"service": "WFS", "version": "2.0.0", "typeNames": layer}
    if bbox:
        if crs_urn is None:
            raise ValueError("a bbox needs its CRS URN, or the server guesses the units")
        common["bbox"] = ",".join(str(v) for v in bbox) + "," + crs_urn
    if fes_filter:
        common["FILTER"] = fes_filter           # OGC Filter Encoding 2.0 — every WFS 2.0 server

    hits = requests.get(base, params={**common, "request": "GetFeature",
                                      "resultType": "hits"},

                        headers=HEADERS, timeout=120)
    match = re.search(r'numberMatched="(\d+)"', hits.text)
    if match is None:                                # an exception report, whatever the status
        message = " ".join(re.sub("<[^>]+>", " ", hits.text).split())
        raise RuntimeError(f"WFS exception (HTTP {hits.status_code}): {message[:200]}")
    expected = int(match.group(1))
    print(f"{layer}: server reports {expected:,} features")

    frames, start = [], 0
    while start < expected:
        r = requests.get(base, params={
            **common, "request": "GetFeature", "outputFormat": "application/json",
            "count": page_size, "startIndex": start, "sortBy": sort_by,
        }, headers=HEADERS, timeout=300)
        if r.text.lstrip().startswith("<"):          # an exception report, whatever the status
            message = " ".join(re.sub("<[^>]+>", " ", r.text).split())
            raise RuntimeError(f"WFS exception (HTTP {r.status_code}): {message[:200]}")
        r.raise_for_status()

        page = gpd.read_file(r.content)
        frames.append(page)
        start += len(page)
        print(f"  {start:>7,} / {expected:,}")
        if len(page) < page_size and start < expected:
            raise RuntimeError(f"short page at {start}: got {len(page)}, asked {page_size}")
        time.sleep(pause)

    gdf = gpd.GeoDataFrame(pd.concat(frames, ignore_index=True), crs=frames[0].crs)
    if len(gdf) != expected:
        raise RuntimeError(f"downloaded {len(gdf):,}, server said {expected:,}")
    return gdf



try:                                                # a page size above the server's cap
    wfs_download(BASE, "wijkenbuurten:buurten", sort_by="buurtcode", page_size=2000)
except RuntimeError as err:
    print(f"RuntimeError: {err}")

buurten = wfs_download(BASE, "wijkenbuurten:buurten", sort_by="buurtcode", page_size=1000)
print(buurten.shape, buurten.crs, f"buurtcode unique: {buurten['buurtcode'].is_unique}")
print(f"aantalInwoners: min {buurten['aantalInwoners'].min():,}, "
      f"{(buurten['aantalInwoners'] < 0).sum()} negative values")
wijkenbuurten:buurten: server reports 14,515 features
    1,000 / 14,515
RuntimeError: short page at 1000: got 1000, asked 2000
wijkenbuurten:buurten: server reports 14,515 features
    1,000 / 14,515
    2,000 / 14,515
    3,000 / 14,515
    4,000 / 14,515
    5,000 / 14,515
    6,000 / 14,515
    7,000 / 14,515
    8,000 / 14,515
    9,000 / 14,515
   10,000 / 14,515
   11,000 / 14,515
   12,000 / 14,515
   13,000 / 14,515
   14,000 / 14,515
   14,515 / 14,515
(14515, 240) EPSG:28992 buurtcode unique: True
aantalInwoners: min -99,997, 94 negative values

The final assertion is the point of the whole function. Every other line exists to make that comparison possible.

Example 2 — fetching only your study area

from shapely.geometry import box

utrecht = wfs_download(
    "https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0",
    "wijkenbuurten:buurten",
    bbox=(132000, 452000, 140000, 460000),
    crs_urn="urn:ogc:def:crs:EPSG::28992",
    sort_by="buurtcode",
)
print(f"{len(utrecht)} neighbourhoods in the box")
print(utrecht[["buurtnaam", "gemeentenaam", "aantalInwoners"]].head(3).to_string(index=False))
wijkenbuurten:buurten: server reports 114 features
      114 / 114
114 neighbourhoods in the box
             buurtnaam gemeentenaam  aantalInwoners
        Groenekan Kern      De Bilt            1305
Groenekan Buitengebied      De Bilt             660
           Weltevreden      De Bilt              30

A hundred-odd instead of 14,515 — one request instead of fifteen, and the hit count reflects the filter, so the assertion still holds. The CRS URN in the bbox is load-bearing: this server happens to default to the Dutch grid, but pass (132000, 452000, …) without it to a server whose default is EPSG:4326 and you get zero features, no error, and a very confusing hour — and if you pass degrees with the EPSG:4326 URN, they must be latitude first (see GeoJSON from an API is empty or truncated).

Note that a bbox returns features intersecting the box, so edge neighbourhoods extend beyond it. Clip afterwards if you need them cut.

Hit counts PDOK returned for the neighbourhood layer under each way of restricting the request. what the server counted for each requestrequesthitswhat happenedno filter14,515the whole layerCQL_FILTER Utrecht14,515ignored, HTTP 200FES FILTER Utrecht111appliedbbox + EPSG:28992 URN114appliedA filter the server does not know is a filter it does not apply — and it says nothing spatialworkflow.io
Ask for the hit count with the filter attached, and a filter that silently did nothing shows up as the full total.

Example 3 — an attribute filter that the server actually applies

The filter most tutorials show is CQL_FILTER. Check what this server makes of it before trusting it:

def hits(base, layer, **extra):
    r = requests.get(base, params={"service": "WFS", "version": "2.0.0", "request": "GetFeature",
                                   "typeNames": layer, "resultType": "hits", **extra},
                     headers=HEADERS, timeout=120)
    return int(re.search(r'numberMatched="(\d+)"', r.text).group(1))


def equals(prop, value):
    """An OGC Filter Encoding 2.0 equality test, the filter every WFS 2.0 server must accept."""
    return ('<fes:Filter xmlns:fes="http://www.opengis.net/fes/2.0"><fes:PropertyIsEqualTo>'
            f"<fes:ValueReference>{prop}</fes:ValueReference><fes:Literal>{value}</fes:Literal>"
            "</fes:PropertyIsEqualTo></fes:Filter>")


print("no filter:          ", hits(BASE, "wijkenbuurten:buurten"))
print("CQL_FILTER Utrecht: ", hits(BASE, "wijkenbuurten:buurten", CQL_FILTER="gemeentenaam='Utrecht'"))
print("FES filter Utrecht: ", hits(BASE, "wijkenbuurten:buurten", FILTER=equals("gemeentenaam", "Utrecht")))
no filter:           14515
CQL_FILTER Utrecht:  14515
FES filter Utrecht:  111

CQL_FILTER is a GeoServer extension. PDOK runs MapServer, which ignores the unknown parameter and returns the whole layer with HTTP 200 — so a "filtered" download would quietly have been all fourteen thousand neighbourhoods. The FES filter is part of the WFS 2.0 standard:

one_city = wfs_download(BASE, "wijkenbuurten:buurten", sort_by="buurtcode",
                        fes_filter=equals("gemeentenaam", "Utrecht"))
residents = one_city["aantalInwoners"]
print(f"{len(one_city)} neighbourhoods; aantalInwoners min {residents.min():,}, "
      f"sum {residents.sum():,}, sum of non-negative {residents.clip(lower=0).sum():,}")
wijkenbuurten:buurten: server reports 111 features
      111 / 111
111 neighbourhoods; aantalInwoners min 0, sum 367,715, sum of non-negative 367,715

Utrecht's neighbourhoods happen to be clean, but look at the minimum before you trust any sum: across the whole national layer, example 1 found 94 negative populations, down to -99,997. Official statistical products use negative sentinels like these for suppressed or unavailable values, and summing them produces a wrong population with no error — the same class of bug as NoData in a raster sum. Check the minimum of any numeric column from an official source before aggregating it.

Why the server cap is invisible

Why the server cap is invisible

A WFS server is configured with a maximum number of features per response. Applying it is spec-compliant behaviour, and the response is a complete, valid document that happens to contain fewer features than exist. There is no status code for "partial", no header, and — in the GeoJSON output of many servers — no numberMatched field either.

The resultType=hits request exists precisely to close this gap. It is one cheap round trip that converts an invisible failure into an assertion, which is why every example here starts with it.

Why sortBy matters more than it looks

startIndex is an offset into a result set. If the server does not impose an order, "the result set" is whatever the underlying query returned this time — and on a database under concurrent write load, that ordering can genuinely differ between requests.

The symptom is horrible to debug: a download that assertion-passes on count but contains 40 duplicates and is missing 40 other features. Sorting on a unique attribute makes the offset meaningful. Verify afterwards:

print(buurten["buurtcode"].duplicated().sum(), "duplicate codes")
0 duplicate codes

Why to read the body before the status

WFS servers return errors as an XML exception report, and they disagree about the status code. Some send it with HTTP 200, on the reasoning that the HTTP request succeeded and only the OGC operation failed; then raise_for_status() passes and gpd.read_file() on the body produces a confusing GDAL error. PDOK sends HTTP 400 — and then raise_for_status() raises a bare 400 Client Error: Bad Request that throws the server's explanation away.

Checking whether the body starts with < before checking the status handles both, and surfaces the real message:

try:
    wfs_download(BASE, "wijkenbuurten:buurtn", sort_by="buurtcode")
except RuntimeError as err:
    print(err)
WFS exception (HTTP 400): msWFSGetFeature(): WFS server error. TYPENAME 'buurtn' doesn't exist in this server. Please check the capabilities and reformulate your request.
Four checks that turn a WFS download from a guess into a verified fetch. four checks, one line eachapplication/json is in the capabilities format listelse GMLresultType=hits gives the real totalbefore fetchingsortBy is set on a unique attributestable pagingdownloaded count equals numberMatchedor raiseThe body starting with "<" is an exception report, not data spatialworkflow.io
Four checks, each one line. Together they are the difference between a download and a guess.

Why to be gentle with the service

These are usually public services with no rate limit published and a real one enforced. The time.sleep(pause) between pages is not superstition — a tight loop of 15 requests will get you throttled on some servers and blocked on others.

For anything scheduled, fetch once and cache the result with an expiry, rather than re-downloading a national dataset every time your pipeline runs.

WFS 1.1 names and zero-based indexes

  • WFS 1.1 uses typeName and maxFeatures, not typeNames and count. Check the version in the capabilities document before assuming 2.0.0 parameter names.
  • startIndex is zero-based in WFS 2.0 but some implementations treat it as one-based. Compare the first feature of page 0 and page 1 once, on a small page size, before trusting a long run.
  • gpd.read_file() accepts bytes directly, which avoids writing temporary files: gpd.read_file(response.content).
  • Negative numbers in official statistics are often suppression sentinels, not values. Check .min() on every numeric column before aggregating.
  • A bbox with no CRS URN is a silent zero-result generator on any server whose default CRS is not the one you meant. Always pass the full URN — and with the EPSG:4326 URN, latitude first.
  • GML output preserves types better than the GeoJSON output on some servers, at the cost of namespaced column names you will need to rename.
  • CQL_FILTER is a GeoServer extension, not part of the WFS standard. PDOK's MapServer ignores it without a word and returns the whole layer; the standard OGC Filter Encoding FILTER parameter works everywhere.
  • Servers clamp count to their cap silently. Set page_size at or below the cap: on this page 2,000 against PDOK's 1,000 tripped the short-page guard on the first page.

FAQ

Why did I get exactly 1,000 features?

That is the server's per-response cap. Ask for the real total with resultType=hits and page through with count and startIndex.

Can I just set count to a huge number?

No. Servers cap it and return the cap without complaint. Paging is the only reliable approach.

Why does resultType=hits return XML when I asked for JSON?

Because the hit count lives in attributes of the WFS response root element, which is a GML construct. Parse numberMatched with a regex rather than calling .json().

My bbox returns nothing. What is wrong?

Almost always units. Append the CRS URN to the bbox parameter — …,urn:ogc:def:crs:EPSG::28992 — so the server knows whether your numbers are metres or degrees.

Why are there duplicate features in my paged download?

The server had no stable ordering. Add sortBy on a unique attribute and check .duplicated() on that column afterwards.

The response has HTTP 200 but GeoPandas cannot read it. Why?

It is probably an <ows:ExceptionReport> — an OGC error delivered with a 200 status. Check whether the body starts with < and surface the message.

Should I download the whole layer or filter it?

Filter, if your study area is stable. A bbox or an FES filter reduces one national download to a handful of pages, and the hit count still verifies it.