vignettes/articles/ghp/scapesClassification_02_3_PKS.Rmd
scapesClassification_02_3_PKS.RmdPeak 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.
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 tableWe 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)## Cell bathymetry local_bpi regional_bpi slope ISU
## 1 1 -1718.056 -4 129 1.410384 NA
## 2 2 -1737.816 -3 107 2.021313 NA
## 3 3 -1755.392 -1 87 2.482018 NA
## 4 4 -1787.728 -16 52 3.149354 NA
## 5 5 -1819.800 -33 18 2.125063 NA
## 6 6 -1806.728 -7 28 1.895863 NAIn 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).
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