This is an introductory tutorial on basic raster analysis carried out for a study site over Attica, Greece. The main tool that will be used is the raster library [1] in R. The objectives of the current tutorial are to:

  1. Get familiar with the handling of a DEM file
  2. Derive slope and aspect surfaces from a DEM
  3. Plot raster layers

The DEM file we will work with is the GTOPO30 DEM and was derived from here [2]. The original .tiff was masked to the study area (Attica).

So, let’s start!

library("maps")
library("rgdal")
## Loading required package: sp
## rgdal: version: 1.2-16, (SVN revision 701)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.0, released 2017/04/28
##  Path to GDAL shared files: F:/R-3.4.3/library/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: F:/R-3.4.3/library/rgdal/proj
##  Linking to sp version: 1.2-7
library('sf')
## Warning: package 'sf' was built under R version 3.4.4
## Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
library('sp')
library("raster")
library("mapdata")
library("colorRamps")
library("maptools")
## Checking rgeos availability: TRUE
library("rgeos")
## rgeos version: 0.3-26, (SVN revision 560)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 r0 
##  Linking to sp version: 1.2-7 
##  Polygon checking: TRUE
library('mapview')
## Warning: package 'mapview' was built under R version 3.4.4
#
# Here you have to set your own path to your working directory
setwd("C:/Users/user/Documents/Tutorials/Tutorial1/")
getwd()
## [1] "C:/Users/user/Documents/Tutorials/Tutorial1"
#
dem_attica = raster("dem_attica.tif")
# Here you can view some basic attributes of the file we work with, concerning its resolution, its spatial extent, its coordinate system and the number of rows and columns that constitute the raster file.
dem_attica 
## class       : RasterLayer 
## dimensions  : 106, 121, 12826  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 23.10833, 24.11667, 37.45833, 38.34167  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : C:\Users\user\Documents\Tutorials\Tutorial1\dem_attica.tif 
## names       : dem_attica 
## values      : 1, 1285  (min, max)
summary(dem_attica)
##         dem_attica
## Min.             1
## 1st Qu.        107
## Median         228
## 3rd Qu.        403
## Max.          1285
## NA's          8634
#
# We retreive the boundary shapefile of all the local municipalities over Greece, corresponding to level=3. Here [3] you can see in more detail the various administrative levels that exist and how they differ.
greece = getData('GADM', country='GR', level=3)
# This is how our shapefile looks like
plot(greece)

#
# Plotting the DEM file
brks = seq(0,1300,100)
plot(dem_attica, breaks=brks, col=terrain.colors(13), legend.shrink=1.1, legend.width=1.8,
main = "Elevation map of Attica", legend.args=list(text='Elevation (m)', side=4, font=2, line=2.5, cex=0.8))
plot(greece, add=T)

#
# Adding contours in the elevation map. 
levels <- seq(0,1300,300)
plot(dem_attica)
contour(dem_attica, levels=levels, add=T)

#
# Let's see how the values are distributed though our domain
hist(dem_attica, las=1, col="grey", breaks="Scott", main="Histogram of DEM values over Attica per pixel", 
     xlab="Elevation")

As it can be seen from the histogram above, the majority of pixels are found over low elevation areas. Thy y-axis of the above plot represents the number of pixels that have a specific elevation value. For instance, approximately 430 pixels have elevation values that range from 0-50 meters. Now, let’s suppose we would like to reclassify our original raster to certain elevation classes that are of interest. This could be done in the following way:

# Define the reclassification rules
m <- c(0, 200, 1,
       200, 400, 2,
       400, 600, 3,
       600, 800, 4,
       800, 1000, 5,
       1000, 1200, 6,
       1200, 1300, 7)
#
rclmat = matrix(m, ncol=3, byrow=TRUE)
dem_attica_rc = reclassify(dem_attica, rclmat)
#
labs = c("0-200", 
       "200-400", 
       "400-600", 
       "600-800", 
       "800-1000", 
       "1000-1200", 
       "1200-1300")
plot(dem_attica_rc, legend=F, legend.shrink=1.1, legend.width=1.8,
     main="Reclassified elevation map", col=rev(terrain.colors(6)))
legend("bottomright", legend=labs, fill=rev(terrain.colors(6)), x.intersp=0.3, cex=1, bg="white")
plot(greece, add=T)

#

In the map above we no longer see the original elevation values, but we see the elevation classes we have defined.

Now, let’s return to our original DEM raster and have a 3D perspective of the elevation over our study region:

library(knitr)
## Warning: package 'knitr' was built under R version 3.4.4
library(rgl)
## Warning: package 'rgl' was built under R version 3.4.4
## 
## Attaching package: 'rgl'
## The following object is masked from 'package:rgeos':
## 
##     triangulate
library(rglwidget)
## Warning: package 'rglwidget' was built under R version 3.4.4
## The functions in the rglwidget package have been moved to rgl.
library(htmltools)
## Warning: package 'htmltools' was built under R version 3.4.4
library(rasterVis)
## Warning: package 'rasterVis' was built under R version 3.4.4
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
#
knit_hooks$set(webgl = hook_webgl)
#
persp(dem_attica, exp=0.2,phi=35, xlab='Longitude', ylab='Latitude', zlab='Elevation', col="green4")

Let’s fly through the DEM! :) Pan with your mouse over the plot, zoom-in, zoom-out.

When looking over the elevation map of a region, we are not only concerned about the absolute values of elevation, but also we care about how abruptly does the elevation map changes in all direction. This quantity is termed “slope” and it is regarded as the degree of change in elevation (first derivative of the elevation map). The terrain [4] algorithm applied in R exploits Horne’s (1981) algorithm [5],[6]. For the calculation of slope and aspect, the 8 neighboring pixels of each pixel were used.

slope = terrain(dem_attica, neighbors=8, opt='slope', unit='degrees')
slope
## class       : RasterLayer 
## dimensions  : 106, 121, 12826  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 23.10833, 24.11667, 37.45833, 38.34167  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : slope 
## values      : 0.03676816, 23.09859  (min, max)
#
brks = seq(0,24,2)
plot(slope, breaks=brks, col=matlab.like(12), legend.shrink=1.1, legend.width=1.8,
     main = "Slope map of Attica", legend.args=list(text='Slope (degrees)', side=4, font=2, line=2.5, cex=0.8))
plot(greece, add=T)

The orientation of the maximum rate of change in elevation is termed as the “aspect” [7].

aspect = terrain(dem_attica, neighbors=8, opt='aspect', unit='degrees')
aspect
## class       : RasterLayer 
## dimensions  : 106, 121, 12826  (nrow, ncol, ncell)
## resolution  : 0.008333333, 0.008333333  (x, y)
## extent      : 23.10833, 24.11667, 37.45833, 38.34167  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : aspect 
## values      : 0, 359.6404  (min, max)
#
brks = seq(0,360,20)
plot(aspect, breaks=brks, col=matlab.like(19), legend.shrink=1.1, legend.width=1.8,
     main = "Aspect map of Attica", legend.args=list(text='Aspect (degrees)', side=4, font=2, line=2.5, cex=0.8))
plot(greece, add=T)

Finally, let’s create an interactive panel of maps with elevation, slope, aspect and the reclassified elevation map.

m1 <- mapview(dem_attica, legend=TRUE, map.types="Esri.WorldImagery")
m2 <- mapview(slope, legend=TRUE, map.types="Esri.WorldImagery")
m3 <- mapview(aspect, legend=TRUE, map.types="Esri.WorldImagery")
m4 <- mapview(dem_attica_rc, legend=TRUE, map.types="Esri.WorldImagery")
sync(m1, m2, m3, m4) # 3 panels synchronised
print("The End")
## [1] "The End"

I hope you enjoyed it! Please send me your feedback at karypidou@geo.auth.gr

[1] https://cran.r-project.org/web/packages/raster/raster.pdf [2] https://earthexplorer.usgs.gov/ [3] http://data-analytics.net/wp-content/uploads/2014/09/geo2.html [4] https://www.rdocumentation.org/packages/raster/versions/2.6-7/topics/terrain [5] https://ieeexplore.ieee.org/abstract/document/1456186/ [6] http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-slope-works.htm [7] https://www.rdocumentation.org/packages/raster/versions/2.6-7/topics/terrain