In Nederland hebben we met het Actueel Hoogtebestand Nederland een schitterend landelijk hoogtebestand. Schitterend omdat het open data is, van goede kwaliteit (5cm stochastisch + 5cm systematisch), hoge resolutie voor een airborne LiDAR bestand, geclassificeerd en regelmatig updates krijgt.
Verschillen tussen deze updates berekenen is bijzonder interessant in een land waar maaivelddaling een belangrijk onderwerp is. Hierbij is het echter belangrijk om rekening te houden met de kwaliteit van het AHN – domweg twee jaargangen van elkaar aftrekken en het resultaat als exacte waarheid te interpreteren is onjuist. Ik heb nu al een paar keer een mooie analyse van Geodelta hierover langs zien komen en was benieuwd of ik die resultaten kon reproduceren en misschien eenvoudig een correctie kan toepassen.
Hiervoor heb ik de 0,5m DTM (alleen maaiveld) en DSM (alle punten) grids van AHN3 en AHN4 gedownload, van de unit 46CN2. Mijn favoriete toolset voor het werken met grids zijn al meer dan 20 jaar de Generic Mapping Tools, en die heb ik dan ook hiervoor gebruikt. GMT bevat veel verschillende command line tools voor het verwerken en visualiseren van geografische data en is uitstekend geschikt voor batch processing. De plaatjes van de actuele weergegevens van het KNMI lijken bijv. met GMT gegenereerd te zijn.
Verschilberekening
De AHN grids worden als GeoTIFF beschikbaar gesteld, GMT werkt liever met het NetCDF formaat. De eerste stap is dus de conversie hiervan:
gmt grdconvert AHN4\M_46CN2.TIF=gd:GTIFF -GAHN4\M_46CN2.nc
gmt grdconvert AHN3\M_46CN2.TIF=gd:GTIFF -GAHN3\M_46CN2.nc
Daarnaast bleek het ene grid gridline registration te gebruiken en de andere pixel registration, maar ook dit is makkelijk op te lossen:
gmt grdedit -T AHN3\M_46CN2.nc
Dan kunnen de verschillen berekend worden:
gmt grdmath AHN4\M_46CN2.nc AHN3\M_46CN2.nc SUB = diff_ahn_34.nc
De modern mode van GMT is nog een beetje wennen voor mij, maar zo kun je het verschilgrid als png visualiseren:
gmt begin diff_ahn_34 png
gmt makecpt -Cpolar -D -T-0.1/0.1/0.01 -Z -I
gmt grdimage diff_ahn_34.nc -JX12c
gmt colorbar -Dx6c/-0.5c+w10c/0.5c+jTC+h -Bxaf+l"Hoogteverschil [m]"
gmt end show
Hier het resultaat (klik voor volle resolutie):
Diagonaal is een strook met sterke “maaivelddaling” te zien die echter duidelijk is gerelateerd aan de vliegrichting. Bijzonder zijn ook de blokachtige artefacten in de noord-west-hoek van de unit die mogelijk gerelateerd zijn aan een bij de verwerking door de inwinner gebruikte blokindeling.
Berekenen correcties
Om wel tot een gedegen uitspraak over de verschillen tussen AHN3 en AHN4 te komen is het dus noodzakelijk om te corrigeren voor de onnauwkeurigheden van beide bestanden. De juiste geodetische aanpak is hiervoor het zoeken van vaste punten en het hieruit bepalen van geschikte correcties. Hiervoor zou je een vereffeningsmodel moeten opstellen zodat via toetsing de instabiele punten gedetecteerd en geëlimineerd kunnen worden. Ik wil hier echter een simpelere methode uitproberen.
Voor de hand liggend is het gebruik van (hopelijk) stabiele gebouwen als vaste punten. Deze zijn in de AHN DSM grid bestanden aanwezig.
Om deze te kunnen gebruiken moet eerst een masker aangemaakt worden waarmee de gebouwen geëxtraheerd kunnen worden. De AHN puntenwolken zijn geclassificeerd, ik heb dus de AHN4 puntenwolk gedownload en met de LAStools de gebouwpunten in een apart bestand gezet:
las2las -i AHN4\C_46CN2.LAZ -o buildings.laz -keep_class 6
Wie over een LAStools licentie beschikt kan met lasgrid hieruit een grid berekenen, maar met de omweg via een tekstbestand gaat dit ook met GMT. Met behulp van grdclip worden alle gridcellen die een waarde hebben op 1 gezet, de overige (dus waar geen gebouw is) krijgen de waarde NaN (not a number). Het resulterende masker kan dan vermenigvuldigd worden met het verschil van de DSM grids.
las2txt -i buildings.laz -o buildings.xyz -parse xyz
gmt nearneighbor buildings.xyz -Gbuildings.nc -I0.5 -Rdiff_ahn_34.nc -S0.5
gmt grdclip buildings.nc -Si-1000/1000/1 -Gmask.nc
gmt grdconvert AHN4\R_46CN2.TIF=gd:GTIFF -GAHN4\R_46CN2.nc
gmt grdconvert AHN3\R_46CN2.TIF=gd:GTIFF -GAHN3\R_46CN2.nc
gmt grdedit -T AHN3\R_46CN2.nc
gmt grdmath AHN4\R_46CN2.nc AHN3\R_46CN2.nc SUB mask.nc MUL = diff_ahn_34_buildings.nc
gmt begin diff_ahn_34_buildings png
gmt makecpt -Cpolar -D -T-0.1/0.1/0.01 -Z -I
gmt grdimage diff_ahn_34_buildings.nc -JX12c
gmt colorbar -Dx6c/-0.5c+w10c/0.5c+jTC+h -Bxaf+l"Hoogteverschil [m]"
gmt end show
Nu beginnen de echte problemen:
- Gebouwen die in het AHN4 aanwezig waren maar in het AHN3 niet resulteren, net als andere mutaties, in grote verschillen.
- Punten op de gevel kunnen in gridhoogtes resulteren die niet met de dakhoogte overeenkomen
- Kleine verschillen in X en Y tussen AHN3 en AHN4 kunnen bij zadeldaken tot hoogteverschillen leiden
De laatste twee punten zouden opgelost kunnen worden door de vergelijking te beperken tot platte daken, wat mogelijk zou zijn door de normaalvector van de punten te berekenen en hierop te filteren. Voor deze analyse heb ik in eerste instantie een andere aanpak gekozen: gewoon over een groter gebied middelen. GMT heeft hiervoor in het geval van tekstbestanden drie opties:
- Met blockmean kun je het aritmetisch gemiddelde berekenen
- Met blockmedian kun je de mediaan berekenen, een robuust criterium
- Met blockmode kun je de mode berekenen – dat is de waarde die het meest voorkomt, dus het histogrambakje met de meeste waarden
Deze verwachten wel allemaal losse punten als input, dus het grid moet eerst naar een tekstbestand omgezet worden. Na het filteren kun je met surface het correctiegrid berekenen, waarbij de T optie bepaald hoe strak het vlak wordt getrokken. De surface berekening heeft even tijd nodig.
gmt begin corrgrid_median png
gmt grd2xyz diff_ahn_34_buildings.nc -s > diff_buildings.xyz
gmt blockmedian diff_buildings.xyz -I100 -Rdiff_ahn_34_buildings.nc > median.xyz
gmt surface median.xyz -Gcorr_median.nc -Rdiff_ahn_34_buildings.nc -I0.5 -T0.35
gmt makecpt -Cpolar -D -T-0.1/0.1/0.01 -Z -I
gmt grdimage corr_median.nc -JX12c
gmt colorbar -Dx6c/-0.5c+w10c/0.5c+jTC+h -Bxaf+l"Correctie [m]"
gmt end show
Mijn eerste poging met 100m filterradius was duidelijk niet genoeg gefilterd:
Met 500m werd het ook niet beter, dus het is duidelijk dat de verschillen eerst beter gefilterd moeten worden. Met grdclip kunnen gridcellen die buiten een bepaald bereik liggen op NaN gezet worden.
gmt grdclip diff_ahn_34_buildings.nc -Gdiff_ahn_34_buildings_clip.nc -Sa0.15/NaN -Sb-0.15/NaN
Dan wordt het een stuk beter:
En hier met blockmode:
Toepassen correcties
Als we een grid hebben waarmee we kunnen leven kan deze met grdmath op de oorspronkelijke verschillen toegepast en gevisualiseerd worden:
gmt begin diff_ahn_34_corr png
gmt makecpt -Cpolar -D -T-0.1/0.1/0.01 -Z -I
gmt grdmath diff_ahn_34.nc corr_mode.nc SUB = diff_ahn_34_corr.nc
gmt grdimage diff_ahn_34_corr.nc -JX12c
gmt colorbar -Dx6c/-0.5c+w10c/0.5c+jTC+h -Bxaf+l"Hoogteverschil [m]"
gmt end show
Hier het resultaat naast het oorspronkelijke verschilplaatje:
Zoals je kunt zien is het gelukt om de stelselmatige verschillen deels te verwijderen, maar aan de oostkant is ook een nieuw fout geïntroduceerd.
Dag zadeldaken
Zadeldaken kun je detecteren door de normaalvector van de punten berekenen, maar GMT biedt nog een tweede optie: met grdgradient kun je de gradient (dus de helling van het gridbestand) berekenen en hieruit met grdclip een tweede masker berekenen.
gmt grdgradient buildings.nc -Ggradient.nc -A0/270
gmt grdclip gradient.nc -Sa0.1/NaN -Sb-0.1/NaN -Si-0.1/0.1/1 -Gmask2.nc
gmt grdmath mask.nc mask2.nc MUL = mask3.nc
Hiermee kun je dan weer zoals boven de berekende verschillen filteren. Het resulterende correctiegrid laat nu veel beter het streeppatroon zien:
En hier is dan het gecorrigeerde verschilgrid:
Conclusie
Is deze simpele methode een zaligmakende oplossing voor het elimineren van stelselmatige verschillen tussen AHN3 en AHN4? Nee, zeker niet. Je eigenlijk een vereffeningsmodel willen gebruiken dat a) rekening houdt met de oorspronkelijke vliegstrookgeometrie en b) een toetsing mogelijk maakt.
Desalniettemin zou je bij zorgvuldig gebruik deze aanpak kunnen gebruiken om puntenwolken uit bijv. dronevluchten op het AHN in te passen of binnen een gemeente stelselmatige verschillen te elimineren. Een bijkomstig voordeel is dat volledig gebruik gemaakt wordt van open source software, dat het goed te automatiseren is en de resultaten makkelijk te visualiseren zijn, ook in QGIS.
