30914_SRTM_Digital_Elevation_Model_for

The SRTM digital elevation model (DEM) for Germany was downloaded 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"

Due 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. , meaning that each new pixel represents the average of a 20×20 block of original pixels. 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 in this context is mask(), which creates a new Raster object by copying the values of the input raster but setting all cells outside a specified area to NA by default. In our case, we use the mask() function to remove areas outside the borders of Germany. We can easily retrieve the national borders of Germany by the gadm() function from the geodata package. Again, depending on your computing resources this step may take some time.

# Retrieve Federal States by the gadm() function from the geodata package
library(geodata)

germany_gadm <- gadm("Germany", level = 1, path = tempdir())
germany_sp <- as(germany_gadm, "Spatial") # Convert SpatVector to Spatial object for compatibility with raster::mask()
srtm_germany_masked <- mask(srtm_germany, germany_sp)

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
)

To save the projected raster for future use, 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 (defining the area of interest), in our case the borders of East Germany. We can easily retrieve the Federal States borders of Germany by the gadm() function from the geodata package.

library(sf)
library(tidyverse)
# Retrieve Federal States using the gadm() function from the geodata package
germany <- gadm("Germany", level = 1, path = tempdir())

# 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_1) %>%
  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 uses geographic coordinates by default. For further analysis, we re-project the data set into the UTM, Zone 33N coordinate reference system providing the EPSG identifier \(32633\).

To ensure compatibility, we also reproject all other relevant data sets to use 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.