31:T53c,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.32:T468,surfcanvas = { const scale = 4; const imgdata = getimgdata(surfraster_gh[0], surfraster_gh.width, surfraster_gh.height,litho_colorscale); const container = d3.select(DOM.element('div')); container.style("position", "relative"); const canvas = container.append('canvas').node(); const context = canvas.getContext('2d'); canvas.width = imgdata.width*scale; canvas.height = imgdata.height*scale; const svg = container.append('svg') .attr("width", canvas.width) .attr("height", canvas.height) .style("position", "absolute") .style("top", '0px') .style("left", '0px'); context.putImageData(imgdata, 0, 0); context.scale(scale,scale); context.drawImage(canvas, 0, 0); //do svg stuff here const mapscale = (extent_gh[2]-extent_gh[0]) / canvas.width; const path = d3.geoPath().projection(null); svg.append("path") .datum(spoortdata_new) .attr("d", path) .attr("fill","none") .attr("stroke-width","2") .attr("stroke","black") //.attr("transform", "translate("+ -extent_gh[0] +","+ -extent_gh[1]+")scale("+mapscale+")") return container.node(); }33:T505,depthcanvas = { const scale = 4; const imgdata = getimgdata(rasters[Math.floor(depth_selected*2)], rasters.width, rasters.height,litho_colorscale); const container = d3.select(DOM.element('div')); container.style("position", "relative"); const canvas = container.append('canvas').node(); const context = canvas.getContext('2d'); canvas.width = imgdata.width*scale; canvas.height = imgdata.height*scale; const svg = container.append('svg') .attr("width", canvas.width) .attr("height", canvas.height) .style("position", "absolute") .style("top", '0px') .style("left", '0px'); context.putImageData(imgdata, 0, 0); context.scale(scale,scale); context.drawImage(canvas, 0, 0); //do svg stuff here const mapscale = (extent_gh[2]-extent_gh[0]) / canvas.width; const path = d3.geoPath().projection(null); svg.append("path") .datum(spoortdata_new) .attr("d", path) .attr("fill","none") .attr("stroke-width","2") .attr("stroke","black") //.attr("transform", "translate("+ -extent_gh[0] +","+ -extent_gh[1]+")scale("+mapscale+")") svg.append("path").datum([{x:0,y:schnitt*scale},{x:canvas.width,y:schnitt*scale}]).attr('id','crossline').attr("d",crossline).attr('stroke','red') return container.node(); }34:T4d0,thickness_canvas = { const scale = 4; const colorscale = d3.scaleSequential().domain([0,20]).interpolator(d3.interpolateViridis); const imgdata = getimgdata(depthraster[soilselect.value], depthraster.width, depthraster.height,colorscale); const container = d3.select(DOM.element('div')); container.style("position", "relative"); const canvas = container.append('canvas').node(); const context = canvas.getContext('2d'); canvas.width = imgdata.width*scale; canvas.height = imgdata.height*scale; const svg = container.append('svg') .attr("width", canvas.width) .attr("height", canvas.height) .style("position", "absolute") .style("top", '0px') .style("left", '0px'); context.putImageData(imgdata, 0, 0); context.scale(scale,scale); context.drawImage(canvas, 0, 0); //do svg stuff here const mapscale = (extent_gh[2]-extent_gh[0]) / canvas.width; const path = d3.geoPath().projection(null); svg.append("path") .datum(spoortdata_new) .attr("d", path) .attr("fill","none") .attr("stroke-width","2") .attr("stroke","black") //.attr("transform", "translate("+ -extent_gh[0] +","+ -extent_gh[1]+")scale("+mapscale+")") return container.node(); }35:[{"id":"49e429144d02d70bbd62014ee642da2e21560554e926463915ae6aeee2ab1c91e40efd23b5a3bb7dec4afab9805f7912b3763239f6108435d1cb19edc92b9ba7","href":"https://static.observableusercontent.com/files/49e429144d02d70bbd62014ee642da2e21560554e926463915ae6aeee2ab1c91e40efd23b5a3bb7dec4afab9805f7912b3763239f6108435d1cb19edc92b9ba7","name":"spoor_alphen@3.geojson","type":"application/geo+json","lastModified":1640030642573.017,"size":1512392}] 2f:{"title":"De GEOTOP als GEOTIFF","cells":[{"id":0,"value":"# De GEOTOP als GEOTIFF","pinned":false,"hidden":false,"mode":"md"},{"id":13,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":764,"value":"## Voorbereiden van de data\n\nEenmaal binnengehaald en uitgepakt kunnen de .IMG files met behulp van GDAL vrij snel worden omgezet naar een COG.\n\nMaak een Virtuele Raster op basis van alle img files\n\n gdalbuildvrt -separate geotopMV.vrt /var/data/geotop/GeoTOP_v01r3/voxelmodel/GeoTOP_v01r3_doorsnede_lithoklasse_MV/grids/*.img\n \nZet 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.\n\n gdal_translate -of COG -co COMPRESS=DEFLATE -co COPY_SRC_OVERVIEWS=YES geotopMV.vrt geotopMV.tiff\n \nPlaats de COG in de cloud:\n\n gsutil cp ./geotopMV.tiff s3://some-s3-bucket/geotopMV.tiff\n \n\n","pinned":false,"hidden":false,"mode":"md"},{"id":6,"value":"## Bekijken van de geotop","pinned":false,"hidden":false,"mode":"md"},{"id":780,"value":"De hele geotop lithologie is nu online beschikbaar als geotiff en we gebruiken de geotiffjs library om de metadata in te laden:","pinned":false,"hidden":false,"mode":"md"},{"id":63,"value":"geotopurl = \"https://ahp-research.storage.googleapis.com/rasters/geotop_MV.tiff\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":616,"value":"geotoptiff = geotiff.fromUrl(geotopurl)","pinned":true,"hidden":false,"mode":"ojs"},{"id":787,"value":"We kunnen zien dat er 5 niveaus (1 origineel en 4 overviews) zijn opgeslagen:","pinned":false,"hidden":false,"mode":"md"},{"id":625,"value":"numlevels = geotoptiff.getImageCount()","pinned":true,"hidden":false,"mode":"ojs"},{"id":796,"value":"We kunnen al een getImage request doen zonder de rasters ook meteen uit te lezen. Hiermee kunnen we zien hoe groot het image is en hoeveel banden (samples) er in zitten.","pinned":false,"hidden":false,"mode":"md"},{"id":618,"value":"geotopimage = geotoptiff.getImage()","pinned":true,"hidden":false,"mode":"ojs"},{"id":648,"value":"numbands = geotopimage.getSamplesPerPixel()","pinned":true,"hidden":false,"mode":"ojs"},{"id":621,"value":"image_bbox = geotopimage.getBoundingBox()","pinned":true,"hidden":false,"mode":"ojs"},{"id":623,"value":"image_width = geotopimage.getWidth()","pinned":true,"hidden":false,"mode":"ojs"},{"id":1197,"value":"image_resolution = geotopimage.getResolution()[0]; //only x direction [0] because we assume cells are square","pinned":true,"hidden":false,"mode":"ojs"},{"id":800,"value":"Om een eerste beeld te krijgen van de data is het niet verstandig om de volledige dataset van 85Mb in te lezen. In plaats daarvan lezen we een kleiner overview (nummer 2) uit als thumbnail.","pinned":false,"hidden":false,"mode":"md"},{"id":628,"value":"geotoptthumbnail = geotoptiff.getImage(3)","pinned":true,"hidden":false,"mode":"ojs"},{"id":977,"value":"We lezen alleen de eerste band uit (eerste 50 cm vanaf maaiveld) met de optie `samples:[0]`","pinned":false,"hidden":false,"mode":"md"},{"id":861,"value":"surfraster = await geotoptthumbnail.readRasters({pool,samples:[0]})","pinned":true,"hidden":false,"mode":"ojs"},{"id":985,"value":"We kunnen eenvoudig uitvinden wat voor waardes er in de data zitten door Set te gebruiken. De 'geen data' waarde is in dit geval `255`","pinned":false,"hidden":false,"mode":"md"},{"id":886,"value":"[...new Set(surfraster[0])].sort();","pinned":true,"hidden":false,"mode":"ojs"},{"id":1004,"value":"Om de data vervolgens weer te geven als kaartje zetten we het om naar een 'Image'. Wat je ziet is de verspreiding van 9 verschillende lithologieen vlak onder het maaiveld. We zien veel zand rond de kust en stuwwal, klei rond de rivieren en plukjes veen er tussenin.","pinned":false,"hidden":false,"mode":"md"},{"id":1291,"value":"legendawaarden = [\n {name: 'Antropogeen',value: 0, color: 'grey'},\n {name: 'Veen',value: 1, color: '#9d4e40'},\n {name: 'Klei', value: 2, color: '#009200'},\n {name: 'Leem', value: 3, color: '#c2cf5c'},\n {name: 'Zand (fijn)', value: 5, color: '#ffff00'},\n {name: 'Zand (med)' , value: 6, color: '#f3e106'},\n {name: 'Zand (grof)', value: 7, color: '#e7c316'},\n {name: 'Grind' , value: 8, color: '#e09c2d'},\n {name: 'Schelpen' , value: 9, color: 'purple'}\n]","pinned":false,"hidden":false,"mode":"ojs"},{"id":1300,"value":"litho_colorscale = d3.scaleOrdinal().domain(legendawaarden.map(d=>d.value)).range(legendawaarden.map(d=>d.color));","pinned":false,"hidden":false,"mode":"ojs"},{"id":1095,"value":"geotop_nl = {\n //const colorscale = d3.scaleSequential().domain([0,d3.max(surfraster[0])]).interpolator(d3.interpolateViridis);\n const imgdata = getimgdata(surfraster[0], surfraster.width, surfraster.height,litho_colorscale);\n let c = DOM.context2d(imgdata.width, imgdata.height, 1);\n c.putImageData(imgdata, 0, 0);\n return c.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":1306,"value":"Plot.legend({\n color: {\n domain: legendawaarden.map(d=>d.name),\n range: litho_colorscale.range() \n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1185,"value":"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.\n","pinned":false,"hidden":false,"mode":"md"},{"id":1213,"value":"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.","pinned":false,"hidden":false,"mode":"md"},{"id":1189,"value":"extent_gh = [98981,455098,114010, 465324]","pinned":true,"hidden":false,"mode":"ojs"},{"id":1206,"value":"rwindow = [\n Math.floor(~~(extent_gh[0] - image_bbox[0]) / image_resolution), \n Math.floor(~~(image_bbox[3] - extent_gh[3]) / image_resolution),\n Math.ceil(~~(extent_gh[2] - image_bbox[0]) / image_resolution), \n Math.ceil(~~(image_bbox[3] - extent_gh[1]) / image_resolution) \n];","pinned":true,"hidden":false,"mode":"ojs"},{"id":1219,"value":"surfraster_gh = await geotopimage.readRasters({pool,samples:[0],window: rwindow})","pinned":true,"hidden":false,"mode":"ojs"},{"id":1402,"value":"Het resulterende raster zetten we weer om in een plaatje, waarbij we ter referentie het spoor bij Alphen aan de Rijn hebben ingetekend:","pinned":false,"hidden":false,"mode":"md"},{"id":1700,"value":"crossline = d3.line().x(d => (d.x)).y(d => (d.y))","pinned":true,"hidden":false,"mode":"ojs"},{"id":1224,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1399,"value":"Plot.legend({\n color: {\n domain: legendawaarden.map(d=>d.name),\n range: litho_colorscale.range() \n }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1011,"value":"## De diepte in\n\nHet 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). \n\nHiervoor hoeven we alleen maar het aantal samples dat we opvragen uit te breiden naar de gewenste diepte:","pinned":false,"hidden":false,"mode":"md"},{"id":640,"value":"rasters = await geotopimage.readRasters({pool,samples:[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19],window: rwindow});","pinned":true,"hidden":false,"mode":"ojs"},{"id":1447,"value":"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.","pinned":false,"hidden":false,"mode":"md"},{"id":1232,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":1237,"value":"viewof depth_selected = Inputs.range([0,9.5],{label: \"Diepte tov. maaiveld\", step: 0.5,value: 1})","pinned":false,"hidden":false,"mode":"ojs"},{"id":1334,"value":"## Doorsnedes","pinned":false,"hidden":false,"mode":"md"},{"id":1331,"value":"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.\n\nHet 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:","pinned":false,"hidden":false,"mode":"md"},{"id":1341,"value":"cross_section = {\n const durchschnitt = new Uint8Array(rasters.length * rasters.width);\n\tfor (let j = 0; jraster.map((d,j)=>{\n if ( d===i ) resultraster[j] = resultraster[j] + 1; //increase raster cell by 1 when soilvalue==1 found\n }));\n returnraster.push(resultraster);\n } \n return returnraster;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":816,"value":"Het resultaat hiervan zetten we om naar een 'image', waarbij we de kleuren meegeven ahv. de gevonden dikte.","pinned":false,"hidden":false,"mode":"md"},{"id":1020,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":1803,"value":"//Plot.legend({opacity: {}, color: })\nviewof thickness_legend = Plot.legend({\n color: {\n interpolate: d3.interpolateViridis,\n domain: [0, 10]\n },\n width: 350,\n ticks: 10,\n label: \"Dikte pakket (m) →\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":445,"value":"viewof soilselect = Inputs.select(legendawaarden, {label: \"Bodemtype\", format: x => x.name, value: legendawaarden[1]})","pinned":false,"hidden":false,"mode":"ojs"},{"id":855,"value":"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.","pinned":false,"hidden":false,"mode":"md"},{"id":1727,"value":"## Conclusie\n\nHoewel 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. \n\nIn 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.","pinned":false,"hidden":false,"mode":"md"},{"id":664,"value":"## Documentatie\n\nMeer informatie over de inhoud van de GEOTOP in: \n\ndinoloket_toelichtingmodellen_20160606_tno_2016_r10133_geotop_v1r3_english.pdf\n\nMeer informatie of Cloud Optimized Geotiffs via: \n\nhttps://www.cogeo.org/","pinned":false,"hidden":false,"mode":"md"},{"id":863,"value":"---","pinned":false,"hidden":false,"mode":"md"},{"id":1059,"value":"## Dependencies","pinned":false,"hidden":false,"mode":"md"},{"id":1786,"value":"width = 600","pinned":false,"hidden":false,"mode":"ojs"},{"id":264,"value":"height = Math.round( ((extent_gh[3] - extent_gh[1])/(extent_gh[2] - extent_gh[0])) * width)","pinned":false,"hidden":false,"mode":"ojs"},{"id":1503,"value":"pathGenerator = function(data,mapscale){\n const coords = data.geometry.coordinates;\n const newcoords = coords.map(c=>{\n return [\n (c[0]-extent_gh[0])/mapscale,\n (c[1]-extent_gh[1])/mapscale\n ];\n });\n return newcoords;\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":1119,"value":"function getimgdata(raster, width, height, colorscale) {\n const data = new Uint8ClampedArray(width * height * 4);\n for (let i = 0; i < raster.length; i++) {\n let color = d3.color(colorscale(raster[i]))|| {r:0,g:0,b:0};\n data[i * 4] = color.r;\n data[i * 4 + 1] = color.g;\n data[i * 4 + 2] = color.b;\n data[i * 4 + 3] = isNaN(raster[i]) || raster[i] == -3.4028234663852886e+38 || raster[i] == 3.4028234663852886e+38 || raster[i] == 255 ? 0 : 255;\n }\n return new ImageData(data, width, height);\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":11,"value":"//geotiff = require( \"geotiff.js/dist/geotiff.bundle.min.js\")\ngeotiff = require('geotiff@2.0.4/dist-browser/geotiff.js')","pinned":true,"hidden":false,"mode":"ojs"},{"id":713,"value":"pool = new geotiff.Pool();","pinned":true,"hidden":false,"mode":"ojs"},{"id":280,"value":"spoordata = FileAttachment(\"spoor_alphen@3.geojson\").json()","pinned":true,"hidden":false,"mode":"ojs"},{"id":1613,"value":"spoortdata_new = {\nconst mapscale = (extent_gh[2]-extent_gh[0]) / width;\n \nreturn {type:\"FeatureCollection\", features: \n spoordata.features.map(d=>{\n const newcoords = d.geometry.coordinates.map(c=>{\n return [\n (c[0]-extent_gh[0])/mapscale,\n -(c[1]-extent_gh[3])/mapscale\n ];\n });\n \n return {type: \"Feature\", geometry: {type: \"LineString\", coordinates:newcoords}};\n})\n}\n\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":83,"value":"d3 = require(\"d3@6\", \"d3-geo-projection@3\")","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}