# Landform Classification in R

Jasiewicz and Stepinkski (2013) provided a novel method for classification and mapping of landform elements from a digital elevation model (DEM) based on the principle of pattern recognition rather than differential geometry. At the core of the method is the concept of a geomorphon (geomorphologic phonotypes) — a simple ternary pattern that serves as an archetype of a particular terrain morphology. From a total of 498 possible combinations a set of ten distinct patterns is formed (see figure 1). A single scan of a DEM assigns an appropriate geomorphon to every cell in the raster using a procedure that self-adapts to identify the most suitable spatial scale at each location. As a result, the method classiﬁes landform elements at a range of different spatial scales with unprecedented computational efﬁciency.

The basis for this case study is the DEM model derived from NASA’s Shuttle Radar Topography Mission (SRTM) flown back in the year 2000. The original data set is scaled-up to an uniform grid size of 0.001 °. R package raster provides a focal function which is applied to the DEM data set. Both, the raster resolution and the size of the focal window determine the magnitude order of landforms which are to be detected. A window of 7 x 7 grid cells is used in this study. It allows for the detection of small-scale landforms in the complex topography in Switzerland.

### Source

Jasiewicz, Jarosław & F Stepinski, Tomasz. (2012). Geomorphons-a pattern recognition approach to classification and mapping of landforms. Geomorphology. 182. 10.1016/j.geomorph.2012.11.005. http://dx.doi.org/10.1016/j.geomorph.2012.11.005

### Example with R Dataset “Volcano”

Following example used the built-in R terrain data set *Volcano*. The script is available for download from Github https://github.com/m-saenger/geomorphon. Package *devtools* allows to source the script directly (see below).

```
library(devtools)
library(raster)
library(tidyverse)
library(reshape2)
# Source function from Github
file.name <- "https://raw.githubusercontent.com/m-saenger/geomorphon/master/geomorphon.R"
devtools::source_url(file.name)
# Definition
focal.window.size <- 7
flatness.thresh <- .5
counter <- 1
# Data
data("volcano")
dat <- volcano/1e2 # Scale z axis to get meaningful slope values
r <- raster(t(dat), xmn=0, xmx=nrow(dat), ymn=0, ymx=ncol(dat))
# Focal function
focal.function <- function(x){
geomorph(x, flatness.thresh, res(r), ncell(r), geomorph.lut.num, verbose = FALSE)
}
# Focal matrix
focal.matrix <- matrix(1, nrow = focal.window.size, ncol = focal.window.size)
# Apply focal function
r.volcano.lf <- suppressWarnings(raster::focal(r, fun = focal.function, w = focal.matrix, pad = T, padValue = NA))
# Melt raster to data frame
df.lf <- r.volcano.lf %>%
flip(direction = 2) %>%
as.data.frame(xy = T)
pl1 <- ggplot(df.lf, aes(x, y, fill = factor(layer, geomorph.def$num_lf))) +
geom_raster() +
scale_fill_manual("Landform", values = geomorph.def$colour, labels = geomorph.def$name, drop = F) +
coord_cartesian(expand = F) +
theme(plot.background = element_rect(colour = "black", size = .3)) +
labs(x = "x", y = "y", title = "Volcano Landform Classes")
pl2 <- ggplot(dat %>% reshape2::melt(), aes(Var1, Var2, fill = value)) +
geom_raster() +
scale_fill_distiller("Elevation", palette = "Spectral") +
coord_cartesian(expand = F) +
theme(plot.background = element_rect(colour = "black", size = .3)) +
labs(x = "x", y = "y", title = "Volcano Elevation")
```

### Appendix

#### Look-up tables

```
geomorph.def <- data.frame(
num_lf = 1:10,
id_lf = c("PK", "RI", "SH", "SP", "SL", "FS", "FL", "HL", "VL", "PT"),
name_en = c("Peak", "Ridge", "Shoulder", "Spur", "Slope", "Footslope", "Flat", "Hollow", "Valley", "Pit"),
colour = c("magenta", "red", "orange", "yellow", "grey40", "grey70", "grey90", "skyblue1", "dodgerblue", "royalblue3"),
stringsAsFactors = F
)
geomorph.lut <- data.frame(
V0 = c("FL", "FL", "FL", "SH", "SH", "RI", "RI", "RI", "PK"),
V1 = c("FL", "FL", "SH", "SH", "SH", "RI", "RI", "RI", NA),
V2 = c("FL", "FS", "SL", "SL", "SP", "SP", "RI", NA, NA),
V3 = c("FS", "FS", "SL", "SL", "SL", "SP", NA, NA, NA),
V4 = c("FS", "FS", "HL", "SL", "SL", NA, NA, NA, NA),
V5 = c("VL", "VL", "HL", "HL", NA, NA, NA, NA, NA),
V6 = c("VL", "VL", "VL", NA, NA, NA, NA, NA, NA),
V7 = c("VL", "VL", NA, NA, NA, NA, NA, NA, NA),
V8 = c("PT", NA, NA, NA, NA, NA, NA, NA, NA)
)
geomorph.lut <- as.matrix(geomorph.lut)
geomorph.lut.num <- matrix(match(geomorph.lut, geomorph.def$id_lf), nrow = nrow(geomorph.lut))
```

#### Function geomorph()

```
geomorph <- function(x, flatness.thresh = NA, res = NA, ncell = NA, geomorph.lut.num, verbose = F){
#' @description Note, that no performance optimisation has been done to this function, yet.
#' @author M. Sänger 2018
#' @source Jasiewicz, Stepinkski 2013
#' @param x vector from raster::focal function
#' @param flatness.thresh Flatness threshold in degrees
#' @param res resolution, same unit as values in r
#' @param ncell total number of cells in raster (used for progress bar only)
#' @param geomorph.lut.num look-up table to derive landform class from ternary patterns
# Breaks for flatness threshold
brks <- c(-Inf, -flatness.thresh, flatness.thresh, Inf)
brks.ind <- c(-1, 0, 1)
# Progress bar (if verbose = TRUE)
if(verbose){
counter <<- counter + 1
if(counter %in% round(seq(0, ncell, length.out = 20)))
cat("=", round(counter/ncell*100)); if(counter == ncell) cat("\n")
}
# Create matrix from incoming vector x
size = sqrt(length(x))
m <- matrix(x, nrow = size)
# Distance from central point to edge (number of cells)
mid <- ceiling(size/2)
# Matrix of all vectors from the central point to the octants
oct <- rbind(
ne = cbind(mid:size, mid:size),
e = cbind(mid:size, mid),
se = cbind(mid:size, mid:1),
s = cbind(mid, mid:1),
sw = cbind(mid:1, mid:1),
w = cbind(mid:1, mid),
nw = cbind(mid:1, mid:size),
n = cbind(mid, mid:size)
)
# Coordinates and cell distance (sqrt(2) for diagonals)
oct.vector <- m[oct]
cell.scaling <- rep(c(sqrt(2), 1), 4) # Horizontal cell distance in all 8 directions
cell.size <- res * cell.scaling
# Matrix octants vs. cell values
m1 <- matrix(oct.vector, nrow = 8, byrow = T)
# z diff from central point
m.diff <- m1[, -1] - m1[, 1]
# Calculate slope angle and transform to degrees
m.slope <- atan(m.diff/(cell.size * 1:ncol(m.diff)))
m.angle <- m.slope * 180/pi
# Calculate zenith and nadir angles for each octant
nadir <- 90 + apply(m.angle, 1, min, na.rm = T)
zenith <- 90 - apply(m.angle, 1, max, na.rm = T)
# Derive ternary pattern
ternary.pattern <- brks.ind[findInterval(nadir - zenith, brks)]
plus.ind <- length(which(ternary.pattern == 1))
neg.ind <- length(which(ternary.pattern == -1))
# Look up ternarity pattern and assign landform class
geomorph.lut.num[neg.ind + 1, plus.ind + 1]
}
```