How to Download Data from a WFS Service in Python
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": "my-project/1.0 ([email protected])"}
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.
Step-by-step solution
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, extra=None):
params = {
"service": "WFS", "version": "2.0.0", "request": "GetFeature",
"typeNames": layer, "outputFormat": "application/json",
"count": size, "startIndex": start,
"sortBy": "id", # 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.
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)
print(gdf.crs, gdf.total_bounds.round(0))EPSG:28992 [116666. 379653. 170907. 510724.]
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.
Code examples
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-example/1.0 ([email protected])"}
def wfs_download(base, layer, *, page_size=1000, sort_by="id", bbox=None,
crs_urn=None, cql=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 cql:
common["CQL_FILTER"] = cql
hits = requests.get(base, params={**common, "request": "GetFeature",
"resultType": "hits"},
headers=HEADERS, timeout=120)
hits.raise_for_status()
match = re.search(r'numberMatched="(\d+)"', hits.text)
if match is None:
raise RuntimeError(f"no numberMatched in hits response: {hits.text[: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)
r.raise_for_status()
if r.text.lstrip().startswith("<"): # an exception report, with HTTP 200
raise RuntimeError(f"WFS exception: {' '.join(r.text.split())[:200]}")
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
buurten = wfs_download(
"https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0",
"wijkenbuurten:buurten",
page_size=2000,
)
print(buurten.shape, buurten.crs)wijkenbuurten:buurten: server reports 14,515 features
2,000 / 14,515
4,000 / 14,515
6,000 / 14,515
8,000 / 14,515
10,000 / 14,515
12,000 / 14,515
14,000 / 14,515
14,515 / 14,515
(14515, 32) EPSG:28992
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",
page_size=1000,
)
print(f"{len(utrecht)} neighbourhoods in the box")
print(utrecht[["buurtnaam", "gemeentenaam", "aantalInwoners"]].head(3).to_string(index=False))wijkenbuurten:buurten: server reports 421 features
421 / 421
421 neighbourhoods in the box
buurtnaam gemeentenaam aantalInwoners
Binnenstad-Noordoost Utrecht 1595
Binnenstad-Zuidwest Utrecht 2170
Wittevrouwen Utrecht 3105
421 instead of 14,515 — one request instead of eight, and the hit count reflects the filter, so the assertion still holds. The CRS URN in the bbox is load-bearing: pass (132000, 452000, …) without it against a server whose default is EPSG:4326 and you will get zero features, no error, and a very confusing hour.
Note that a bbox returns features intersecting the box, so edge neighbourhoods extend beyond it. Clip afterwards if you need them cut.
Example 3 — an attribute filter with CQL
one_city = wfs_download(
"https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0",
"wijkenbuurten:buurten",
cql="gemeentenaam='Utrecht'",
page_size=1000,
)
print(f"{len(one_city)} neighbourhoods, "
f"{one_city['aantalInwoners'].clip(lower=0).sum():,} residents")wijkenbuurten:buurten: server reports 133 features
133 / 133
133 neighbourhoods, 361,742 residents
.clip(lower=0) is there because this dataset — like many official statistical products — uses negative sentinels such as -99999 for suppressed values. Summing without handling them produces a large negative population, which is 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.
Explanation
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 exception reports arrive with HTTP 200
Many WFS servers return errors as an <ows:ExceptionReport> XML document with a 200 status, on the reasoning that the HTTP request succeeded — it is the OGC operation that failed. raise_for_status() therefore passes, and gpd.read_file() on the body produces a confusing GDAL error rather than the server's actual message.
Checking whether the body starts with < when you asked for JSON catches this in one line, and surfaces the real message:
WFS exception: <?xml version="1.0"?><ows:ExceptionReport … Could not locate
{http://wijkenbuurten.geonovum.nl}buurtn. Check the capabilities document …
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.
Edge cases or notes
- WFS 1.1 uses
typeNameandmaxFeatures, nottypeNamesandcount. Check the version in the capabilities document before assuming 2.0.0 parameter names. startIndexis 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. Always pass the full URN.
- 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_FILTERis a GeoServer extension, not part of the WFS standard. On other implementations you need an XML<Filter>in a POST body.- Some servers ignore
countabove their cap silently. Setpage_sizeto something modest — 1,000 to 5,000 — rather than trying to defeat the limit.
Internal links
- Spatial web services explained: WFS, WMS and OGC API Features — what the standards are and which to pick
- GeoJSON downloaded from an API is empty or truncated — diagnosing the failure this article prevents
- How to download administrative boundaries in Python — the layer most WFS work is after
- How to cache downloaded GIS data so you fetch it once — not re-downloading a nation every run
- GIS data sources explained — licences and provenance for service data
- How to clip spatial data in Python with GeoPandas — cutting bbox results to a real boundary
- How to retry a flaky GIS step in Python — backoff for a service that occasionally 503s
- How to calculate zonal statistics in Python — where negative sentinels cause the same damage
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 CQL_FILTER reduces one national download to a handful of pages, and the hit count still verifies it.