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 classifies landform elements at a range of different spatial scales with unprecedented computational efficiency.

Ternary patterns and geomorphons after Jasiewicz and Stepinkski (2013)

Figure 1: Ternary patterns and geomorphons after Jasiewicz and Stepinkski (2013)

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.

Landform Classification in Eastern Switzerland (Mount Bachtel). Data: NASA SRTM Mission.

Figure 2: Landform Classification in Eastern Switzerland (Mount Bachtel). Data: NASA SRTM Mission.

Landform Classification in Eastern Switzerland. Data: NASA SRTM Mission.

Figure 3: Landform Classification in Eastern Switzerland. Data: NASA SRTM Mission.

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")
"Volcano" data set: elevation. Source: author

Figure 4: “Volcano” data set: elevation. Source: author

"Volcano" data set: landform classes Source: author

Figure 5: “Volcano” data set: landform classes Source: author

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]  
}