Samenvatting
ik heb de R package rAHNextract gemaakt om snel en automatisch meerdere AHN hoogtes of AHN gebieden op te halen via R, inclusief point clouds. Deze post beschrijft wat ik gedaan heb en hoe het werkt.
Zie ook de Readme op GitHub en R documentatie voor een volledig uitleg hoe de package precies werkt. Mochten jullie er gebruiken van willen maken, ontvang ik graag jullie feedback als er opmerkingen, vragen of wensen zijn.
— —
Ik heb de package rAHNextract geschreven en gepubliceerd om in R automatisch individuele hoogtes of gebieden uit de AHN1, AHN2 of AHN3 te halen. Ook is het mogelijk om point clouds van je gebied op te halen.
Tijdens mijn GIMA master scriptie vorig jaar heb ik veel gezocht naar voorgeschreven oplossingen online. Deze zijn er wel hoewel je nog steeds veel handmatige stappen moet ondergaan, vooral als je meerdere hoogtes of specifiek gebied(en) nodig hebt. Daarom wilde ik daar een mooie en complete oplossing voor maken de je kunt integreren in een script.
Dit is een volledig pakket dat maximaal gebruik maakt van alle AHN geo services om gegevens te downloaden. Hoewel je R nodig hebt, is zover ik kan zoeken momenteel de makkelijkste open source manier (al zeg ik het zelf ;)) om snel meerdere individuele hoogte of eigen gebieden op te halen, zowel voor hoogtes als point clouds.
Deze package werkt voor alle beschikbare versies van de AHN, resoluties en type hoogtes die online via geos services beschikbaar zijn gemaakt door AHN en beschikbaar gesteld via PDOK om gegevens te downloaden. Voor de point clouds zijn de (uit)gefilterde beschikbaar voor AHN1 en AHN2 net als de AHN3 geclassificeerde puntenwolken.
Installatie
De package is in eerste instantie beschikbaar gesteld via Github. In R kun je hem installeren via devtools::install_github("jellest/rAHNextract"). Package devtools is vereist voor de instalatie en de package raster voor het gebruik van rAHNextract.
library(raster)
library(devtools)
devtools::install_github("Jellest/rAHNextract")
library(rAHNextract)
Hoogtes ophalen
WCS vs. Kaartbladen methode
Voor het ophalen van hoogtes van individuele punten of hele gebieden kan beide op twee manieren opgehaald worden: via de WCS methode of via de beschikbare kaartbladen van PDOK. Via de WCS methode haalt hij de data direct van het internet. Dit is een snelle methode voor als je maar een paar AHN hoogtes nodig hebt. Als je veel punten moet ophalen uit een klein gebied, kan het gebruik maken van de kaartbladen methode sneller zijn.
Zet de parameter sheets.method = FALSE in onderstaande functies (uitgezonderd point clouds) om de WCS methode te gebruiken. Dit is tevens standaard zo ingesteld voor individuele hoogtes.
Individuele hoogte(s)
ahn_point()
WCS methode
Bij het ophalen van individuele hoogtes via de WCS geo service, haalt hij de volledige pixel van je locatie op net als de volledige 9 omringende pixels. Dit zorgt ervoor dat de juiste hoogte wordt opgehaald. De uitkomst is een getal dat de (bilinear) hoogte op de locatie van je punt representeert.
ahn_point(X = 136550, Y = 456060)
#> [1] "Download raster image succeeded."
#> [1] "Intersecting raster. Getting elevation..."
#> [1] "Elevation of AHNelevation: 7.38 m."
#> [1] 7.38
Kaartblad(en) methode
Van je punt maakt hij een BBOX gebied van 9 pixels. Nadat de bijbehorende kaartblad(en) zijn gedownload, haalt hij de (bilinear) pixel waarde op van je opgegeven punt via een intersect: de hoogte. Vallen je 9 pixels precies op een grens? Dan download hij alle benodigde kaartbladen.
Voor beide methoden kun je de functie for loops of mapply() gebruiken om voor meerdere hoogtes tegelijkertijd op te halen via een lijst met X/Y coördinaten.
AHN raster gebied(en)
ahn_area()
Je kunt een punt opgeven met een radius afstand om een cirkel als gebied te maken, een BBOX gebied opgeven (via een radius of BBOX coördinaten), of een heel eigen geometrie inlezen. Zie de R documentatie hoe je dit precies moet doen.
WCS methode
Via de WCS methode maakt hij in eerste instantie een BBOX van je gebied en download hij de volledige pixels (van je ingestelde resolutie) waar je gebied in valt. Daarna clipt hij het tot jouw gekozen gebied (cirkel, bbox of eigen gebied). De uitkomst is een raster geotiff FLOAT32 bestand.
Utrecht_circleWCS <- ahn_area(name = "Utrecht circle", X = 136550, Y = 456060, resolution = 0.5,
radius = 100, output.dir = "C:/myProject")
#> [1] "Creating circle from radius input."
#> [1] "Destination directory of output AHN area: C:/myProject"
#> [1] "Download raster image succeeded."
#> [1] "C:/myProject/Utrechtcircle_100m_rAHN3_05m_DSM.tif"
plot(Utrecht_circleWCS, xlab = "RD X", ylab = "RD Y", main = "AHN Elevation (m)")
Kaartblad(en) methode
Via het ophalen van kaartbladen, download hij eerst alle kaartbladen die binnen je gebied vallen. Daarna wordt elk kaartblad verkleint tot het gedeelte van jouw gekozen gebied. Als laatst worden deze delen van kaartbladen samengevoegd tot een nieuw raster bestand. Dit betekent dat je uitkomst één .tif bestand is dat de kaartbladen heeft samengevoegd en geclipt tot jouw gekozen gebied (cirkel, bbox of eigen gebied). Je kunt ervoor kiezen om de volledige kaartbladen weer te verwijderen uit je geheugen of dat je deze ook wilt bewaren (voor vervolg analyses)
AHN puntwolken gebied(en)
ahn_pc()
Het is mogelijk om de puntwolken (point clouds) van je eigen gebied op te halen van de AHN1, AHN2 of AHN3. Via de paramater ‘gefilterd = TRUE’ geef je aan of je de gefilterde of uitgefilterde bestanden wilt (alleen van toepassing voor AHN1 en AHN2).
Voor het ophalen van puntwolken gebruikt hij de kaartbladen die beschikbaar zijn gesteld op de website van de PDOK. Je eindresultaat is een .laz bestand van je gekozen gebied. Het downloaden en het lezen van deze puntwolken kaartbladen kost wel veel geheugen en tijd in beslag. Het downloaden van de kaartbladen duurt even en het gebied vervolgens kleiner maken tot bijv. 100m2 duurt al minimaal 5 minuten voor de AHN3.
library(sf)
Utrecht.shp <- sf::st_read("C:/myProject/Utrecht_oudegracht.shp")
#> Reading layer `Utrecht_oudegracht' from data source `C:\myProject\Utrecht_oudegracht.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 1 feature and 1 field
#> geometry type: POLYGON
#> dimension: XY
#> bbox: xmin: 136422 ymin: 455924.1 xmax: 136762 ymax: 456208.3
#> proj4string: +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs
Utrecht_shape_pc <- ahn_pc(name = "Utrecht polygon pc", polygon = Utrecht.shp)
#> [1] "The AHN sheets are loaded from or downloaded in: C:/Users/jelle/Documents/R/rAHNextract/AHN_sheets/AHN3/PC"
#> [1] "Creating area from custom geometry."
#> [1] "Destination directory of output point clouds area: AHN_output"
#> [1] "Found 1 sheet(s) with name(s):"
#> [1] "31hz2"
#> [1] "Downloading point cloud sheets..."
#> [1] "https://download.pdok.nl/rws/ahn3/v1_0/laz/C_31HZ2.LAZ"
#> [1] "Filter string: -keep_xy 136421 455924 136762 456209"
#> [1] "writing .laz"
Documentatie en problemen
De documentatie in R bevat uitleg over de verschillende functies inclusief alle parameters.
Mochten er nog vragen, opmerkingen of nieuwe ideeën zijn, dan hoor ik ze graag van jullie. Problemen kunnen ook op mijn Github pagina geplaatst worden.
. Alleen bij de AHN 1 krijg ik waardes in de duizenden meters terwijl die bij AHN2/3 rond de 7 meter zit. Enig idee hoe dit zit, vermoedelijk doe ik iets fout?