R/find_secDiff.R

In nbelouard/slfjumps: Identify jump dispersal in biological invasions

#'Find secondary diffusion 
#'
#'Find secondary diffusion points after having applied find_jumps()
#'
#'@export
#'
#'@param potDiffusion object exported from find_jumps
#'@param Jumps object exported from find_jumps
#'@param Dist object exported from find_jumps
#'@param gap_size Distance necessary for a point to be considered a jump, 
#'in kilometers (default: 15)
#'@param crs Coordinate system for spatial layers
#'
#'@return A data frame \code{SecDiff} containing the secondary diffusion points, a data frame \code{NewDist} containing the final convex hull vertex of the diffusive spread.
#' 
#'@examples
#' secDiff <- find_secDiff()


find_secDiff <- function(potDiffusion = potDiffusion,
                         Jumps = Jumps,
                         diffusers = diffusers, 
                         Dist = preDist, 
                       gap_size = 15,
                       crs = 4326) {
  
  cat(paste(Sys.time(), "Start finding secondary diffusion... "))
  
  potDiffusion %<>% arrange(year) # sort potential secondary diffusion points by year
  yearsDiffusion = unique(potDiffusion$year) # years where potential secondary diffusion points are identified
  newDist = data.frame()

    # for each year, we will look for points coming from secondary diffusion
    # for each year where there is secDiff, we take the year before (because if no secDiff, no need to)
    for (i in 1:length(yearsDiffusion)){
      cat(paste("Year", yearsDiffusion[i], "..."))
      # select jumps up to that year
      jumps_i = Jumps %>% filter(between(year, yearsDiffusion[1]-1, 
                                         yearsDiffusion[i]-1))
  
      # if there are jumps in this selection, continue, else go to the next year
      if (dim(jumps_i)[1] > 0){
        # select potential secDiff points that year
        potDiffusion_i = potDiffusion %>% filter(year == yearsDiffusion[i]) %>% mutate(DistToSLF = NA)
        newThresholds_i = data.frame(1)
        
        while(dim(newThresholds_i)[1] > 0 & dim(potDiffusion_i)[1] > 0){
          # convert to a layer
          potDiffusion_layeri <- st_as_sf(x = potDiffusion_i, coords = c("longitude", "latitude"), 
                                       crs = crs, remove = F)
          # select diffusion points up to that year
          diffusers_layeri <- st_as_sf(x = diffusers %>% filter(between(year, yearsDiffusion[1]-1, 
                                                            yearsDiffusion[i])), 
                                  coords = c("longitude", "latitude"), 
                      crs = crs, remove = F)
          
        for (point in 1:length(potDiffusion_layeri$DistToIntro)){ # for each potential diffusion
          rangeDist <- potDiffusion_layeri[point,] %>% pull(DistToIntro)
          diffusers_pointi <- diffusers_layeri %>% filter(between(DistToIntro, rangeDist - gap_size, rangeDist + gap_size))
          # from all diffusers, we subselect those that are less than a gap size away when looking at DistToIntro to accelerate the analysis
          
          if (dim(diffusers_pointi)[1] > 0){
            pairwise_dist <- st_distance(x = potDiffusion_layeri[point,], 
                                         y = diffusers_pointi) # calculate the distance to close diffusers 
            potDiffusion_layeri$DistToSLF[point] <- min(pairwise_dist) # add the minimal distance to the table
          } else {
            potDiffusion_layeri$DistToSLF[point] <- gap_size*1000 # if no point is close, indicate the gap size
          } 
        }
        
        # Select points less than <gap size> away from all other points, these are new thresholds
        st_geometry(potDiffusion_layeri) <- NULL
        newThresholds_i = potDiffusion_layeri %>% filter(DistToSLF < gap_size*1000) # potDiffusion that are close to Thresholds
        # add these points to the diffusers' dataset
        diffusers <- bind_rows(diffusers, newThresholds_i)
        newDist = bind_rows(newDist, newThresholds_i)
        
        potDiffusion_i = potDiffusion_layeri %>% filter(DistToSLF >= gap_size*1000) # potDiffusion pruned of newThresholds
  
        }
      } else { # if there are no jumps the year before, there cannot be any secondary diffusion, all points are diffusion
        newThresholds_i = potDiffusion %>% filter(year == yearsDiffusion[i]) 
        diffusers <- bind_rows(diffusers, newThresholds_i)
        newDist = bind_rows(newDist, newThresholds_i)
        }
    }
  
    Dist = bind_rows(Dist, newDist)
    secDiff = dplyr::setdiff(potDiffusion, newDist %>% select(-DistToSLF))

    # select the results to return
    results <- list("secDiff" = secDiff, "Dist" = Dist)
    
    cat("Analysis of secondary diffusion done. \n")
    
    return(results)
}