The raster package simplifies the access to large raster data sets and extends the analytic tools available for both raster and vector data. Used with the rasterVis package, it also provides enhanced visualization.

library(raster)
library(rasterVis)
## load Berlin spatial data
load(url("https://userpage.fu-berlin.de/soga/data/r-data/berlin_districts_sp_wgs84.RData"))
load(url("https://userpage.fu-berlin.de/soga/data/r-data/berlin_geocampus.RData"))

To showcase some of the features provided by the raster and rasterVis packages we download the SRTM digital elevation model (DEM) for Germany. The data was retrieved from OPENDEM on June 25, 2017. Please note 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.

## read the data
srtm <- raster("srtm_germany_dtm.tif")
## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
## Path to GDAL shared files: C:/Users/johan/AppData/Local/R/win-library/4.3/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
## Path to PROJ shared files: C:/Users/johan/AppData/Local/R/win-library/4.3/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

Printing the RasterLayer object summarizes the raster data set.

srtm
## class      : RasterLayer 
## dimensions : 10801, 13201, 142584001  (nrow, ncol, ncell)
## resolution : 0.0008333333, 0.0008333333  (x, y)
## extent     : 4.999583, 16.00042, 46.99958, 56.00042  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : srtm_germany_dtm.tif 
## names      : srtm_germany_dtm

When working with large data sets it is recommendable to check the size of the file in advance. The processing time of large data may be quite long.

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 the aggregate() function. The aggregate() function performs aggregation by grouping adjacent pixels into new ones, given an input raster and the aggregation factor fact. For example, using fact = 5 means that each set of 5 \(\times\) 5 cells in the original raster becomes a single cell in the new raster. In addition, the fun argument specifies which function will be used to calculate the new values based on existing ones (the default is mean).

For demonstration purposes we apply the aggregate() function on the SRTM raster, aggregating 5 \(\times\) 5 sets of cells.

srtm_small <- aggregate(srtm, 5, fun = mean)
plot(srtm_small)

Further, in order to clean up the data set we exclude the area that covers the ocean. We can achieve this by setting all values that are equal to 0 to NA.

srtm_small[srtm_small == 0] <- NA

Now we plot the aggregated and cleaned up data set.

library(viridisLite)

plot(srtm_small,
  main = "DEM-SRTM Germany",
  col = viridis(1000)
)

Well, still not perfect but we keep it like this for now.

Now we want to subset the raster data set to reflect our area of interest (AOI). Therefore, we apply the crop() function and feed it with a sp object, in our case the district border of the city of Berlin (berlin.sp.wgs84).

srtm_berlin <- crop(srtm, berlin.sp.wgs84)

We may plot the resulting subset using the generic plot() function. In order to add additional layers we simply call the plot() function once again, this time, however, adding the argument add = TRUE to the function call. To showcase this functionality we plot our AOI (the area of Berlin) and add the SpatialPoints object geocampus.sp, representing the location of the GeoCampus Lankwitz.

plot(srtm_berlin,
  main = "DEM-SRTM Berlin and GeoCampus Lankwitz",
  cex.main = 0.7,
  col = viridis(10)
)
plot(geocampus.sp, add = TRUE)

Another 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 use the mask() function to mask out the area outside the district of Berlin. We use the objects srtm_berlin and berlin.sp.wgs84 and plot the result.

srtm_berlin_masked <- mask(srtm_berlin, berlin.sp.wgs84)
plot(srtm_berlin_masked,
  main = "DEM-SRTM Berlin and GeoCampus Lankwitz",
  cex.main = 0.7,
  col = viridisLite::viridis(10)
)
plot(geocampus.sp, add = T)

The rasterVis package defines additional methods for spatial raster data access, manipulation and visualization (for more information on the rasterVis package go here).

Here, we apply the levelplot() function which is optimized for displaying raster data and provides some very nice additional graphic features. Type ??rasterVis::levelplot in your console for further information.

m <- 20
levelplot(srtm_berlin,
  margin = FALSE,
  contour = TRUE,
  at = seq(0, maxValue(srtm_berlin), m),
  main = paste("DEM-SRTM, contour lines", m, "m"),
  par.settings = rasterTheme()
)

Another useful function in the raster package is the terrain() function. This function computes the slope, aspect and other terrain characteristics from a raster with elevation data.

slope <- terrain(srtm_berlin,
  opt = "slope",
  unit = "degrees"
)
levelplot(slope,
  margin = FALSE,
  par.settings = BuRdTheme(),
  main = "Slope"
)

The levelpot() function returns a trellis object, which relates to the lattice graphics package. In order to overlay different trellis objects we make use of the layer() function of the latticeExtra package. For example, let us add the administrative border of Berlin. We overlay trellis objects with the + sign.

library(latticeExtra)

levelplot(slope,
  margin = FALSE,
  par.settings = BuRdTheme(),
  main = "Slope"
) + layer(sp.polygons(berlin.sp.wgs84))

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.

Creative Commons License
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.