11140_the_raster_and_rastervis

The terra package provides bindings to the Geospatial Data Abstraction Library (GDAL) and access to projection and transformation operations from the PROJ library. This software is the basis for spatial software applications and geodata analysis within R and is the successor of the raster package. Its bindings allow fast I/O (input/output) operations for GDAL-supported raster formats, OGR-supported vector formats and numerous utility functions to convert geographic longitude and latitude coordinates into Cartesian coordinates (and vice versa). Furthermore, the terra package simplifies the access to large raster data sets and extends the analytic tools available for both raster and vector data. Additionally, the rasterVis package provides enhanced visualization opportunities.

library(terra)
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"))

If you have not installed the required packages yet, you can do so by install.packages("terra", "rasterVIS")

To showcase some of the features provided by the terra 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 rast() function out of the terra package.

## read the data
srtm <- rast("srtm_germany_dtm.tif")

Printing the SpatRaster object summarizes the raster data set.

srtm

## class       : SpatRaster 
## dimensions  : 10801, 13201, 1  (nrow, ncol, nlyr)
## resolution  : 0.0008333333, 0.0008333333  (x, y)
## extent      : 4.999583, 16.00042, 46.99958, 56.00042  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source      : srtm_germany_dtm.tif 
## name        : srtm_germany_dtm

When working with large data sets, checking the file size in advance is recommended. 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 downscaled 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, main = "DEM-SRTM with default colormap")

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

srtm_small[srtm_small == 0] <- NA

Now, we plot the aggregated and cleaned-up data set. Therefore, we want to create two plots that show the same visualization with different implementations. For the first one, we only use functions provided directly by the base R and the terra package. We use the additional library viridisLite for the second plot, which provides several colorblind-friendly colour maps for visualization purposes.

par(mfrow = c(1,2))

# creating an own colorRamp for DEM-Visualization
mycols <- colorRampPalette(c("#4a235a","slateblue4", "lightseagreen", "seagreen3",
                             "chartreuse2", "yellowgreen", "lightgoldenrod1"))
plot(srtm_small,
     main = "DEM-SRTM Germany with own colorramp",
     col=mycols(100), cex = 1.2)

# Create the same plot by using the viridisLite library
library(viridisLite)
plot(srtm_small,
  main = "DEM-SRTM Germany",
  col = viridis(1000), cex = 1.2)

Well, still there could be done some further adaptations, but we keep it 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 call the plot() function once again; this time, however, we add 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, which represents 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 with the same values as the input Raster object, except for the masked cells, which are set to NA by default. Therefore, we need to transform the sp object to the equivalent terra package representation by using the vect() function.

We use the mask() function to mask the area outside the district of Berlin. We plot the result using the objects srtm_berlin and geocampus.sp.

berlin_vect <- vect(berlin.sp.wgs84)

srtm_berlin_masked <- mask(srtm_berlin, berlin_vect)
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, minmax(srtm_berlin)[2], m),
          main = paste("DEM-SRTM, contour lines", m, "m"),
          par.settings = rasterTheme()
)

Another useful function in the terra 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,
                 v = "slope", neighbors=8,
                 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)
library(sp)

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.

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.