Create a 2d array of CRFD hotspots - using base R only

crfd_hotspots_baseR.R

# Use the Cumulative Relative Frequency Curve hotspot method
# to identify hotspots of high values

# Ref:
## Bartolino, V., Maiorano, L., & Colloca, F. (2011).
## A frequency distribution approach to hotspot identification.
## Population Ecology, 53(2), 351-359. doi:10.1007/s10144-010-0229-2

# example array data
# could do with raster/terra, but wanted example in base r
## rT is an exact copy of r used to store hotspot locations

r <- rT <- array(data = NA, dim = c(100, 50, 10))
r[] <- rT[] <- rpois(n = (dim(r)[1]*dim(r)[2]) * dim(r)[3], lambda = 15)

# empty vector & dataframe to store thresholds
thresh <- NULL
threshDF <- data.frame()

# Iterate through each layer and identify hotspots
for (i in 1:dim(r)[3]) {
  z <- as.vector(na.omit(r[,,i]))
  A <- table(round(z, 0))
  mat <- cbind(as.numeric(A), as.numeric(names(A)))
  x <- mat[, 2] / max(mat[, 2])
  y <- cumsum(mat[, 1]) / sum(mat[, 1])
  vec <- NULL
  for (j in 1:length(y)) {
    s.ang <- (y[j + 1] - y[j]) / sqrt((x[j + 1] - x[j])^2 + (y[j + 1] - y[j])^2)
    vec <- c(vec, s.ang)
  }
  deg <- (asin(vec) * 90) / (pi / 2)
  aa <- deg[deg > 45]
  M <- na.omit(data.frame(x, y, deg, mat))
  M$Year <- i
  M$TF <- ifelse((M$deg) > 45, TRUE, FALSE)
  thr <- max(M$X2[M$TF])
  thresh <- c(thresh, thr)
  M$Thr <- thr
  message("Threshold for layer ", i, " is ", thr)
  int <- max(M$y[M$TF]) - max(M$x[M$TF])
  {plot(1, 1,
      type = "n", xlim = c(0, 1), ylim = c(0, 1),
      xlab = "Relative abundance", ylab = "Cumulative Frequency",
      main = paste0("CRFD Curve - (Layer: ", i, ")"),
      cex.lab = 1.3
    )
    lines(x, y, lwd = 2)
    abline(int, 1, lty = 2, lwd = 2)
    xint <- thr / max(z)
    abline(v = xint, lty = 2)
    text(0.5, 0.25,
      lab = paste0("hotspot threshold = ", thr), cex = 1.1
    )}
  rT[,,i][] <- ifelse(rT[,,i][] >= thr, 1, 0)
  threshDF <- rbind(threshDF, M)
  rm(z, A, mat, x, y, vec,s.ang,deg, aa,M, thr, int,j,xint)
}

# Determine the proportion of time a cell is considered a hotspot
propHS <- apply(X = rT, MARGIN = c(1:2), FUN = function(x, ...) {
  sum(x >= 1, na.rm = TRUE)/length(x)
  })
image(propHS, col = rainbow(101), axes = FALSE)