De geotop van Nederland is onderdeel van de 'BasisRegistratie Ondergrond' (BRO). Het is een bestand waarin de bodemlagen van een groot deel van Nederland tot een diepte van 50 meter worden gegeven. Grofweg zijn er twee typen data: stratigrafische data (wanneer en hoe zijn de lagen afgezet) en lithologische data (van wat voor materiaal is het). In dit voorbeeld ga ik in op de lithologische data.
De geotop is hier te downloaden van het PDOK. Je vind hier een behoorlijk uitgebreide set aan data maar wij concentreren ons hier op het voxelmodel met doorsnedes van de lithoklassen. In de praktijk zijn dit 200 rasterbestanden met waardes van 1 tot 9 die de lithologie beschrijven. Tezamen is dat 2.9 Gigabyte aan bestanden.
Hoewel de data prima bruikbaar is in lokale analyses op goede computers is het niet zo geschikt om online weer te geven. Om te onderzoeken wat de mogelijkheden zijn om dit wel te kunnen doen heb ik de data omgezet naar een Cloud Optimized Geotiff (COG) waar alle data bij elkaar zit en welke te raadplegen is via een simpele http interface.
Het resultaat van dit korte onderzoek heb ik hier beschreven.
Voorbereiden van de data
Eenmaal binnengehaald en uitgepakt kunnen de .IMG files met behulp van GDAL vrij snel worden omgezet naar een COG.
Maak een Virtuele Raster op basis van alle img files
gdalbuildvrt -separate geotopMV.vrt /var/data/geotop/GeoTOP_v01r3/voxelmodel/GeoTOP_v01r3_doorsnede_lithoklasse_MV/grids/*.img
Zet de Virtuele Raster om naar een grote Cloud-Optimized-Geotiff (COG). Merk op dat de oorspronkelijke IMG bestanden zeer goed comprimeerbaar zijn, uiteindelijk houden we een COG bestand van 85Mb aan data over.
gdal_translate -of COG -co COMPRESS=DEFLATE -co COPY_SRC_OVERVIEWS=YES geotopMV.vrt geotopMV.tiff
Plaats de COG in de cloud:
gsutil cp ./geotopMV.tiff s3://some-s3-bucket/geotopMV.tiff
Bekijken van de geotop
De hele geotop lithologie is nu online beschikbaar als geotiff en we gebruiken de geotiffjs library om de metadata in te laden:
Zoals je ziet zit niet heel Nederland in de geotop, het zijn vooral de rivier en kustgebieden die zijn opgenomen. Het plaatje wordt interessanter als we wat inzoomen op zo'n stuk. In dit geval gaan we kijken naar het groene hart, waar een mix zit van veen, klei en zanderige rivierbeddingen.
Hiervoor moeten we eerst de 'extent' van het groene hart in rijksdriekhoekcoordinaten omrekenen naar de 'extent' in pixelwaarden van de geotiff. Dit is goed te doen omdat we de volledige extent van het image weten en de schaal van de pixels.
Het resulterende raster zetten we weer om in een plaatje, waarbij we ter referentie het spoor bij Alphen aan de Rijn hebben ingetekend:
De diepte in
Het wordt natuurlijk pas interessant als we meerdere lagen van de geotop gaan bekijken. In totaal zitten er 100 lagen in de geotop (dat is 50 meter want elke laag beschrijft een halve meter sediment) maar voor nu zijn we alleen geinteresseerd in het bovenste gedeelte, de eerste 10 meter (dus de eerste 20 banden).
Vervolgens kunnen we voor elk van deze dieptes het bijbehorende plaatje tekenen. Gebruik de slider onder het plaatje om snel door de lagen heen te gaan. De rode lijn in het plaatje is de locatie van de doorsnede waar we hieronder op terugkomen.
Als je een echt goed beeld wilt hebben van wat er in de bodem zit dan is een doorsnede van de data noodzakelijk. Tot nu toe hebben we gekeken naar 'horizontale' plakjes geotop maar voor een doorsnede zullen we daar verticale plakjes van moeten maken.
Het simpelste is om een rechte lijn van west naar oost te trekken voor de doorsnede, omdat de pixelwaarden in het raster ook van links naar rechts lopen. Je hoeft voor elke diepte rasters.length alleen maar de juiste reeks uit het raster te vissen.
De dikte van het pakket
Omdat we ook geinteresseerd zijn in de totale dikte van bepaalde grondsoorten dan alleen de samenstelling van een bepaalde laag of doorsnede gaan we de eerste 10 meter aan data bij elkaar optellen. We maken hiervoor een raster met 8 nieuwe banden aan waarin voor elke grondsoort de totale dikte zit. Hiervoor lopen we door de 20 rasterset heen en tellen we 0.5 meter bij de totale dikte als we de gezochte lithologie tegenkomen.
Het resultaat hiervan zetten we om naar een 'image', waarbij we de kleuren meegeven ahv. de gevonden dikte.
In het resultaat zien we een beeld dat we wel zo'n beetje verwachten. Rond de Oude Rijn vooral zand, de Veenpakketten liggen een stukje van de oude Rijn af en daaromheen de polder waar het veen is weggehaald en alleen nog klei over is.
Conclusie
Hoewel bovenstaande voorbeelden relatief simpel zijn en ook al vaker uitgevoerd, is dit naar mijn weten de eerste keer dat de geotop als 'cloud optimized geotiff' wordt ingezet. Het voordeel hiervan is dat de data eenvoudig online beschikbaar is, relatief klein van formaat (60mb) en, dankzij http2 range requests, eenvoudig in nog kleinere stukjes op te vragen zodat de uiteindelijke analyse toe kan met een minimale hoeveelheid bytes.
In een volgende post zal ik wat meer voorbeelden geven van de praktische toepassing, onder meer door te kijken naar de dataset op werkelijke hoogte in combinatie met het AHN en of er een verband zit tussen de lithologie en gemeten verzakking op bepaalde plekken.
Documentatie
Meer informatie over de inhoud van de GEOTOP in:
dinoloket_toelichtingmodellen_20160606_tno_2016_r10133_geotop_v1r3_english.pdf
Meer informatie of Cloud Optimized Geotiffs via: