Het Actueel Hoogtebestand Nederland (AHN) voorziet ons al 20 jaar lang van puntenwolken van Nederland, inmiddels in versie 5. Maar Caribisch Nederland ontbrak tot nu toe. Aan dit gemis is nu een einde gekomen – Bonaire, Sint Eustatius en Saba zijn de afgelopen twee jaar ingewonnen en de data is nu voor iedereen beschikbaar. Een uitgelezen kans dus om niet alleen van de tropen te dromen, maar met deze data ook twee kaarten te maken. Hiervoor gebruik ik de Generic Mapping Tools.
Het AHN van Caribisch Nederland wordt in tegels van 1x1km geleverd, was aanzienlijk aangenamer is om mee te werken dan de gewone AHN units van 5×6,25km. De puntenwolken zijn zoals gewoonlijk geclassificeerd, met klasse 2 voor maaiveld en klasse 1 voor overig, waaronder vegetatie. De eerste stap is het splitsen van maaiveld en vegetatie. Voor de verwerking met GMT is het handig om daarbij meteen een conversie naar ASCII uit te voeren en een downsampling toe te passen, wat met de LAStools eenvoudig te doen is:
las2txt -i CN2023_C_07000_32000.LAZ -o ground.xyz -keep_class 2 -keep_every_nth 10 -parse xyz las2txt -i CN2023_C_07000_32000.LAZ -o veg.xyz -keep_class 1 -keep_every_nth 10 -parse xyz
Hiermee wil ik twee kaarten maken:
- Een topografische kaart met hoogtelijnen.
- Een kaart met vegetatiehoogtes.
Topografische kaart
gmt begin topo png
gmt info ground.xyz -I1 > region.txt
set /p region=<region.txt
gmt blockmedian ground.xyz %region% -I10 > median.xyz
gmt surface median.xyz %region% -Gground.nc -I1
gmt grd2cpt ground.nc -Crelief -Z
gmt basemap %region% -JX15c -B100
gmt grdimage ground.nc -I+a0
gmt colorbar -DJRM+o1c/0+mc -I0.3 -Bx20+lHoogte -By+lm
gmt grdcontour ground.nc -C10 -L0/1000 -A20 -Gd5c
gmt end show
del region.txt median.xyz ground.nc
Het batch script doet de volgende dingen in de “nieuwe” GMT mode:
- Het bestand zal topo heten en in PNG formaat opgeslagen worden.
- De omtrek van het invoerbestand wordt bepaald.
- Er wordt een mediaanfilter toegepast om een mooi resultaat te krijgen.
- Uit de gefilterde hoogtes wordt een grid met 1m resolutie berekend. Eventuele gaten in de data worden daarbij dichtgerekend. Het standaard gridformaat van GMT is NetCDF, deze kan ook in QGIS geopend worden.
- Het grid wordt gebruikt om een kleurenschaal aan te maken.
- De basemap wordt aangemaakt, deze legt afmetingen en layout zoals de annotaties vast.
- Het grid wordt geplot.
- De kleurenschaal wordt toegevoegd.
- De contourlijnen worden berekend en geplot.
- De plotfunctie wordt beëindigd en de plot getoond.
- Tijdelijke bestanden worden gewist.
GMT biedt daarbij ontzettend veel mogelijkheden voor de opmaak van de plots en is ideaal voor het automatiseren van plots (kijk maar eens naar naar de Actuele waarnemingen van het KNMI). Het zou weinig moeite kosten om op deze manier met behulp van een for loop voor alle AHN tegels een plot te maken.
Vegetatiehoogtes
gmt begin veg png
gmt info ground.xyz -I1 > region.txt
set /p region=<region.txt
gmt blockmedian ground.xyz %region% -I1 > median.xyz
gmt nearneighbor median.xyz %region% -Gground.nc -I2 -S4
gmt blockmedian veg.xyz %region% -I1 > median.xyz
gmt nearneighbor median.xyz %region% -Gtemp.nc -I2 -S4
gmt grdfilter temp.nc -Gveg.nc -D0 -Fu5
gmt grdmath veg.nc ground.nc SUB = vegheight.nc
gmt grd2cpt vegheight.nc -Chaxby -I -D -E -Z
gmt basemap %region% -JX15c -B100
gmt grdimage vegheight.nc
gmt colorbar -DJRM+o1c/0+mc -I0.3 -Bx1+lVegetatiehoogte -By+lm
gmt end show
del region.txt median.xyz ground.nc veg.nc vegheight.nc temp.nc
Er zijn aanvullende stappen om de vegetatiehoogte te berekenen:
- Berekenen maaiveldgrid.
- Berekenen ruw vegetatiegrid.
- Filtering hiervan, waarbij binnen een zoekstraal van 5m de hoogste waarde gebruikt wordt.
- Verschilberekening.
- Plot van het resulterende grid.
De resulterende plot laat zijn dat de vegetatie grotendeels rond de 4m hoog is, in sommige gevallen tot 7m. Er zijn uitschieters tot 17m aanwezig, maar deze kunnen ook het gevolg van grote hoogteverschillen in het terrein zijn. Ook wegen worden duidelijk zichtbaar.
Meer GMT
Dit zijn maar twee voorbeelden van wat er mogelijk is met de Generic Mapping Tools. De voorbeelden op de GMT website laten nog veel meer mogelijkheden zien.
