In ghTaranto/scapesClassification: User-Defined Classification of Raster Surfaces

knitr::opts_chunk$set(out.width = "100%", echo = TRUE)

Study case and input data

Our study area is located in the 'Triangle' of the Azores (NE Atlantic), around the islands of Faial, Pico and São Jorge (Figure 1). Geomorphic management units (GMUs) will be identified using a SpatRaster(terra package) of 4 layers that includes bathymetry and bathymetric derivatives (Walbridge et al., 2018):

Bathymetry: depth values;
Local BPI: bathymetric position index (BPI) computed with an outer radius of 10px (ca. 5 km);
Regional BPI: bathymetric position index (BPI) computed with an outer radius of 40px (ca. 20 km);
Slope: slope values;

We will start by loading the required libraries and data into the workspace:

# LOAD LIBRARIES
library(scapesClassification)
library(terra)

# LOAD DATA
grd <- list.files(system.file("extdata", package = "scapesClassification"), full.names = T)
grd <- grd[grepl("\\.grd", grd)]
grd <- grd[!grepl("hillshade", grd)]
rstack <- rast(grd)

library(mapview)
library(leaflet)
library(leafem)

# set gerenal options 
mapviewOptions(basemaps = c("Esri.OceanBasemap", "OpenStreetMap"), na.color = "grey88", homebutton = FALSE, 
               query.position = "bottomright", legend.opacity = 0.9)

# set palettes
palRYB <- c("#2892c7", "#fafa64", "#e81014")
palYGB <- c("#ffff80" ,"#3dba65", "#0d1178")
palBYR <- c("#4575b5", "#ffffbf", "#d62f27")
gm_col <- c("grey70", "#380000", "#c80000", "#ff5014", "#fad23c", "#ffff3c", "#b4e614", "#3cfa96", "#0000ff", "#000038")

We can plot the raster stack as an interactive map using the R packages mapview and leaflet:

m <- mapview(raster::raster(rstack[[1]]), layer.name = c("Bathymetry (m)"), col.regions = palRYB) +
  mapview(raster::raster(rstack[[2]]), layer.name = c("Local BPI"), col.regions = palBYR, hide = TRUE) + 
  mapview(raster::raster(rstack[[3]]), layer.name = c("Regional BPI"), col.regions = palBYR, hide = TRUE) +
  mapview(raster::raster(rstack[[4]]), layer.name = c("Slope (deg)"), col.regions = palYGB, hide = TRUE)

m@map %>% 
  hideGroup(c("Local BPI", "Regional BPI", "Slope (deg)")) %>% 
  setView(lng = -28.55, lat = 38.45, zoom = 8) %>%
  addLayersControl(baseGroups = c("Esri.OceanBasemap", "OpenStreetMap"),
                   overlayGroups = c("Bathymetry (m)", "Local BPI", "Regional BPI", "Slope (deg)"),
                   position = "topleft",
                   options = layersControlOptions(collapsed = FALSE))

**FIGURE 1 -** Raster layers used for the classification of Geomorphic Management Units (GMUs). The geomorphon classes correspond to (1) flat, (2) summit, (3) ridge, (4) shoulder, (5) spur, (6) slope, (7) hollow, (8) footslope, (9) valley and (10) depression.

Plots

In the working example articles we will show how class vectors are computed. However, in order to improve the reading experience, the plots' code is hidden. It can be accessed in the *.RMD files used to generate the html files.

The plotting procedure is to (i) convert a class vectors into a raster using the function cv.2.rast() and to (ii) visualize the raster using R or an external software. In our examples we will use the R package terra to create static maps and the R packages mapview and leaflet to create interactive maps (note that mapview do not support terra raster objects yet, therefore, they have to be converted into raster raster objects before plotting).

Format input data

The classification process in scapesClassification begins with the computation of two elements (see format inputs):

Attribute table: the raster object converted into a data.frame. Attribute tables include only raster cells having no missing values. Function attTbl().

Classification rules are built accessing by name the variables stored in the attribute table (see ?conditions).

# COMPUTE ATTRIBUTE TABLE
atbl <- attTbl(rstack, var_names = c("bathymetry", "local_bpi", "regional_bpi", "slope"))

# VIEW THE TOP 3 ROWS OF `atbl`
# Each row corresponds to a raster cell as indicated in the column `atbl$Cell`

head(atbl, 3)

List of neighborhoods: raster cell neighborhoods. Neighborhoods are computed only for cells included in the attribute table. Function ngbList().

# COMPUTE NEIGHBORHOOD LIST
nbs <- ngbList(rstack, rNumb = TRUE, attTbl = atbl) # the neighbors are identified by row numbers (see ?ngbList)

# VIEW THE TOP ELEMENT OF `nbs`
nbs[1]

# nbs[1] reads as:
# the cell in row $`1` of `atbl` has 
# cells in `atbl` rows 2, 291 and 292 as neighbors