Inleiding

Hieronder laten wij een script zien waarmee meerdere rasterbestanden - met de raster package in R - snel en eenvoudig (èn geheugenefficiënt) kunnen worden samengevoegd tot één nieuw rasterbestand. In dit voorbeeld maken we gebruik van een aantal kaartbladen van het Actueel Hoogtebestand Nederland (AHN), en wel van het AHN2 5 meter maaiveldraster.

Waar?

Via de PDOK viewer hebben wij 18 kaartbladen in het midden van het land uigezocht: het gebied in de driehoek tussen de A27, de A2 en de A15. Hier kronkelt de Linge lieflijk tussen de Lek en de Waal. En ja, het geoefende oog zal ook het Geofort op de kaart ontwaren.

De 18 bladnummers staan in dit tekstbestand. Het kaartbladnummer komt zowel terug in de URL van de downloadlocatie, als in de naam van het zip- en het tifbestand.

De volgende stappen worden geautomatiseerd:

Het volledige script staat hier.

Werkwijze

Laten we beginnen:

Het script wordt hieronder regel voor regel toegelicht.

Eerst moet de benodigde bibliotheek raster geladen worden:

# Het laden van de benodigde package
library(raster)

Vervolgens lezen we het invoerbestand uit. De inhoud wordt opgeslagen in wat in R een vector heet, een soort array, een reeks waarden.

# Het inlezen van de kaartbladnummers in de vector kaartbladen
kaartbladen <- readLines("kaartbladen.txt")
# Bekijk de inhoud van de vector kaartbladen (ter informatie)
kaartbladen
##  [1] "38ez1" "38ez2" "38fz1" "38fz2" "39az1" "39az2" "38gn1" "38gn2"
##  [9] "38gz1" "38gz2" "38hn1" "38hn2" "38hz1" "38hz2" "39cn1" "39cn2"
## [17] "39cz1" "39cz2"

Op basis van de waarden in ‘kaartbladen’ kunnen we de namen van de tifbestanden aanmaken:

# Het aanmaken van de vector bestanden.tif
bestanden.tif <- sprintf("ahn2_5_%s.tif", kaartbladen)
# Bekijk de inhoud van de vector bestanden.tif (ter informatie)
bestanden.tif
##  [1] "ahn2_5_38ez1.tif" "ahn2_5_38ez2.tif" "ahn2_5_38fz1.tif"
##  [4] "ahn2_5_38fz2.tif" "ahn2_5_39az1.tif" "ahn2_5_39az2.tif"
##  [7] "ahn2_5_38gn1.tif" "ahn2_5_38gn2.tif" "ahn2_5_38gz1.tif"
## [10] "ahn2_5_38gz2.tif" "ahn2_5_38hn1.tif" "ahn2_5_38hn2.tif"
## [13] "ahn2_5_38hz1.tif" "ahn2_5_38hz2.tif" "ahn2_5_39cn1.tif"
## [16] "ahn2_5_39cn2.tif" "ahn2_5_39cz1.tif" "ahn2_5_39cz2.tif"

En dan de URL’s voor het daadwerkelijk downloaden:

# Het aanmaken van de vector urls
urls <- sprintf("http://geodata.nationaalgeoregister.nl/ahn2/extract/ahn2_5m/%s.zip", bestanden.tif)

Bekijk desgewenst de lijst met webadressen:

# Bekijk de inhoud van de vector urls (ter informatie)
urls

Definiëer een functie voor het downloaden. Deze functie retourneert de naam van het zipbestand (d.w.z. het laatste deel - de bestandsnaam - van de URL).

# Deze functie download een bestand (als dat aanwezig is)
download.indienmogelijk <- function(url, refetch=FALSE, path=".") {
  dest <- file.path(basename(url))
  if (refetch || !file.exists(dest))
    download.file(url, dest)
  dest
}

Pas deze functie toe op alle URL’s: het downloaden gaat beginnen!

# Pas de functie download.indienmogelijk toe op alle urls
# en vul meteen de vector bestanden.zip
bestanden.zip <- sapply(urls, download.indienmogelijk)

Kijk maar in de werkmap: daar staan nu 18 zipbestanden. Deze bestanden (iets meer dan 2 MB per stuk) zijn samen ongeveer 36 MB groot. En er is in de voorgaande stap meteen een vector aangemaakt met de namen van de zipbestanden:

# Bekijk de inhoud van de vector bestanden.zip (ter informatie)
# en check - optioneel - de aanwezigheid van de gedownloade bestanden
bestanden.zip

Pak de zipbestanden uit:

# Pas de functie unzip toe op alle zipbestanden
lapply(bestanden.zip, unzip)

Kijk maar in de werkmap: daar staan nu 18 tifbestanden. Deze bestanden (bijna 5 MB per stuk) zijn samen ongeveer 72 MB groot.

Verwijder nu de zipbestanden (dat scheelt weer schijfruimte):

# Met unlink wordt een bestand weggegooid - de zipbestanden zijn niet meer nodig
lapply(bestanden.zip, unlink)

Laad de rasterbestanden in R:

# Laad alle rasters in een lijst
input.rasters <- lapply(bestanden.tif, raster)
names(input.rasters)[1:2] <- c('x', 'y')

Bedenk een mooie naam voor het uitvoerbestaand. Als je deze naam meegeeft bij het samenvoegen, dan kan R het rasterbestand meteen opslaan. (Dus dat hoef je dan niet achteraf zelf te doen…)

# Geef een naam waarmee het rasterbestand direct op schijf opgeslagen kan worden
# (anders moet het volledig in het geheugen geladen worden, met alle gevolgen vandien...)
input.rasters$filename <- "rastertotaal.tif"

En dan volgt het eigenlijke samenvoegen:

# Voeg de rasters samen
raster.totaal <- do.call(merge, input.rasters)

Kijk maar in de werkmap: daar staat nu het nieuwe bestand rastertotaal.tif, met een grootte van ongeveer 60 MB. Bekijk de beschrijving van het bestand in R:

# Toon de beschrijving
raster.totaal
## class       : RasterLayer 
## dimensions  : 3750, 6000, 22500000  (nrow, ncol, ncell)
## resolution  : 5, 5  (x, y)
## extent      : 120000, 150000, 425000, 443750  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +units=m +no_defs 
## data source : C:\RasterImport\rastertotaal.tif 
## names       : rastertotaal 
## values      : -4.308, 19.916  (min, max)

De databron is het bestand op schijf - daarom is er slechts weinig geheugenruimte nodig om het object in R op te slaan. Dit kun je nagaan met de functie object.size; deze geeft een schatting van het gebruikte geheugen:

# Toon de gebruikte geheugenruimte
object.size(raster.totaal)
## 11792 bytes

Toon het raster in het plotscherm:

# Plot de kaart
plot(raster.totaal)

(Deze afbeelding staat bovenaan de pagina.)

En voeg desgewenst een laag met gemeentegrenzen toe - ter oriëntatie:

plot(gem2014, add=TRUE, border="blue")

(Voor meer informatie over, en het bronbestand van de laag 'gem2014', zie deze blog: readOGR(): het inlezen van geografische data in R)


En ten slotte kunnen we nog even een beetje opruimen:

# Ruim op:
# Verwijder alle oorspronkelijke tifbestanden
lapply(bestanden.tif, unlink)
# Verwijder alle variabelen en functies die niet meer nodig zijn uit het geheugen van R
rm(kaartbladen, bestanden.tif, urls, download.indienmogelijk, bestanden.zip, input.rasters)

Omdat de raster() module zeer efficiënt met het geheugen omgaat, kan er via bovenstaande werkwijze - gegeven voldoende schijfruimte - een veel groter aantal rasterbestanden samengevoegd worden, dan het luttele aantal dat in deze test is gebruikt.



Egge-Jan Pollé