R/calc_dist_kernel_shp.R

In qleclerc/efficientspatial: Efficient Spatial Modelling of Epidemics

#' @title Calculates the spatial kernel value (Shapefile)
#'
#' @description Calculates the spatial kernel values between all areas in a Shapefile using their distance and
#'              population. This uses a gravity kernel structure (either matched or smooth depending on the threshold
#'              distance).
#'
#' @param dist_mat The distance matrix between all areas.
#' @param dist_c The threshold distance for the matched kernel. Note that setting this to an extremely high value
#'               will effectively make the gravity kernel a smooth kernel with a single function since the threshold
#'               distance will never be crossed.
#' @param shp_data The dataframe object containing the population data extracted from the Shapefile object
#'                 (generated by the \code{\link{prep_simulation_shp}} function)
#' @param alpha The destination population power.
#' @param p The distance power.
#' @param p2 The secondary distance power for the matched kernel.
#' @param aa The offset distance.
#' @param delta The home cell modifier (to compensate the distance within cells equal to 0).
#'
#' @return Returns one matrix object containing the spatial kernel values.
#'
#'
#' @export


calc_dist_kernel_shp = function(dist_mat, dist_c, shp_data, alpha, p, p2, aa, delta=0.3){

  N = shp_data$Population

  #set up matrix to fill in:
  dist_kernel = matrix(0, nrow=nrow(dist_mat), ncol=ncol(dist_mat))

  #value at threshold distance to ensure a continuous function:
  A = (dist_c^p2)/((1+dist_c/aa)^(p))

  #pre-threshold values:
  K = 1/((1+dist_mat/aa)^(p))

  #only adds delta if i == j (i.e. if calculating transmission kernel within an area)
  delta_diag = matrix(1, nrow=nrow(dist_mat), ncol=ncol(dist_mat))
  diag(delta_diag) = 1+delta
  K = K*delta_diag

  #post-threshold values:
  K2 = A/(dist_mat^p2)

  #CHANGE THIS: useless calculations being performed, currently calculating K and K2 for ALL distances, better to
  #only calculate the relevant one based on the distance value


  for (i in 1:length(N)) {

    for (j in 1:length(N)) {

      if(dist_mat[i,j] > dist_c){

        dist_kernel[i,j] = ((N[j]^alpha))*(K2[i,j])

      } else {

        dist_kernel[i,j] = ((N[j]^alpha))*(K[i,j])

      }

    }

  }


  #normalise so that rowSums = 1:
  dist_kernel = dist_kernel/rowSums(dist_kernel)

  return(dist_kernel)

}