Geomorphometric variables calculation
Ivan Lizarazo
October 19, 2016
Introduction
This document illustrates calculation of geomorphometric variables using raster library (Hijmans & van Etten, 2016). Furthermore, it shows an example of landform classification.
ESRI ASCII raster format
The ESRI ASCII raster format can be used to transfer information to or from other cell-based or raster systems. When an existing raster is output to an ESRI ASCII format raster, the file will begin with header information that defines the properties of the raster such as the cell size, the number of rows and columns, and the coordinates of the origin of the raster. The header information is followed by cell value information specified in space-delimited row-major order, with each row seperated by a carraige return. In order to convert an ASCII file to a raster, the data must be in this same format. The parameters in the header part of the file must match correctly with the structure of the data values. The basic structure of the ESRI ASCII raster has the header information at the beginning of the file followed by the cell value data:
The syntax of the header information is a keyword paired with the value of that keyword. The definitions of the kewords are:
## Warning: package 'pander' was built under R version 3.2.5| Parameter | Description | Requirements |
| NCOLS | Number of cell columns. | Integer greater than 0. |
| NROWS | Number of cell rows. | Integer greater than 0. |
| XLLCENTER or XLLCORNER | X coordinate of the origin | Match with Y coordinate type. |
| YLLCENTER or YLLCORNER | Y coordinate of the origin | Match with X coordinate type. |
| CELLSIZE | Cell size | Greater than 0. |
| NODATA_VALUE | The input NoData values | Optional. Default is -9999. |
Download data
Before starting, download DEMs provided by Prof Jantiene Baartman at Wageningen University and Research Centre from this link.
Unzip the DEMs in your working directory.
List available DEM’s
As with any project, the starting point is to load the data weâll be using and to understand its properties.
# define working directory
# setwd("./documentos/rpubs")
# load the packages we'll be using
library(rgdal)## Warning: package 'rgdal' was built under R version 3.2.5## Loading required package: sp## Warning: package 'sp' was built under R version 3.2.5## rgdal: version: 1.1-10, (SVN revision 622)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: C:/Users/IVAN/Documents/R/win-library/3.2/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: C:/Users/IVAN/Documents/R/win-library/3.2/rgdal/proj
## Linking to sp version: 1.2-3library(sp)
library(raster)## Warning: package 'raster' was built under R version 3.2.5#
# list files:
list.files(path="./DEMs/", pattern=".asc")## [1] "dem1_deathvalley2_filled0.asc" "dem10_agrit_upperhalf.asc"
## [3] "dem11_agrit_full.asc" "dem2_oregon_filled0.asc"
## [5] "dem2_v_without_fp.asc" "dem3_england_400x400_filled0.asc"
## [7] "dem3_v_with_fp.asc" "dem4_u_without_fp.asc"
## [9] "dem5_maas_filled0.asc" "dem6_1t_morph.asc"
## [11] "dem6_athabasca_filled0.asc" "dem7_2t_morph.asc"
## [13] "dem7_swiss_filled0.asc" "dem8_3t_morph.asc"
## [15] "dem9_agrit_lowerhalf.asc"Select the Death Valley DEM
Let’s load the selected DEM dem1_deathvalley2_filled0.asc
#
# read municipios data:
valley <- raster('./DEMs/dem1_deathvalley2_filled0.asc')
# what is valley?
valley## class : RasterLayer
## dimensions : 400, 400, 160000 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 0, 12000, 0, 12000 (xmin, xmax, ymin, ymax)
## coord. ref. : NA
## data source : C:\Users\IVAN\Documents\rpubs\DEMs\dem1_deathvalley2_filled0.asc
## names : dem1_deathvalley2_filled0
## values : -2147483648, 2147483647 (min, max)# define projection
proj4string(valley) <- "+proj=tmerc +lat_0=0 +lon_0=0 +k=9996 +x_0=500000 +y_0=10000000"
# what has change in the DEM?
valley## class : RasterLayer
## dimensions : 400, 400, 160000 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 0, 12000, 0, 12000 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=0 +lon_0=0 +k=9996 +x_0=500000 +y_0=10000000 +ellps=WGS84
## data source : C:\Users\IVAN\Documents\rpubs\DEMs\dem1_deathvalley2_filled0.asc
## names : dem1_deathvalley2_filled0
## values : -2147483648, 2147483647 (min, max)# plot DEM
plot(valley)
# add DEM values
#text(valley, col="darkred",
#cex=0.6, font=2)
# add a plot title and legend
title("Death Valley DEM")
#The terrain function
The terrain function compute slope, aspect and other terrain characteristics from a raster with elevation data. The elevation data should be in map units (typically meter) for projected (planar) raster data. They should be in meters when the coordinate reference system (CRS) is longitude/latitude.
Usage
terrain(x, opt=‘slope’, unit=‘radians’, neighbors=8, filename=’’, …)
Arguments
x is a RasterLayer object with elevation values. Values should have the same unit as the map units, or in meters when the crs is longitude/latitude
opt is a character vector containing one or more of these options: slope, aspect, TPI, TRI, roughness, flowdir
unit is a Character. ‘degrees’, ‘radians’ or ‘tangent’. Only relevant for slope and aspect. If ‘tangent’ is selected that is used for slope, but for aspect ‘degrees’ is used (as ‘tangent’ has no meaning for aspect)
neighbors is an integer indicating how many neighboring cells to use to compute slope for any cell. Either 8 (queen case) or 4 (rook case). Only used for slope and aspect
filename is a character indicating output filename (optional)
Details
When neighbors=4, slope and aspect are computed according to Fleming and Hoffer (1979) and Ritter (1987). When neigbors=8, slope and aspect are computed according to Horn (1981). The Horn algorithm may be best for rough surfaces, and the Fleming and Hoffer algorithm may be better for smoother surfaces (Jones, 1997; Burrough and McDonnell, 1998). If slope = 0, aspect is set to 0.5*pi radians (or 90 degrees if unit=‘degrees’). When computing slope or aspect, the CRS (projection) of the RasterLayer x must be known (may not be NA), to be able to safely differentiate between planar and longitude/latitude data.
flowdir returns the ‘flow direction’ (of water), i.e. the direction of the greatest drop in elevation (or the smallest rise if all neighbors are higher). They are encoded as powers of 2 (0 to 7). The cell to the right of the focal cell ‘x’ is 1, the one below that is 2, and so on:
32 64 128 16 x 1 8 4 2
If two cells have the same drop in elevation, a random cell is picked. That is not ideal as it may prevent the creation of connected flow networks. ArcGIS implements the approach of Greenlee (1987) and I might adopt that in the future.
The terrain indices are according to Wilson et al. (2007), as in gdaldem.
TRI (Terrain Ruggedness Index) is the mean of the absolute differences between the value of a cell and the value of its 8 surrounding cells.
TPI (Topographic Position Index) is the difference between the value of a cell and the mean value of its 8 surrounding cells.
Roughness is the difference between the maximum and the minimum value of a cell and its 8 surrounding cells.
Such measures can also be computed with the focal function:
f <- matrix(1, nrow=3, ncol=3)
TRI <- focal(x, w=f, fun=function(x, …) sum(abs(x[-5]-x[5]))/8, pad=TRUE, padValue=NA)
TPI <- focal(x, w=f, fun=function(x, …) x[5] - mean(x[-5]), pad=TRUE, padValue=NA)
rough <- focal(x, w=f, fun=function(x, …) max(x) - min(x), pad=TRUE, padValue=NA, na.rm=TRUE)
Slope and Aspect
Let’s calculate slope and aspect:
x <- terrain(valley, opt = c("slope", "aspect"), unit = "degrees")
plot(x)
Hill shading
The hillShade function compute hill shade from slope and aspect layers (both in radians). A hill shade layer is often used as a backdrop on top of which another, semi-transparent, layer is drawn.
pendiente <- terrain(valley, opt = "slope")
aspecto <- terrain(valley, opt = "aspect")
hill <- hillShade(pendiente, aspecto, 40, 270)
plot(hill, col = grey(0:100/100), legend = FALSE, main = "Death Valley")
plot(valley, col = rainbow(25, alpha = 0.35), add = TRUE)
Topographic Indices
Let’s calculate TPI using a 5x5 vindow:
# TPI for different neighborhood size:
# first step: define customized function
tpiw <- function(x, w=5) {
m <- matrix(1/(w^2-1), nc=w, nr=w)
m[ceiling(0.5 * length(m))] <- 0
f <- focal(x, m)
x - f
}
# second step: apply the function
tpi5 <- tpiw(valley, w=5)
tpi5## class : RasterLayer
## dimensions : 400, 400, 160000 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 0, 12000, 0, 12000 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=0 +lon_0=0 +k=9996 +x_0=500000 +y_0=10000000 +ellps=WGS84
## data source : in memory
## names : layer
## values : -31.91667, 34.33333 (min, max)tpi5n <- setMinMax(tpi5)
col <- rainbow(20)
plot(tpi5n, col=col, main="Topographic Position Index in Death Valley")
Have you understood how the tpiw function works? If no, it’s time now to review it. If yes, try writing a function to compute TRI and TPI using a 3x3 window. Compare your results with those of the raster library.
A 3D visualization
The plot3D() function in the RasterVis package can be useful some times:
library(rasterVis)## Loading required package: lattice## Loading required package: latticeExtra## Loading required package: RColorBrewerplot3D(valley) 
Another plot of slope
The vectorplot function in the RasterVis package can be useful for displaying the slope vector field as a set of arrows:
vectorplot(x, par.settings=RdBuTheme())
Contour
The levelplot function in the RasterVis package can be useful for displaying a DEM in several ways:
library(colorspace)##
## Attaching package: 'colorspace'## The following object is masked from 'package:raster':
##
## RGBmyTheme <- rasterTheme(region=sequential_hcl(10, power=2.2))
levelplot(valley, par.settings = myTheme, contour = TRUE)
levelplot(valley, par.settings = RdBuTheme, contour = TRUE)
levels <- levels <- seq(100, 1500, 100)
contourplot(valley, at = levels, margin = FALSE)## Warning in .local(x, data, ...): FUN.margin is deprecated. Use margin as a
## list instead.## Warning in .local(x, data, ...): axis.margin is deprecated. Use margin as a
## list instead.
Landform classification
Let’s examine TPI and use it to classify landforms using Robby R. Marrotte’s code:
# Some stats
quantile(tpi5[],na.rm=T)## 0% 25% 50% 75% 100%
## -31.91666667 -0.75000000 -0.04166667 0.58333333 34.33333333hist(tpi5[])
# Get the standard deviation of the TPI
SD <- sd(tpi5[],na.rm=T)
# Make landform classes
#Morphologic class De Reu et al. 2013; Weiss (2001)
landform <- reclassify(tpi5, matrix(c(-Inf, -SD, 1,
-SD, -SD/2, 2,
-SD/2, 0, 3,
0, SD/2, 4,
SD/2, SD, 5,
SD, Inf, 6),
ncol = 3, byrow = T),
right = T)
# Turn it into categorical raster
landform <- as.factor(landform)
rat <- levels(landform)[[1]]
rat[["landform"]] <- c('Valley', 'Lower Slope',
'Flat Area','Middle Slope',
'Upper Slope', 'Ridge')
levels(landform) <- rat
# Plot the classification
x11(12,12)
levelplot(landform, col.regions = rev(brewer.pal(6,'RdYlBu')),
labels = rat$landcover,
main = "Landform Classification",
colorkey=list(labels=list(at=1:6, labels=rat[["landform"]])))
require(rasterVis)
require(dismo)## Loading required package: dismo## Warning: package 'dismo' was built under R version 3.2.5require(GISTools) ## Loading required package: GISTools## Warning: package 'GISTools' was built under R version 3.2.5## Loading required package: maptools## Warning: package 'maptools' was built under R version 3.2.5## Checking rgeos availability: TRUE## Loading required package: MASS##
## Attaching package: 'MASS'## The following objects are masked from 'package:raster':
##
## area, select## Loading required package: rgeos## Warning: package 'rgeos' was built under R version 3.2.5## rgeos version: 0.3-20, (SVN revision 535)
## GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084
## Linking to sp version: 1.2-3
## Polygon checking: TRUErequire(spdep)## Loading required package: spdep## Warning: package 'spdep' was built under R version 3.2.5## Loading required package: Matrixrequire(RColorBrewer)
# Extract ridges from the classification
ridges <- landform
ridges[] <- ifelse(ridges[] != 6, NA, 6)
# Reset the levels
ridges <- ratify(ridges)
# Plot it
x11(16,12)
plot(valley)
plot(ridges, col = adjustcolor("red", alpha.f = 0.4) , add = T, legend = F, main="Ridges in Death Valley")
Now, it’s your turn to code
Do not forget to replicate this exercise. Furthermore, compute geomorphometric variables for the remaining DEMs in the folder you downloaded.
Remember: Do not copy & paste. If you are serious about learning R, get your hands dirty and write your code.
Need more?
Try calculation of geomorphometric variables using your favorite GIS software (i.e. ArcGIS, QGIS, SAGA GIS).
REFERENCES
Burrough, P., and R.A. McDonnell, 1998. Principles of Geographical Information Systems. Oxford University Press.
Fleming, M.D. and Hoffer, R.M., 1979. Machine processing of landsat MSS data and DMA topographic data for forest cover type mapping. LARS Technical Report 062879. Laboratory for Applications of Remote Sensing, Purdue University, West Lafayette, Indiana.
Greenlee, D.D., 1987. Raster and vector processing for scanned linework. Photogrammetric Engineering and Remote Sensing 53:1383-1387
Hijmans, R.J & van Etten, J., 2016. raster: Geographic analysis and modeling with raster data. R package version 2.5-16. http://CRAN.R-project.org/package=raster
Horn, B.K.P., 1981. Hill shading and the reflectance map. Proceedings of the IEEE 69:14-47
Jones, K.H., 1998. A comparison of algorithms used to compute hill slope as a property of the DEM. Computers & Geosciences 24: 315-323
Ritter, P., 1987. A vector-based slope and aspect generation algorithm. Photogrammetric Engineering and Remote Sensing 53: 1109-1111
Wilson, M.F.J., O’Connell, B., Brown, C., Guinan, J.C., Grehan, A.J., 2007. Multiscale terrain analysis of multibeam bathymetry data for habitat mapping on the continental slope. Marine Geodesy 30: 3-35.