GDAL en ogr2ogr gebruiken met BAG en AHN

GDAL is een vertaalbibliotheek voor ruimtelijke data. Rasterformaten behandelt het onder de naam GDAL, vectorformaten onder de naam OGR, maar het is één project met één commandoregeltoolset. QGIS, PostGIS (voor rasters), MapServer en de Python-pakketten rasterio en Fiona gebruiken GDAL of de OGR-laag eronder. Wie een gemeentelijke dataset wil omzetten, controleren of publiceren, komt er vroeg of laat mee in aanraking.

De broncode valt onder een vrije licentie. De projectdocumentatie zegt hierover:

In general, the GDAL/OGR source code is licensed under the MIT license

Bron: GDAL-documentatie, pagina License. Dat betekent dat u GDAL zonder licentiekosten en ook in gesloten software mag gebruiken, zolang de licentietekst meegaat. Let wel: afzonderlijke drivers of afhankelijkheden (bijvoorbeeld voor bepaalde propriëtaire formaten) kunnen eigen voorwaarden hebben.

Controleer uw eigen versie met gdalinfo --version. Veel Linux-distributies leveren een oudere GDAL mee dan de actuele release; bepaalde drivers en opties in dit artikel (bijvoorbeeld het commando -if OAPIF) vragen een recente versie.

Met de OAPIF-driver van GDAL opent u de PDOK OGC API-endpoints als een gewone databron. De documentatie omschrijft de driver als volgt:

This driver can connect to a OGC API - Features service.

Bron: GDAL-documentatie, pagina OGC API - Features. De naamsyntax is OAPIF:<url van de landingspagina>. De BAG heeft bij PDOK zes collecties: adres, ligplaats, pand, standplaats, verblijfsobject en woonplaats. Zie ook de pagina over de BAG voor de achtergrond van de registratie.

Het volgende commando haalt de panden in een kader van 200 bij 150 meter rond het Stationsplein in Amersfoort op en schrijft ze naar een GeoPackage. Het kader staat in RD-meters, vandaar -spat_srs.

ogr2ogr -f GPKG bag.gpkg "OAPIF:https://api.pdok.nl/kadaster/bag/ogc/v2" pand \
  -spat 154100 462650 154300 462800 -spat_srs EPSG:28992 -limit 50

Dit duurde in twee tests (11 en 12 oktober 2026) tussen 1,6 en 4,5 seconden en leverde een bestand van 176 kB. Daarna controleert u het resultaat:

ogrinfo -so bag.gpkg pand

De relevante regels uit de uitvoer op 11 oktober 2026:

Geometry: Polygon
Feature Count: 32
Extent: (153984.979000, 462637.676000) - (154403.432000, 462972.915000)
...
aantal_verblijfsobjecten: Integer (0.0)
bouwjaar: Integer (0.0)
documentdatum: Date
gebruiksdoel: String (0.0)
identificatie: String (0.0)
status: String (0.0)
verblijfsobject.href: String(JSON) (0.0)

Twee dingen vallen op. Het aantal features (32) is lager dan de opgegeven -limit 50, dus het limiet bepaalde hier niet de uitkomst. De extent is groter dan het opgegeven kader, omdat een pand dat het kader raakt als geheel in de laag komt. De laag staat in RD New (EPSG:28992), het stelsel waarin de PDOK-collectie de geometrie aanleverde.

Na het ophalen werkt u lokaal verder, zonder nog de PDOK-service te belasten. Dat is ook de netste werkwijze richting de bron: haal eenmaal op, bewerk lokaal. De regels voor gebruik en caching staan op de pagina over PDOK fair use en caching.

Een SQL-query met het SQLite-dialect gebruikt de ruimtelijke functies van de GeoPackage. Zo vindt u de drie grootste panden in het kader, met oppervlakte in vierkante meters (RD is een metrisch stelsel, dus ST_Area geeft direct m²):

ogrinfo -ro -dialect SQLite \
  -sql "SELECT identificatie, status, bouwjaar, ROUND(ST_Area(geom)) AS m2 FROM pand ORDER BY m2 DESC LIMIT 3" bag.gpkg

De uitvoer, ingekort tot de eerste twee records:

identificatie (String) = 0307100000365462
status (String) = Pand in gebruik
bouwjaar (Integer) = 2003
m2 (Real) = 10099

identificatie (String) = 0307100000580175
status (String) = Pand in gebruik
bouwjaar (Integer) = 1998
m2 (Real) = 5646

Voor een webkaart zet u dezelfde laag om naar GeoJSON in WGS84. De optie RFC7946=YES zorgt dat de uitvoer de standaard volgt, inclusief de volgorde lengtegraad-breedtegraad en WGS84 als enig stelsel. Met -select bewaart u alleen de velden die de kaart nodig heeft, wat het bestand klein houdt:

ogr2ogr -f GeoJSON panden.geojson bag.gpkg pand -t_srs EPSG:4326 \
  -lco RFC7946=YES -select identificatie,bouwjaar,status

Het resultaat begint (opnieuw gedraaid op 12 oktober 2026 met GDAL 3.13.3) met {"type":"Feature","properties":{"bouwjaar":1986,"identificatie":"0307100000308303","status":"Pand in gebruik"},"geometry":{"type":"Polygon","coordinates":[[[5.3753709,52.1523555],.... Wilt u losse punten omrekenen in plaats van bestanden, gebruik dan de RD naar WGS84-tool. Voor grotere bewerkingen op veel data past een database beter bij het werk dan een bestand per bewerking.

Het AHN (Actueel Hoogtebestand Nederland) is te benaderen via een WCS van PDOK. De capabilities noemen de coverages dsm_05m (oppervlaktemodel) en dtm_05m (terreinmodel) met een rastercel van 0,5 meter. Als kostenregel vermeldt de service Fees: NONE. Achtergrond over de dataset staat op de pagina over het AHN.

Een WCS 2.0.1 GetCoverage-aanroep met SUBSET-parameters in RD-meters levert een GeoTIFF voor het opgegeven gebied. De test haalde hetzelfde kader als hierboven op:

curl -s -o ahn.tif "https://service.pdok.nl/rws/ahn/wcs/v1_0?SERVICE=WCS&VERSION=2.0.1&REQUEST=GetCoverage&COVERAGEID=dtm_05m&FORMAT=image/tiff&SUBSET=x(154100,154300)&SUBSET=y(462650,462800)"
gdalinfo -stats ahn.tif

Het resultaat is een raster van 400 bij 300 cellen. De relevante regels uit gdalinfo:

Size is 400, 300
Origin = (154100.000000000000000,462800.000000000000000)
Pixel Size = (0.500000000000000,-0.500000000000000)
Band 1 Block=400x5 Type=Float32, ColorInterp=Gray
  Minimum=5.094, Maximum=12.549, Mean=10.226, StdDev=1.374
  NoData Value=3.4028235e+38

De waarden zijn hoogtes in meters; de coverage dtm_05m bevat alleen terrein. Na -stats meldt gdalinfo ook STATISTICS_VALID_PERCENT=67.85: bijna een derde van de cellen in dit kader heeft geen waarde, omdat het terreinmodel onder gebouwen leeg is. Hoe u daar per pand omheen werkt, staat in het rasterio-voorbeeld bij Python voor GIS. Let op de NoData-waarde: dat is het grootste 32-bits getal, niet 0 of -9999. Software die deze waarde niet kent, rekent er gewoon mee en verknoeit gemiddelden en legenda's. gdalinfo -stats leest de NoData-waarde wel uit het bestand en houdt er rekening mee.

Daarna zet u het bestand om. Een Cloud Optimized GeoTIFF (COG) is geschikt om later op een webserver te zetten, een hillshade is een snelle visuele controle, en gdalwarp herprojecteert naar WGS84:

gdal_translate -of COG -co COMPRESS=DEFLATE ahn.tif ahn_cog.tif
gdaldem hillshade ahn.tif hs.tif -compute_edges
gdalwarp -t_srs EPSG:4326 -r bilinear ahn.tif ahn4326.tif

De COG-controle met gdalinfo ahn_cog.tif toont LAYOUT=COG. De hillshade kreeg in de test waarden van 1 tot en met 255. Na gdalwarp is de celgrootte uitgedrukt in graden: ongeveer 0,0000064 graden, in noord-zuidrichting circa 0,7 meter. Dat is een gevolg van de herprojectie; bij rasters blijft RD voor analyses de betere keuze, omdat de cellen daar gelijk blijven in meters.

De onderstaande tabel koppelt een taak aan het commando dat daarvoor bedoeld is. Alle commando's zijn opgenomen in de GDAL-documentatie onder Programs; de kolom met de test geeft aan wat op deze pagina is uitgevoerd.

De laatste twee regels zijn bewust als niet getest aangemerkt: ze staan in de officiële documentatie maar zijn hier niet gedraaid. Houd GDAL ook los van de vraag welke viewer u erbij kiest. GDAL bereidt data voor en zet haar om; een kaart voor eindgebruikers vraagt daarnaast een viewer.

De meeste problemen zitten niet in GDAL zelf maar in aannames over het coördinatenstelsel en over de omvang van het verzoek. De volgende punten kwamen in de tests naar voren of staan in de documentatie.

• Het coördinatenstelsel van het kader. Zonder -spat_srs leest GDAL het kader in het stelsel van de laag. Geef de kadercoördinaten daarom altijd samen met -spat_srs op, zodat RD-meters en WGS84-graden niet door elkaar lopen.
• Lijstvelden in GeoPackage. De BAG-collectie pand bevat het lijstveld verblijfsobject.href. De GPKG-driver kan dat niet native opslaan en zet het om naar JSON-tekst, met een waarschuwing. Wilt u dat anders, dan stuurt -mapFieldType de omzetting.
• NoData bij rasters. Zoals hierboven: 3.4028235e+38 is de NoData-waarde van de AHN-WCS. Ga niet van 0 of -9999 uit.
• Resolutie bij herprojectie. Een raster dat u naar graden omzet krijgt cellen in graden. Het ziet er op de kaart goed uit, maar oppervlakte- of afstandsberekeningen gaan dan fout.
• Bulkverzoeken aan PDOK. OAPIF pagineert het antwoord. Een groot gebied kost veel opeenvolgende verzoeken en kan tegen de gebruiksvoorwaarden aanlopen. Voor heel gemeentelijke of landelijke volumes biedt PDOK downloadservices en bestanden; zie het overzicht van open geodata en PDOK voor de services.
• Versieverschillen. De geünificeerde gdal-opdracht (bijvoorbeeld gdal raster info) is volgens de documentatie toegevoegd in versie 3.11; oudere installaties kennen alleen de losse programma's. Schrijf scripts voor de losse programma's als ze op verschillende machines moeten draaien.

GDAL is een batchhulpmiddel en geen database of kaartserver. Voor dagelijkse analyse met veel joins en filters op grote tabellen werkt u beter in een database; de aanpak staat bij PostGIS. Voor het publiceren van data aan anderen als dienst gebruikt u een server, zoals GeoServer of pygeoapi, die voor rasterformaten kunnen steunen op GDAL. Voor interactieve verkenning door een medewerker is QGIS handiger dan een commandoregel.

Voor gemeenten die hun beheer van kaartdata willen inrichten, onder andere de periodieke omzetting van bronbestanden naar een publicatievorm, staat de aanpak op de pagina GIS- en kaartoplossingen. Een GDAL-script in een geplande taak is een gebruikelijk bouwblok voor zo'n keten, maar het hoort bij een beheerde omgeving met logging en controle op de uitkomst, niet op een enkele laptop.