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 
## 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 
## source     : memory
## names      : layer 
## values     : 127, 2528  (min, max)

Welche Eigenschaften weist das Raster noch auf

srtmmosaic@crs #Cordinate Reference System

## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

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
## Bounding box:  xmin: 9.506647 ymin: 48.9927 xmax: 10.07162 ymax: 49.35342
## Geodetic CRS:  WGS 84

#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)

Daten laden

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

## Warning in getData("worldclim", var = "bio", res = 0.5, lon = 5, lat = 45): getData will be removed in a future version of raster
## . You can use access these data, and more, with functions from the geodata package instead

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).

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")

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps unknown in Proj4 definition: +proj=stere +lat_0=90
## +lat_ts=60 +lon_0=10 +x_0=0 +y_0=0 +R=6370040 +units=m +no_defs +type=crs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum unknown in Proj4 definition

#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(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html

## Warning in `proj4string<-`(`*tmp*`, value = new("CRS", projargs = "+proj=longlat +ellps=WGS84 +no_defs")): A new CRS was assigned to an object with an existing CRS:
## +proj=longlat +datum=WGS84 +no_defs
## 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))