9. Raster data in R

Foreword on working directory, data and packages

The Data for this tutorial are provided via Github. In order to reproudce the code with your own data, replace the url with your local filepath

If you work locally, you may define the working directory

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

This should be the place where your data and other folders are stored. more Information can be found here

Code chunks for saving data and results are disabeled in this tutorial (using a #). Please replace (“YOUR/FILEPATH”) with your own path.

The code is written in such a way, that packages will be automatically installed (if not already installed) and loaded. Packages will be stored in the default location of your computer.

Intro

Um mit Rastern in R zu arbeiten, benötigen wir die drei Schlüsselpakete, “sp”, “raster” und “rgdal”:

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 `C:\Users\Nils Riach\Github\spatial_thinking_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