R/dispersal.R

In matthewkling/stranger: Demographic species range modeling with dispersal

# some alternative dispersal kernel functions

#' @export
dlognormal <- function(x, L, S) (1/((2*pi)^1.5 * S * x^2)) * exp(-(log(x/L))^2 / (2 * S^2))

#' @export
d2Dt <- function(x, L, S) S / (pi * L * (1 + (x^2 / L))^(S + 1))

#' @export
dexponential <- function(x, L) dexp(x, rate = L) / (2*pi*x)



# function to aggregate matrices, for internal use
agg <- function(x, res, fun = sum, ...){
  y <- matrix(NA, nrow(x) / res, ncol(x) / res)
  yi <- rep(1:nrow(y), each = res)
  yj <- rep(1:ncol(y), each = res)
  for(i in 1:nrow(y)){
    for(j in 1:ncol(y)){
      y[i,j] <- fun(x[which(yi == i), which(yj == j)], ...)
    }
  }
  return(y)
}


#' Neighborhood dispersal probability matrix.
#'
#' This function uses numerical integration of a dispersal kernel to construct
#' a matrix representing the probability that dispersal originating in a given
#' grid cell will arrive at each cell across its neighborhood.
#'
#' @param kernel A dispersal kernel object, e.g. generated by \code{species_template()$kernel}.
#' Note that the kernel function must already be area-adjusted by incorporating a denominator of 2*pi,
#' as has been done for the included functions \code{dlognormal}, \code{d2Dt}, and \code{dexponential}.
#' Distance units are assumed to be in meters.
#' @param diameter Neighborhood size, in grid cells (odd integer).
#' @param cell_res Grid cell size, in meters.
#' @param method Either "area" (default), "area-centroid", or "centroid"; see Chipperfield et al. (2011).
#' @param res Resolution for numerical integration (odd integer).
#' @return A matrix of dispersal probabilities.
#' @export
#' @importFrom stats integrate
neighborhood <- function(kernel, diameter = 7, cell_res, method = "area", res = 11){

  # kernel density function
  kdf <- function(x){
    kernel$params$x <- x
    do.call(kernel$fun, kernel$params)
  }

  radius <- (diameter - 1) / 2 # radius of neighborhood (cells)
  d <- expand.grid(x = -radius:radius, y = -radius:radius)
  d$d <- sqrt(d$x^2 + d$y^2)
  for(z in unique(d$d)) d$u[d$d == z][1] <- z
  crs <- d[!is.na(d$u),]

  r <- (res - 1) / 2 # radius of coarse cell, in fine units
  if(method == "centoid") r <- 0

  oo <- (-r:r) / res # offset
  od <- (-r:r) / res # offset
  if(method == "mixed") oo <- 0

  for(i in 1:nrow(crs)){
    x <- expand.grid(xo = oo, yo = oo, xd = od + crs$x[i], yd = od + crs$y[i])
    dist <- sqrt((x$xd - x$xo) ^ 2 + (x$yd - x$yo) ^ 2)
    cdist <- sqrt(x$xd ^ 2 + x$yd ^ 2)
    dist <- ifelse(cdist > diameter / 2, Inf, dist) # clip to circle
    ci <- d$d == crs$d[i]
    d$p[ci] <- sum(kdf(dist * cell_res), na.rm = T)
    d$edge[ci] <- ! mean(dist == Inf) %in% 0:1
  }

  # reassign probability mass that falls out of bounds to boundary cells
  # (even allocation, which is imperfect; arc length would be better)
  ib <- integrate(function(x) kdf(x) * 2 * pi * x,
                  0, diameter / 2 * cell_res)$value
  d$ob <- 0
  d$ob[d$edge] <- (1 - ib) / sum(d$edge)

  # normalize probabilities, sum inbounds and outbounds, format
  d$ib <- d$p / sum(d$p) * ib
  d$total <- d$ib + d$ob
  matrix(d$total, diameter, diameter)
}