R/csd_estimate.R

Defines functions csd_estimate

Documented in csd_estimate

#' csd_estimate: Estimate the connectivity suitability and dispersal plot
#' @description csd_plot gives an estimate of the number of geographic clusters
#' given a set of dispersal hypothesis and a suitability raster
#'
#' @param model A raster model or a setA object representing the suitability model
#' @param dispersal_steps A numeric vector with elements representing the dispersal
#' hypothesis to test.
#' @author Luis Osorio-Olvera & Jorge Soberon
#' @examples
#' \dontrun{
#' model_path <- system.file("extdata/Lepus_californicus_cont.tif",
#'                           package = "bam")
#' model <- raster::raster(model_path)
#' model <- model > 0.7
#' csd_plot <- csd_estimate(model,
#'                          dispersal_steps=c(2,4,8,16,32,45))
#'
#'
#' }
#' @export

csd_estimate <- function(model,dispersal_steps=c(2,4,8,16,32,64)){

  if(methods::is(model,"RasterLayer")){
    model <- bam::model2sparse(model)
  }
  dispersal_steps <- sort(dispersal_steps)
  ds <- seq_along(dispersal_steps)
  pb <- utils::txtProgressBar(min = 0, max = max(ds), style = 3)
  testworks <- TRUE
  csd <- ds %>% purrr::map(function(x){
    if(testworks){
      bclust <- try({
        r <- bam::bam_clusters(model,ngbs = dispersal_steps[x])
      },silent = TRUE)

      if(!methods::is(bclust,"csd")){
        warning("There is not enough memory to calculate",
                "the adjacency matrix for dispersal step =",
                dispersal_steps[x],"\n returning ",
                paste(dispersal_steps[1:(x-1)],collapse = " "))
        testworks <- FALSE
      }
    }
    utils::setTxtProgressBar(pb, x)
    return(bclust)
  })



  ndisp <- length(csd)
  d_all <- 1:ndisp %>% purrr::map_df(function(x){
    if(methods::is(csd[[x]],"csd")){
      d1 <- data.frame(csd[[x]]@connections,
                       d=dispersal_steps[x])
      return(d1)
    }
  })

  #d_clust <- d_all %>%
  #  dplyr::group_by(d) %>%
  #  dplyr::summarise(n_clusters=max(clusterID))
  d <- NULL; clusterID <- NULL; nclust <-NULL
  d_clust <- d_all  %>% dplyr::group_by(d,clusterID) %>%
    dplyr::summarise(nclust=dplyr::n()) %>% dplyr::group_by(d) %>%
    dplyr::summarise(Clusters=max(clusterID),
                     mean_area=mean(nclust))

  s1 <- as.factor(d_clust$mean_area)
  nlevs <- length(unique(d_clust$mean_area))
  levels(s1) <- seq(1,2,length.out = nlevs)
  #data.frame(f1=s1,
  #           f2=d_clust$mean_area)
  sizes <- as.numeric(as.character(s1))

  graphics::plot(d_clust$d,d_clust$Clusters,frame=TRUE,
                 type="l",pch=19,lwd=2,cex=2,
                 xlab="Neighbors",ylab="Clusters",
                 cex.lab=1.25, cex.axis=1.25)
  graphics::points(d_clust$d,d_clust$Clusters,pch=19,
                   cex=sizes*1.5,col="brown")


  #sizes <- ( d_clust$mean_area - min(d_clust$mean_area))/
  # (max(d_clust$mean_area)-min(d_clust$mean_area))
  #sizes[1] <- sizes[2]/2
  #sizes <- sizes*5


  graphics::legend("topright",
                   legend = paste(round(d_clust$mean_area)),
                   col = "brown",
                   pch = 19,
                   bty = "n",
                   pt.cex = sizes,
                   cex = 1,
                   title="MNCC",
                   text.col = "black",
                   horiz = F ,
                   inset = c(0.1, 0.1))

  p <- grDevices::recordPlot()
  #graphics::plot.new() ## clean up device
  #
  names(csd) <- paste0("dispersal_step:_",dispersal_steps)
  csd <- list(csd=csd,plot_data=d_clust,plot=p)


  return(csd)

}
luismurao/bam documentation built on Nov. 28, 2022, 3:02 p.m.