Create a raster of CRFD hotspots

crfd_hotspots.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


library(raster)

# example raster data
## rT is an exact copy of r used to store hotspot locations
r <- rT <- brick(nl = 10, ncols = 100, nrows = 50)
r[] <- rT[] <- rpois(n = ncell(r) * nlayers(r), lambda = 15)
r

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

# Iterate through each raster layer and identify hotspots
for (i in 1:nlayers(r)) {
  # extract raster values
  z <- as.vector(na.omit(values(r[[i]])))
  # build a frequency table
  A <- table(round(z, 0))
  mat <- cbind(as.numeric(A), as.numeric(names(A)))
  # cumulative frequency
  x <- mat[, 2] / max(mat[, 2])
  y <- cumsum(mat[, 1]) / sum(mat[, 1])
  # Calculate the tangent to the curve at each point
  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)
  }
  # in degrees
  deg <- (asin(vec) * 90) / (pi / 2)
  # attempt to identify the 45°
  aa <- deg[deg > 45]
  # Convert to dataframe
  M <- na.omit(data.frame(x, y, deg, mat))
  M$Year <- as.numeric(sub("layer.", "", names(r)[i]))
  M$TF <- ifelse((M$deg) > 45, TRUE, FALSE)
  # find value where tangent line intersects
  # Identifies highest value where angle >45deg
  thr <- max(M$X2[M$TF])
  thresh <- c(thresh, thr)
  M$Thr <- thr
  message("Threshold for layer ", sub("layer.", "", names(r)[i]), " is ", thr)
  # and intercept
  int <- max(M$y[M$TF]) - max(M$x[M$TF])
  # Plot
  {plot(1, 1,
      type = "n", xlim = c(0, 1), ylim = c(0, 1),
      xlab = "Relative abundance", ylab = "Cumulative Frequency",
      main = paste0("CRFD Curve - (Layer: ", sub("layer.", "", names(r)[i]), ")"),
      cex.lab = 1.3
    )
    lines(x, y, lwd = 2)
    abline(int, 1, lty = 2, lwd = 2)
    # identify the x intercept
    xint <- thr / max(z)
    abline(v = xint, lty = 2)
    text(0.5, 0.25,
      lab = paste0("hotspot threshold = ", thr), cex = 1.1
    )}
  # Create raster of hotspots
  rT[[i]][] <- ifelse(values(rT[[i]]) >= thr, 1, 0)
  threshDF <- rbind(threshDF, M)
  rm(z, A, mat, x, y, vec,deg, aa,M, thr, int)
}

# Determine the proportion of time a cell is considered a hotspot
propHS <- calc(rT, fun = function(x, ...) {
  sum(x >= 1, na.rm = TRUE)/length(x)
  })
propHS; plot(propHS, col = bpy.colors(101))