30914_SRTM_Digital_Elevation_Model_for

The SRTM digital elevation model (DEM) for Germany was retrieved from OPENDEM on November 12, 2022. Please note that the zipped file is about 86.4 MB in size.

## build a temporary folder on disk
temp <- tempfile()
download_url <- "https://userpage.fu-berlin.de/soga/data/raw-data/spatial/"
zipfile <- "srtm_germany_dtm.zip"
## download the file
download.file(paste0(download_url, zipfile), temp, mode = "wb")
## unzip the file
unzip(temp)
unlink(temp)

To load the raster data set we apply the raster() function.

library(raster)
## read in the data
srtm_germany <- raster("srtm_germany_dtm.tif")

When working with large data sets it is recommendable to check the size of the file in advance. The processing time increases with the size of the data set.

fs <- file.size("srtm_germany_dtm.tif")
print(paste("The file size is", as.numeric(fs) / 1e6, "MB"))

## [1] "The file size is 285.254776 MB"

Owing to the large file size we downscale the raster data set using he aggregate() function. The aggregate() function performs aggregation by grouping adjacent pixels into new ones, given an input raster and the aggregation factor fact. In addition, the aggregate() function accepts the fun argument, which specifies the function that will be used to calculate the new values based on existing ones (the default is mean). Here we set fact = 20. Depending on your computing resources this step may take some time.

srtm_germany <- aggregate(srtm_germany,
  fact = 20,
  fun = mean
)
plot(srtm_germany)

A useful function is the mask() function, which creates a new Raster object that has the same values as the input Raster object, except for the masked cells, which are set to NA by default. We apply the mask() function to mask out the area outside the borders of Germany. We can easily retrieve the national borders of Germany by the getData() function from the raster package. Again, depending on your computing resources this step may take some time.

# Retrieve Federal States by the getData() function from the raster package
germany <- getData(country = "Germany", level = 1)
srtm_germany_masked <- mask(srtm_germany, germany)

plot(srtm_germany_masked,
  main = "DEM-SRTM Germany",
  cex.main = 0.85
)

Please note that the SRTM data set is projected in geographic coordinates. For further usage we re-project the data set into the ETRS89/LAEA coordinate reference system (European Terrestrial Reference System 1989/Lambert Azimuthal Equal-Area projection coordinate reference system) providing the EPSG identifier \(3035\).

newproj <- CRS("+init=epsg:3035")
srtm_germany_masked_ETRS89 <- projectRaster(srtm_germany_masked, crs = newproj)

plot(srtm_germany_masked_ETRS89,
  main = "DEM-SRTM Germany - ETRS89",
  cex.main = 0.85
)

For later usage we write the projected raster object srtm_germany_masked_ETRS89 to disk.

writeRaster(srtm_germany_masked_ETRS89,
  filename = "srtm_germany_ETRS89.tif",
  format = "GTiff", overwrite = TRUE
)

SRTM Digital Elevation Model for East Germany

In this section we subset the SRTM raster data set to reflect a specific area of interest, East Germany. Therefore, we apply the crop() function and feed it with an sp object, in our case the borders of East Germany. We can easily retrieve the Federal States borders of Germany by the getData() function from the raster package.

library(sf)
library(tidyverse)
# Retrieve Federal States by the getData() function from the raster package
germany <- getData(country = "Germany", level = 1)

# transform to sf object for convenience
germany_sf <- st_as_sf(germany)
states <- c(
  "Sachsen", "Sachsen-Anhalt",
  "Mecklenburg-Vorpommern", "Brandenburg",
  "Thüringen", "Berlin"
)
east_germany_states <- germany_sf[germany_sf$NAME_1 %in% states, ]
east_germany <- east_germany_states %>%
  group_by(NAME_0) %>%
  summarize()

# back-transform to Spatial object to plug into the crop function
east_germany_sp <- as(east_germany, "Spatial")
# crop raster
srtm_east_germany <- raster::crop(srtm_germany, east_germany_sp)

# plot
plot(srtm_east_germany,
  main = "DEM-SRTM East Germany",
  cex.main = 0.85
)

Please note that SRTM data set is projected in geographic coordinates. For further usage we re-project the data set into the UTM, Zone 33N coordinate reference system providing the EPSG identifier \(32633\).

Further, we make sure that all other data sets of interest share the same coordinate reference system.

# specify CRS UTM 33N
newproj <- CRS("+init=epsg:32633")

# reproject raster object
srtm_east_germany_utm <- projectRaster(srtm_east_germany,
  crs = newproj
)

# reproject spatial object
east_germany_sp <- spTransform(east_germany_sp,
  CRSobj = newproj
)

# plot
plot(srtm_east_germany_utm,
  main = "DEM-SRTM East Germany - UTM 33N",
  cex.main = 0.85
)
plot(east_germany_sp, add = TRUE)

Citation

The E-Learning project SOGA-R was developed at the Department of Earth Sciences by Kai Hartmann, Joachim Krois and Annette Rudolph. You can reach us via mail by soga[at]zedat.fu-berlin.de.

You may use this project freely under the Creative Commons Attribution-ShareAlike 4.0 International License.

Please cite as follow: Hartmann, K., Krois, J., Rudolph, A. (2023): Statistics and Geodata Analysis using R (SOGA-R). Department of Earth Sciences, Freie Universitaet Berlin.