Python voor GIS met GeoPandas, Shapely, pyproj en PDOK

Vijf bibliotheken vormen de kern, en ze vullen elkaar aan. GeoPandas beschrijft zichzelf zo:

GeoPandas is a project to add support for geographic data to pandas objects.

Bron: README van GeoPandas op GitHub, gelezen op 11 oktober 2026. GeoPandas rekent met Shapely-geometrieën, leest en schrijft bestanden via pyogrio (dat bij de installatie vanzelf meekwam) en gebruikt pyproj voor coördinaatstelsels.

Let op de minimale Python-versie: pyproj, rasterio en OWSLib vragen 3.12. De wheels van rasterio en pyogrio brengen elk hun eigen GDAL mee; in de testomgeving waren dat GDAL 3.12.2 en 3.12.4, naast GDAL 3.13.3 op de commandoregel. Een oudere omgeving op een beheerserver kan de installatie blokkeren. Over beheer en hosting van zulke omgevingen staat meer op GIS-beheer en hosting.

U zoekt het adres met de Locatieserver, haalt de panden eromheen op via de BAG OGC API en analyseert ze in GeoPandas. Het script hieronder doet dat in vijf stappen en schrijft het resultaat naar een GeoPackage. Installeer eerst geopandas, pyproj, owslib en requests met pip.

"""Panden rond een adres ophalen bij PDOK en analyseren met GeoPandas."""
import requests
import geopandas as gpd
from pyproj import Transformer
from shapely.geometry import Point, box
from owslib.ogcapi.features import Features

LOCATIESERVER = "https://api.pdok.nl/bzk/locatieserver/search/v3_1/free"
BAG = "https://api.pdok.nl/kadaster/bag/ogc/v2"
RD = "EPSG:28992"

# 1. Adres naar RD-coördinaat (Locatieserver geeft RD en WGS84 terug)
r = requests.get(
    LOCATIESERVER,
    params={"q": "Stationsplein 12 Amersfoort", "fq": "type:adres", "fl": "weergavenaam,centroide_rd,centroide_ll", "rows": 1},
    timeout=30,
)
r.raise_for_status()
doc = r.json()["response"]["docs"][0]
punt_rd = Point(*map(float, doc["centroide_rd"].removeprefix("POINT(").removesuffix(")").split()))
print(doc["weergavenaam"], punt_rd.wkt)

# 2. Zelf omrekenen met pyproj en vergelijken met wat PDOK teruggeeft
lon, lat = map(float, doc["centroide_ll"].removeprefix("POINT(").removesuffix(")").split())
naar_rd = Transformer.from_crs("EPSG:4326", RD, always_xy=True)  # always_xy: lengte, breedte
x, y = naar_rd.transform(lon, lat)
print(f"pyproj: {x:.3f} {y:.3f}; verschil met PDOK: {abs(x - punt_rd.x):.3f} m, {abs(y - punt_rd.y):.3f} m")

# 3. Panden binnen een vierkant van 400 m via de BAG OGC API Features (OWSLib)
vlak = punt_rd.buffer(200, cap_style="square")
bag = Features(BAG)
data = bag.collection_items(
    "pand",
    bbox=[round(c) for c in vlak.bounds],
    **{"bbox-crs": "http://www.opengis.net/def/crs/EPSG/0/28992", "crs": "http://www.opengis.net/def/crs/EPSG/0/28992"},
    limit=500,
)
if any(link["rel"] == "next" for link in data["links"]):
    raise SystemExit("Meer panden dan limit: volg de next-link of verklein het vlak")
panden = gpd.GeoDataFrame.from_features(data["features"], crs=RD)
print(len(panden), "panden, kolommen:", list(panden.columns)[:6])

# 4. Analyse: alleen panden in gebruik, oppervlak, afstand tot het adres
gebruik = panden[panden["status"] == "Pand in gebruik"].copy()
gebruik["m2"] = gebruik.geometry.area.round(1)
gebruik["afstand_m"] = gebruik.geometry.distance(punt_rd).round(1)
print(gebruik.nsmallest(3, "afstand_m")[["identificatie", "bouwjaar", "m2", "afstand_m"]].to_string(index=False))

# 5. Opslaan als GeoPackage (RD blijft behouden)
gebruik.to_file("panden.gpkg", layer="pand", driver="GPKG")
print("opgeslagen:", gpd.read_file("panden.gpkg", layer="pand").crs.to_string())

Uitvoer op 11 oktober 2026, met Python 3.13.1 en de versies uit de tabel:

Stationsplein 12, Amersfoort POINT (154201.236 462720.566)
pyproj: 154201.236 462720.567; verschil met PDOK: 0.000 m, 0.001 m
244 panden, kolommen: ['geometry', 'aantal_verblijfsobjecten', 'bouwjaar', 'documentdatum', 'documentnummer', 'gebruiksdoel']
   identificatie  bouwjaar     m2  afstand_m
0307100000595635      2024   10.3        0.0
0307100000341963      1970  455.0       12.0
0307100000347585      1973 2436.6       23.5
opgeslagen: EPSG:28992

De Locatieserver en de BAG leveren hier hun gegevens zonder sleutel. De uitvoer verandert mee met de BAG, dus een ander bouwjaar of aantal op een latere dag is geen fout. Meer over de Locatieserver staat bij adressen zoeken met PDOK, over de BAG op BAG, en over de API-standaard bij WFS en OGC API Features.

U maakt een Transformer tussen twee stelsels en geeft always_xy=True mee. In het script hierboven liet dat de eigen berekening binnen een millimeter uitkomen op de coördinaten die PDOK zelf teruggeeft: een verschil van 0,000 meter in oostrichting en 0,001 meter in noordrichting.

from pyproj import Transformer

# always_xy=True: volgorde (lengte, breedte) in en (oost, noord) uit, voor elk stelsel
naar_rd = Transformer.from_crs("EPSG:4326", "EPSG:28992", always_xy=True)
print(naar_rd.transform(5.37553215, 52.15266021))
# (154201.23591069333, 462720.56677658844)

# zonder always_xy verwacht EPSG:4326 (breedte, lengte); omgekeerd geeft onzin, geen fout
zonder = Transformer.from_crs("EPSG:4326", "EPSG:28992")
print(zonder.transform(5.37553215, 52.15266021))
# (6349628.560937857, -3622470.4804081977)

De tweede aanroep laat zien waarom de optie bestaat. Het stelsel EPSG:4326 verwacht volgens de standaard eerst de breedtegraad en dan de lengtegraad. Wie (lengte, breedte) aanlevert zonder die optie, krijgt geen foutmelding maar het getal 6349628 en -3622470, ver buiten Nederland.

Pyproj downloadt in zijn standaardinstelling geen rasterbestanden van het internet: pyproj.network.is_network_enabled() gaf False. Voor dit RD-voorbeeld was dat niet nodig. Sommige omrekeningen kunnen rasterbestanden gebruiken, dus controleer bij nauwkeurig werk welke transformatie pyproj kiest. Voor een snelle controle zonder code is er de RD naar WGS84-tool.

U haalt het oppervlaktemodel (DSM) en het terreinmodel (DTM) van het AHN op via de WCS van PDOK, leest ze met rasterio als gemaskeerde arrays en trekt per pand het maaiveld van de dakhoogte af. Het script bouwt voort op panden.gpkg uit het vorige voorbeeld.

import geopandas as gpd, numpy as np, requests
from rasterio.features import geometry_mask
from rasterio.io import MemoryFile

WCS = "https://service.pdok.nl/rws/ahn/wcs/v1_0"
panden = gpd.read_file("panden.gpkg", layer="pand")
xmin, ymin, xmax, ymax = (round(c) for c in panden.total_bounds)

def raster(coverage):
    r = requests.get(WCS, timeout=60, params=[
        ("SERVICE", "WCS"), ("VERSION", "2.0.1"), ("REQUEST", "GetCoverage"),
        ("COVERAGEID", coverage), ("FORMAT", "image/tiff"),
        ("SUBSET", f"x({xmin},{xmax})"), ("SUBSET", f"y({ymin},{ymax})")])
    r.raise_for_status()
    with MemoryFile(r.content) as mem, mem.open() as src:
        return src.read(1, masked=True), src.transform  # NoData telt niet mee

dsm, transform = raster("dsm_05m")
dtm, _ = raster("dtm_05m")
print("DTM zonder waarde:", f"{dtm.mask.mean():.0%}")

def mediaan(arr, geom):
    w = arr[~geometry_mask([geom], arr.shape, transform)].compressed()
    return float(np.median(w)) if w.size else np.nan

# dak: DSM binnen het pand; maaiveld: DTM in een rand van 3 m eromheen
panden["hoogte_m"] = [round(mediaan(dsm, g) - mediaan(dtm, g.buffer(3).difference(g)), 1)
                      for g in panden.geometry]
print(panden.nlargest(3, "hoogte_m")[["identificatie", "bouwjaar", "hoogte_m"]].to_string(index=False))

Uitvoer op 12 oktober 2026, met rasterio 1.5.2:

DTM zonder waarde: 33%
   identificatie  bouwjaar  hoogte_m
0307100000308550      1999      38.3
0307100000316048      1992      27.8
0307100000338371      1995      27.7

De eerste regel verklaart de opzet. Het DTM heeft onder gebouwen geen waarde; een eerste versie die DSM en DTM op één punt binnen elk pand vergeleek, gaf daardoor voor 205 van de 211 panden geen hoogte. Het maaiveld uit een rand rond het pand lost dat op. De uitkomst is een schatting, geen gemeten bouwhoogte.

De eerste twee valkuilen kwamen bij het testen van dit script echt voor. De rest volgt uit de documentatie.

• Coördinaatvolgorde. Zonder always_xy=True keert de volgorde om voor stelsels met breedte eerst. Zie de meting hierboven.
• OWSLib wil een lijst voor bbox. Een tekenreeks als "154000,462520,154400,462920" gaf de foutmelding dat de waarde van bbox leeg was. OWSLib voegt de elementen van een lijst zelf samen.
• Het limiet en de paginering. Het script vraagt maximaal 500 panden. Zijn het er meer, dan staat er een next-link in het antwoord, en het script stopt dan met een melding. Zonder die controle werkt u stilzwijgend met een deel van de data.
• Het stelsel van de uitvoer. Een GeoJSON-antwoord bevat geen vaste stelselaanduiding in de features. Geef crs= mee aan from_features, anders heeft de tabel geen stelsel en rekent to_crs niet.
• Cartesische bewerkingen. De GeoPandas-README zegt daarover:

GeoPandas geometry operations are cartesian.

Bron: README van GeoPandas, gelezen op 11 oktober 2026. In RD zijn oppervlakten en afstanden daardoor in vierkante meters en meters. In WGS84 zijn ze in graden, en dus niet bruikbaar. Reken eerst om naar RD, zoals het script doet met crs=RD.

Python is niet de plek voor grote bestanden die meerdere mensen tegelijk bevragen. Een script laadt data in het geheugen van één machine, terwijl een database zoals PostGIS indexen, rechten en gelijktijdige gebruikers regelt. Voor het publiceren van data als API is pygeoapi een logischer keuze, dat zelf in Python is geschreven maar als dienst draait.

Voor het visueel inrichten en controleren van kaarten is een bureaubladprogramma als QGIS geschikter. Voor bestandsconversie met een vast commando is GDAL korter en sneller in te stellen dan een eigen script.

Een script is ook een stuk software. Het heeft een eigenaar nodig die bij een nieuwe Python-versie of bibliotheekversie kan ingrijpen, en een vaste plek waar het draait. Voor wie dat niet intern wil regelen, bestaat API-koppelingen als dienst, en voor ontsluiting van data het onderwerp dataplatforms. TechinFlow levert die diensten zelf, dus weeg deze pagina mee als een bron met een belang.