9. Raster data in R

Arbeitsverzeichnis, Daten und Pakete

Die Daten für dieses Tutorial werden über GitHub bereitgestellt Github. Um den Code mit Ihren eigenen Daten nachzubilden, ersetzen Sie die URL durch Ihren lokalen Dateipfad

Wenn Sie lokal arbeiten, können Sie das Arbeitsverzeichnis entsprechnd anpassen

setwd("YOUR/FILEPATH") 
# this will not work. if you want to use it, define a path on your computer

Dies sollte ihre Arbeitsumgebung darstellen, d.h. der Ort an dem ihre Ordner und dazugehörigen Daten abgespeichert sind. Weitere Informationen finden Sie unter hier

Codeabschnitte zum Speichern von Daten und Ergebnissen sind in diesem Tutorial (mit einer #) gekennzeichnet. Bitte ersetzen Sie (“YOUR/FILEPATH”) durch Ihren eigenen lokalen Dateipfad.

Der Code ist so geschrieben, dass Pakete automatisch installiert (wenn nicht bereits installiert) und geladen werden. Die Pakete werden am Standardspeicherort Ihres Computers abgelegt.

Um mit Rastern in R zu arbeiten, benötigen wir die drei Schlüsselpakete, “sp”, “raster” und “rgdal”. Laden Sie zunächst diese und weitere Pakete:

if (!require(raster)){install.packages("raster"); library(raster)}
if (!require(rgdal)){install.packages("rgdal"); library(rgdal)}
if (!require(sp)){install.packages("sp"); library(sp)}
if (!require(sf)){install.packages("sf"); library(sf)}
if (!require(rgeos)){install.packages("rgeos"); library(rgeos)}
if (!require(topmodel)){install.packages("topmodel"); library(topmodel)}

Raster laden

Raster als R-Objekt “SRTM” von Github laden und dessen Eigenschaften überprüfen. Die Daten kommen ursprünglich von der Shuttle Radar Topography Mission (SRTM) der NASA.

urlfile1 <- "https://raw.githubusercontent.com/RafHo/teaching/master/angewandte_geodatenverarbeitung/datasource/srtm_39_03_extract.tif" # get url from github repositury

srtm1 <- raster(urlfile1)

urlfile2 <- "https://raw.githubusercontent.com/RafHo/teaching/master/angewandte_geodatenverarbeitung/datasource/srtm_38_03_extract.tif" # get url from github repositury
srtm2 <- raster(urlfile2)

#Mosaic/merge srtm tiles
srtmmosaic <- mosaic(srtm1, srtm2, fun=mean)

#Attribute der Raster erkunden
srtmmosaic
## class      : RasterLayer 
## dimensions : 2422, 1924, 4659928  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : 9.08375, 10.68708, 47.37542, 49.39375  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 127, 2528  (min, max)
cellStats(srtmmosaic, min) #Statistiken auslesen
## [1] 127
cellStats(srtmmosaic, max) #Statistiken auslesen
## [1] 2528
srtmmosaic <- setMinMax(srtmmosaic) #MinMax als MetaInfo definieren

srtmmosaic
## class      : RasterLayer 
## dimensions : 2422, 1924, 4659928  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : 9.08375, 10.68708, 47.37542, 49.39375  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 127, 2528  (min, max)

Welche Eigenschaften weist das Raster noch auf

srtmmosaic@crs #Cordinate Reference System
## CRS arguments:
##  +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
srtmmosaic@extent #Ausdehnung
## class      : Extent 
## xmin       : 9.08375 
## xmax       : 10.68708 
## ymin       : 47.37542 
## ymax       : 49.39375
#the distribution of values in the raster
hist(srtmmosaic, main="Distribution of elevation values", col= "blue")

Raster zuschneiden und plotten

#load subset shapefile from Github
download.file("https://raw.githubusercontent.com/RafHo/teaching/master/angewandte_geodatenverarbeitung/datasource/subset.zip", destfile = "subset.zip")

# unzip shapefile
unzip("subset.zip")

# unzip shapefile
subset <- st_read("subset.shp")
## Reading layer `subset' from data source `V:\github\nils_website\blogdown\content\courses\angewandte_geodatenverarbeitung\subset.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 8 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: 9.506647 ymin: 48.9927 xmax: 10.07162 ymax: 49.35342
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
#Transform subset data to other CRS for cropping
subset_transformed <- st_transform(subset, "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")

extent(subset_transformed) 
## class      : Extent 
## xmin       : 9.506647 
## xmax       : 10.07162 
## ymin       : 48.9927 
## ymax       : 49.35342
#Crop raster data by extent of subset_transform
srtm_croped <- crop(srtmmosaic, extent(subset_transformed)) 
#alternativ: extent(9.506647, 10.07162, 48.9927, 49.35342 )

#hist(srtm_croped, main="Distribution of elevation values", col= "blue")
#extent(srtmmosaic)

plot(srtm_croped, main="Elevation (SRTM)")
plot(subset, add=TRUE)
## Warning in plot.sf(subset, add = TRUE): ignoring all but the first attribute

Morphometrische Analysen

#Computate slope and aspect
dem_aspect<- terrain(srtm_croped, opt = "aspect")
dem_slope<- terrain(srtm_croped, opt = "slope")
dem_flowdir<- terrain(srtm_croped, opt = "flowdir")

dem_terrain <- terrain(srtm_croped, opt = c("slope", "aspect","flowdir"))
plot(dem_terrain)

Speichern der erstellten Plots als PDF

# saving only when locally
# as pdf
pdf(paste("YOUR/FILEPATH","dem_aspect.pdf", sep=""))
plot(dem_aspect, main="Aspect", col=terrain.colors(100))
dev.off()

# # saving only when locally
pdf(paste("YOUR/FILEPATH","/dem_slope.pdf", sep=""))
plot(dem_slope, main="Slope", col=terrain.colors(100))
dev.off()

Schummerungskarte und Höhenlinien

#Hillshade
hill <- hillShade(dem_slope, dem_aspect, 40, 270)
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Study Area")
plot(srtm_croped, col = terrain.colors(100, alpha = 0.50), add = TRUE)

#Countourlines
cont <- rasterToContour(srtm_croped, levels= c(200,300,400,500,600))

#plot on top of the hillshade
plot.new()
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Study Area")
plot(srtm_croped, col = terrain.colors(100, alpha = 0.50), add = TRUE)
plot(cont, add=TRUE)

Karte in 3D erkunden

#view in 3D
if (!require(rasterVis)){install.packages("rasterVis"); library(rasterVis)}
if (!require(rgl)){install.packages("rgl"); library(rgl)}

plot3D(srtm_croped)

Mit WORLDCLIM-Daten arbeiten

WorldClim ist eine Datenbank mit globalen Wetter- und Klimadaten mit hoher räumlicher Auflösung. Diese Daten können für Kartierungen und räumliche Modellierung verwendet werden. In R lassen sich die Worldclim Daten über die Funktion ‘getData’ direkt laden. Im aufgeführten Beispiel wird aus der Datenbank ‘worldclim’ die Variable ‘bio’ in der angegebenen Auflösung für die entsprechende Länge und Breite gewählt. Die bioklimatischen Variablen wurden aus Temperatur und Niederschlag errechnet.

Daten laden

#Derive raster data by worldclim
bio <- getData('worldclim', var='bio', res=0.5, lon=5, lat=45)

Rasterdaten speichern

#write Raster into working directory
writeRaster(bio, filename=names(bio), bylayer=TRUE, format="GTiff", overwrite=TRUE)

Auswählen und zuschneiden

#Crop raster data by extent of state subset
prec_wettest_m_sub <- crop(bio$bio13_16, extent(subset)) #The precipitation of the wettest month.
prec_driest_m_sub <- crop(bio$bio14_16, extent(subset)) #The precipitation of the driest week or month.
prec_annual <- crop(bio$bio12_16, extent(subset)) #The sum of all the monthly precipitation estimates.
prec_seasonality  <- crop(bio$bio15_16, extent(subset)) #The Coefficient of Variation (C of V) is the standard deviation of the monthly precipitation estimates expressed as a percentage of the mean of those estimates (i.e. the annual mean).

Mit ASCI-Daten des DWD arbeiten

Daten laden, transformieren und zuschneiden

#Load data
urlfile3 <- "https://raw.githubusercontent.com/RafHo/teaching/master/angewandte_geodatenverarbeitung/datasource/RW_20160529-1650.asc" # get url from github repositury

dwd <- raster(urlfile3)

#Assign a CRS to raster spatial data
projection(dwd) <- CRS("+proj=stere +lat_0=90.0 +lon_0=10.0 +lat_ts=60.0 +a=6370040 +b=6370040 +units=m")

#Transform rsater spatial data to other CRS
dwd_wgs84 <- projectRaster(dwd, crs="+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")  # geographical, datum WGS84

#Crop raster data by extent of subset
croped_raster <- crop(dwd_wgs84, extent(subset))

#Mask raster data by extent of state subset
croped_raster_mask <- mask(croped_raster, subset)

#Get shp-file of Germany
germany0 <- getData('GADM', country='DE', level=0)
#Set the projection of the shp-file
proj4string(germany0) <- CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
## Warning in `proj4string<-`(`*tmp*`, value = new("CRS", projargs = "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")): A new CRS was assigned to an object with an existing CRS:
## +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## without reprojecting.
## For reprojection, use function spTransform
#Crop raster data by extent of state subset
croped_raster_ov <- mask(dwd_wgs84, germany0)

Daten umrechnen und plotten

#Calculate 1/10-value of RADOLAN-Data
c_10 <- croped_raster_mask/10

#Calculate 1/10-value of RADOLAN-Data
 c_10_ov <- croped_raster_ov/10

#Define Breaks
brk_g <- c(0,10,20,30,40,50,60,70,80,90)
brk <- c(0,10,20,30,40,50,60,70)

#Plot maps
plot(c_10, main="2016-05-29-16:50", breaks=brk, col=terrain.colors(8))

plot(c_10_ov, main="2016-05-29-16:50", breaks=brk_g, col=terrain.colors(10))

Mehrere Daten via For-Schleife einladen

For-Schleifen eignen sich für wiederkehrende Prozesse, auch Batchprozesse genannt. Mit Ihnen lassen sich Operationen automatisiert anwenden. Hier Werden aus dem working directory alle Dateien mit der Dateiendung .tif eingeladen.

#change working directory
#setwd("your/Filepath")

#Jahresmittel der monatlich gemittelten t?glichen Lufttemperatur 
#in 2 m H?hhe, in 1/10?C
#metadaten: ftp://ftp-cdc.dwd.de/pub/CDC/grids_germany/annual/air_temperature_mean/BESCHREIBUNG_gridsgermany_annual_air_temperature_mean_de.pdf

#Liste erstellen aufgrund von Pattern
rlist=list.files(getwd(), pattern="tif$", full.names=F)

# ###load all rasters
for(i in 1:length(rlist)){ 
assign(unlist(strsplit(rlist[i], "[.]"))[1], raster(rlist[i])) } 
Previous
Next