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

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

Peak cells are local maxima on a raster surface (Figure 1). In scapesClassification, a cell constitutes a local maximum if its elevation value is larger than the values of all the cells in its neighborhood.

Load data

Libraries and data are loaded and processed as explained in format input data.

# 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)

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

# COMPUTE NEIGHBORHOOD LIST
nbs <- ngbList(rstack, rNumb = TRUE, attTbl = atbl) # neighbors are identified by their row number in the attribute table

Class vector filters

We will also load the class vector of the island shelf unit (ISU), add it to the attribute table and use it as a classification filter to exclude ISU areas from peak search areas.

# ADD ISU FILTER
ISU <- list.files(system.file("extdata", package = "scapesClassification"), full.names = T,
                  pattern = "ISU_Cells\\.RDS")

atbl$ISU <- readRDS(ISU)
atbl$ISU[!is.na(atbl$ISU)] <- 1 # Class 1 identifies all ISU cells 

head(atbl)

Plots

In the following example we will show how the class vector of peak cells is computed. However, in order to improve the reading experience, the plots' code is hidden. It can be accessed in the *.RMD file used to generate the html file.

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 example we will use 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).

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)

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

# 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")

Peak cell

We will identify peak cells on our bathymetric surface using the function peak.cell(). In a second step, we will use the function cond.reclass() to filter out (i) peaks within the island shelf unit and (ii) peaks that are not on prominent relieves.

• peak.cell: finds local maxima or minima on a raster surface. In this example, we will find local maxima in the bathymetry column of the attribute table. Cells on the edge of the raster and cells adjacent to NA-cell are excluded from the peak search area (p_edge = FALSE).

• cond.reclass: evaluates conditions for cells of a class and reclassifies them if conditions are true. In this examples peak cells not meeting the filter condition are reclassified as NA-cells and excluded from the peak class vector.

# PEAK CELL
atbl$pc <- peak.cell(atbl, nbs, rNumb = TRUE,
                     p_col = "bathymetry", # column on which local maxima are searched
                     p_edge = FALSE)       # local maxima or minima are not searched on edge cells

# FILTER PEAK CELL

# Remove:
# (i)  peak cell within the ISU            -> !is.na(ISU)
# (ii) peak cell not on prominent features -> BPI < 100

cond <- "regional_bpi < 100 | local_bpi < 100 | !is.na(ISU)"

# Filter
atbl$pc <- cond.reclass(atbl, nbs, rNumb = TRUE, classVector = atbl$pc, # filter class vector atbl$pc
                        cond = cond,  # filter condition
                        class = 1,    # peak cells meeting conditions...
                        reclass = NA) # are reclassified as NA-cells

# Land cells
lc <- readRDS(list.files(system.file("extdata", package = "scapesClassification"), 
                         full.names = T, pattern = "Land_Cells.RDS"))

rL <- rstack[[1]]; rL[]  <- NA; rL[lc] <- 1
rL <- as.polygons(rL, dissolve = TRUE)

# Peak cells
rPC <- cv.2.rast(rstack, atbl$pc, atbl$Cell)
rISU<- cv.2.rast(rstack, atbl$ISU, atbl$Cell)
rPC  <- as.polygons(rPC, dissolve = TRUE)
rISU <- as.polygons(rISU, dissolve = TRUE)

# Map
rL$dummy   <- 1
rPC$dummy  <- 1
rISU$dummy <- 1

m <-
  mapview(raster::raster(hs), col.regions = c("gray0", "gray50", "grey100"), alpha.regions=1, 
          legend = FALSE, label = FALSE) +
  mapview(as(rL,"Spatial"), alpha.regions = 1, col.regions = "black", popup = NULL, 
          legend = FALSE, label = FALSE) +
  mapview(as(rISU,"Spatial"), alpha.regions = 1, 
          col.regions = "gray", layer.name = "ISU cell", popup = NULL, label = "ISU cell") + 
  mapview(as(rPC,"Spatial"), alpha.regions = 0.8, 
          col.regions = "black", layer.name = "Peak cell", popup = NULL, label = "Peak cell") + 
  mapview(raster::raster(rstack[["bathymetry_exm"]]), layer.name = "Depth (m)", col.regions = palRYB) 

m@map %>% 
  addLayersControl(overlayGroups = c("ISU cell", "Peak cell", "Depth (m)"),
                   position = "topleft",
                   options = layersControlOptions(collapsed = TRUE))

**Figure 1 - Peak cells,** _local maxima_ found on the bathymetric raster layer. Bathymetry values are provided in background (due to plot re-projection, displayed and original values might be slightly different).

Final class vector

The peak class vector is saved in the file 'PKS_Cells.RDS' and shared as part of the package data.

# # PEAK CELLS
# unique(atbl$PKS)
# 
# # atbl$PKS == 1,  Peak cells
# # atbl$PKS == NA, Unclassified cells

ghTaranto/scapesClassification documentation built on April 8, 2022, 1:43 p.m.