Makurhini: Makurhini: Analyzing landscape connectivity

Documented in MK_RMCentrality

#' Estimate radial and medial centrality measures.
#'
#' Use this function to calculate radial and medial centrality measure under one or several distance thresholds.
#' @param nodes Object containing nodes (e.g., habitat patches or fragments) information. It can be of the following classes:\cr
#' -   \code{Data.frame} with at least two columns: the first for node IDs and the second for attributes. If the `restoration` argument is used, the data frame must include a third column for restoration values.\cr
#' -   Spatial data of type vector (class \code{sf, SpatVector, SpatialPolygonsDataFrame}). It must be in a projected coordinate system.\cr
#' -   Raster (class \code{RasterLayer, SpatRaster}). It must be in a projected coordinate system. The values must be integers representing the ID of each habitat patch or node, with non-habitat areas represented by NA values (see \link[raster]{clump} or \link[terra]{patches}).
#' @param distance  A \code{matrix} or \code{list} to establish the distance between each pair of nodes. Distance between nodes may be Euclidean distances (straight-line distance) or effective distances (cost distances) by considering the landscape resistance to the species movements. If it is a matrix, then the number of columns and rows must be equal to the number of nodes. This distance matrix could be generated by the \link[Makurhini]{distancefile} function.  If it is a list, then it must contain the distance parameters necessary to calculate the distance between nodes. For example, two of the most important parameters: \code{“type”} and \code{“resistance”}. For \code{"type"} choose one  of the distances:  \bold{"centroid" (faster), "edge", "least-cost" or "commute-time"}. If the type is equal to \code{"least-cost"} or \code{"commute-time"}, then you must use the \code{"resistance"} argument. To see more arguments see the \link[Makurhini]{distancefile} function.
#' @param distance_thresholds A \code{numeric} indicating the dispersal distance or distances (meters) of the considered species. If \code{NULL} then distance is estimated as the median dispersal distance between nodes. Alternatively, the \link[Makurhini]{dispersal_distance} function can be used to estimate the dispersal distance using the species home range.
#' @param binary \code{logical}. Binary metrics, it only considers the distance thresholds to establish if a pair of nodes is (1) or not connected (0). Probability argument is not necessary.
#' @param probability A \code{numeric} value indicating the probability that corresponds to the distance specified in the \code{distance_threshold}. For example, if the \code{distance_threshold} is a median dispersal distance, use a probability of 0.5 (50\%). If the \code{distance_threshold} is a maximum dispersal distance, set a probability of 0.05 (5\%) or 0.01 (1\%). If \code{probability = NULL}, then a probability of 0.5 will be used.
#' @param rasterparallel \code{logical}. If nodes is "raster" then you can use this argument to assign the metrics values to the nodes raster. It is useful when raster resolution is less than 100 m<sup>2</sup>.
#' @param write \code{character}. Write output \code{sf} object. It is necessary to specify the "Folder direction"
#'  + "Initial prefix",  for example, \code{"C:/ejemplo"}.
#' @param intern \code{logical}. Show the progress of the process, \code{default = TRUE}. Sometimes the advance process does not reach 100 percent when operations are carried out very quickly.
#' @details This function implements Patch-Scale Connectivity or Centrality Measures. Radial measures: degree, strength (using probability argument, for weighted graphs),
#' eigenvector centrality (eigen), and closeness centrality (close). Medial measures: betweenness centrality (BWC),
#' node memberships (cluster), and modularity (modules, using probability argument).
#' The function builds on functions out of Csardi’s ’igraph’ package.
#' @references
#' Borgatti, S. P., & Everett, M. G. (2006). A Graph-theoretic perspective on centrality. Social Networks, 28(4), 466–484. https://doi.org/10.1016/j.socnet.2005.11.005\cr
#' Hanski, I. and Ovaskainen, O. 2000. The metapopulation capacity of a fragmented landscape. Nature 404: 755–758.
#' @export
#' @examples
#' \dontrun{
#' library(Makurhini)
#' library(sf)
#' data("habitat_nodes", package = "Makurhini")
#' nrow(habitat_nodes) # Number of patches
#' #Two distance threshold,
#' centrality_test <- MK_RMCentrality(nodes = habitat_nodes,
#'                                 distance = list(type = "centroid"),
#'                                  distance_thresholds = c(10000, 100000),
#'                                  probability = 0.05,
#'                                  write = NULL)
#' centrality_test
#' plot(centrality_test$d10000["degree"], breaks = "jenks")
#' plot(centrality_test$d10000["BWC"], breaks = "jenks")
#' plot(centrality_test$d10000["cluster"], breaks = "jenks")
#' plot(centrality_test$d10000["modules"], breaks = "jenks")
#' }
#' @importFrom magrittr %>%
#' @importFrom sf write_sf st_as_sf
#' @importFrom igraph graph_from_adjacency_matrix strength eigen_centrality closeness betweenness cluster_louvain degree components
#' @importFrom raster as.matrix raster stack
#' @importFrom future multicore multisession plan availableCores
#' @importFrom furrr future_map
#' @importFrom purrr map
#' @importFrom terra rast unique writeRaster crs crs<- classify
MK_RMCentrality <- function(nodes,
                            distance = list(type = "centroid"),
                            distance_thresholds = NULL,
                            binary = TRUE,
                            probability = NULL,
                            rasterparallel = FALSE,
                            write = NULL, intern = TRUE){
  if(missing(nodes)) {
    stop("error missing file of nodes")
  } else {
    if (is.numeric(nodes) | is.character(nodes)) {
      stop("error missing file of nodes")
    }
  }

  if(is.null(distance_thresholds)) {
    stop("error missing numeric distance threshold(s)")
  }

  if(!is.null(probability) & isTRUE(binary)){
    binary <- FALSE
  }

  if(!is.null(write)) {
    if (!dir.exists(dirname(write))) {
      stop("error, output folder does not exist")
    }
  }

  if(class(nodes)[1] == "SpatialPolygonsDataFrame"| class(nodes)[1] == "sf"){
    if(nrow(nodes) < 2) {
      stop("error, you need more than 2 nodes")
    } else {
      nodes <- st_as_sf(nodes); nodes$IdTemp <- 1:nrow(nodes); idT <- "IdTemp"
    }
    nodes.2 <- nodes
  } else {
    rcl <- class(nodes)[1]
    if( rcl == "RasterLayer"){
      nodes <- terra::rast(nodes)
    }
    t_uniq <- terra::unique(nodes)
    if(nrow(t_uniq) < 2){
      stop("error, you need more than 2 nodes")
    } else {
      idT <- NULL
    }
    nodes.2 <- t_uniq; names(nodes.2) <- "IdTemp"
  }
  if(isTRUE(intern)){
    if(!is.null(distance$resistance)){
      message("Estimating distances. This may take several minutes depending on the number of nodes and raster resolution")
    } else {
      if(nrow(nodes) > 1000){
        message("Estimating distances. This may take several minutes depending on the number of nodes")
      }
    }
  }
  dist <- distancefile(nodes = nodes,  id = idT, type = distance$type,
                       distance_unit = distance$distance_unit, keep = distance$keep,
                       resistance = distance$resistance,
                       resist.units = distance$resist.units,
                       CostFun = distance$CostFun, ngh = distance$ngh,
                       mask = distance$mask,
                       threshold = distance$threshold,
                       geometry_out = distance$geometry_out,
                       bounding_circles = distance$bounding_circles,
                       parallel = distance$parallel,
                       least_cost.java = distance$least_cost.java,
                       cores.java = distance$cores.java, ram.java = distance$ram.java,
                       pairwise = FALSE,
                       write = NULL)
  if(!is.null(idT)){
    rownames(dist) <- nodes$IdTemp; colnames(dist) <- nodes$IdTemp
  } else {
    rownames(dist) <- t_uniq[[1]]; colnames(dist) <- t_uniq[[1]]
  }

  centrality_result <- map(1:length(distance_thresholds), function(x){
    x <- distance_thresholds[x]
    if(isTRUE(binary)){
      Adj_matr <- dist*0; Adj_matr[dist < x] <- 1; diag(Adj_matr) <- 0
      graph_nodes <- tryCatch(igraph::graph_from_adjacency_matrix(Adj_matr, mode = "undirected"), error = function(err) err)
    } else {
      Adj_matr <- dist*0
      if(is.null(probability)){
        k =(1 / x); Adj_matr <- exp(-k * dist)
      } else {
        Adj_matr <- exp((dist * log(probability))/x)
      }
      diag(Adj_matr) <- 0
      graph_nodes <- tryCatch(igraph::graph_from_adjacency_matrix(Adj_matr, mode = "undirected", weighted = TRUE), error = function(err) err)
    }

    if (inherits(graph_nodes, "error")) {
      stop("graph adjacency error")
    }

    if(isFALSE(binary)){
      metric.strength <- strength(graph_nodes, weights = 1 / E(graph_nodes)$weight)
      metric.eigen <- eigen_centrality(graph_nodes, weights = 1 / E(graph_nodes)$weight)
      metric.close <- closeness(graph_nodes, weights = 1 / E(graph_nodes)$weight, normalized = TRUE)
      metric.between <- betweenness(graph_nodes, weights = 1 / E(graph_nodes)$weight)
      metric.Membcomponents <- components(graph_nodes)$membership
      metric.modularity <- cluster_louvain(graph_nodes, weights = 1 / E(graph_nodes)$weight)
      modules <- rep(0, nrow(dist))
      for(i in 1:length(metric.modularity)){
        n <- metric.modularity[[i]] %>% as.numeric()
        modules[which(nodes.2$IdTemp %in% n)] <- i
      }

      metric_conn <- cbind(rownames(dist), metric.strength, metric.eigen$vector, metric.close,
                           metric.between, metric.Membcomponents,
                           modules) %>% as.data.frame()
      names(metric_conn) <- c("id", "strength", "eigen", "close", "BWC", "cluster", "modules")
    } else {
      metric.degree <- degree(graph_nodes); metric.eigen <- eigen_centrality(graph_nodes)
      metric.close <- closeness(graph_nodes); metric.between <- betweenness(graph_nodes)
      metric.Membcomponents <- components(graph_nodes)$membership; metric.modularity <- cluster_louvain(graph_nodes)
      modules <- rep(0, nrow(dist))
      for(i in 1:length(metric.modularity)){
        n <- metric.modularity[[i]] %>% as.numeric()
        modules[which(nodes.2$IdTemp %in% n)] <- i
      }
      metric_conn <- cbind(rownames(dist), metric.degree, metric.eigen$vector, metric.close,
                           metric.between, metric.Membcomponents,
                           modules) %>% as.data.frame()
      names(metric_conn) <- c("id", "degree", "eigen", "close", "BWC", "cluster", "modules")
    }

    if(!is.null(idT)){
      if(isFALSE(binary)){
        nodes.2$strength <- metric.strength
      } else {
        nodes.2$degree <- metric.degree
      }
      nodes.2$eigen <- metric.eigen$vector
      nodes.2$close <- metric.close; nodes.2$BWC <- metric.between
      nodes.2$cluster <- metric.Membcomponents; nodes.2$modules <- modules
      nodes.2 <- nodes.2[moveme(names(nodes.2), "geometry last")]
      nodes.2$IdTemp <- NULL
      if(!is.null(write)){
        write_sf(nodes.2, paste0(write, "_", "d", x,  ".shp"), delete_layer = TRUE)
      }
    } else {
      rp <- unique(nodes)[[1]]; rp <- as.vector(rp); rp <- rp[which(!is.na(rp))]
      if(isTRUE(rasterparallel)){
        if(class(nodes)[1] == "SpatRaster"){
          nodes <- raster(nodes)
        }
        m <- matrix(nrow = nrow(dist), ncol = 2); m[,1] <- rownames(dist) %>% as.numeric()
        works <- as.numeric(availableCores())-1
        if(.Platform$OS.type == "unix") {
          strat <- future::multicore
        } else {
          strat <- future::multisession
        }
        plan(strategy = strat, gc = TRUE, workers = works)
        r_metric <- tryCatch(future_map(2:ncol(metric_conn), function(c){
          x1 <- metric_conn[, c(1, c)]
          for(i in rp){
            n <- x1[[which(x1[,1]== i), 2]] %>% as.character() %>% as.numeric()
            m[which(m[,1] == i),2] <- n
          }
          x1 <- reclassify(nodes, rcl = m)
          return(x1)}, .progress = TRUE), error = function(err) err)
        close_multiprocess(works)
        r_metric <- lapply(r_metric, rast); nodes <- rast(nodes)
      } else {
        m <- matrix(nrow = nrow(dist), ncol = 2); m[,1] <- rownames(dist) %>% as.numeric()
        r_metric <- lapply(2:ncol(metric_conn), function(c){
          x1 <- metric_conn[, c(1, c)]
          for(i in rp){
            n <- x1[[which(x1[,1]== i), 2]] %>% as.character() %>% as.numeric()
            m[which(m[,1] == i),2] <- n
          }
          x1 <- classify(nodes, rcl = m)
          return(x1)})
      }

      nodes.2 <- list()
      for(i in 2:(length(r_metric)+1)){
        nodes.2[[i]] <- r_metric[[i-1]]
      }

      nodes.2[[1]] <- nodes; nodes.2 <- rast(nodes.2)
      names(nodes.2) <- names(metric_conn)

      if (!is.null(write)){
        n <- names(metric_conn)
        n <- lapply(as.list(2:length(n)), function(w){
          x1 <- nodes.2[[w]]; crs(x1) <- crs(nodes.2)
          writeRaster(x1, filename = paste0(write, "_", n[w], "_",  x, ".tif"), overwrite = TRUE, options = c("COMPRESS=LZW", "TFW=YES"))
        })
      }

      if(rcl != "SpatRaster"){
        nodes.2 <- lapply(nodes.2, raster) %>% stack()
      }
    }
    return(nodes.2) }, .progress = intern)

  if (length(distance_thresholds) == 1){
    centrality_result <- centrality_result[[1]]
  } else {
    names(centrality_result) <- paste0("d", distance_thresholds)
  }
  if(isTRUE(intern)){
    message("Done!")
  }
  return(centrality_result)
}