connectivity_matrix(x, y, ...)

## S4 method for signature 'Spatial,Raster'
connectivity_matrix(x, y, ...)

## S4 method for signature 'Spatial,character'
connectivity_matrix(x, y, ...)

## S4 method for signature 'sf,character'
connectivity_matrix(x, y, ...)

## S4 method for signature 'sf,Raster'
connectivity_matrix(x, y, ...)

## S4 method for signature 'Raster,Raster'
connectivity_matrix(x, y, ...)

# load data
data(sim_pu_raster, sim_pu_sf, sim_features)

# create connectivity matrix using raster planning unit data using
# the raster cost values to represent conductance
## extract 9 planning units
r <- crop(sim_pu_raster, c(0, 0.3, 0, 0.3))

## extract conductance data for the 9 planning units
cd <- crop(sim_features, r)

## make connectivity matrix using the habitat suitability data for the
## second feature to represent the planning unit conductance data
cm_raster <- connectivity_matrix(r, cd[[2]])

## plot data and matrix
## Not run: 
par(mfrow = c(1,3))
plot(r, main = "planning units (raster)", axes = FALSE, box = FALSE)
plot(cd[[2]], main = "conductivity", axes = FALSE, box = FALSE)
plot(clamp(raster(as.matrix(cm_raster)), lower = 1e-5, useValues = FALSE),
     main = "connectivity", axes = FALSE, box = FALSE)

## End(Not run)
# create connectivity matrix using polygon planning unit data using
# the habitat suitability data for the second feature to represent
# planning unit conductances
## subset data to 9 polygons
ply <- sim_pu_sf[c(1:2, 10:12, 20:22), ]

## make connectivity matrix
cm_ply <- connectivity_matrix(ply, sim_features[[2]])

## plot data and matrix
## Not run: 
par(mfrow = c(1, 2))
plot(st_geometry(ply), main = "planning units (sf)")
plot(clamp(raster(as.matrix(cm_ply)), lower = 1e-5, useValues = FALSE),
     main = "connectivity", axes = FALSE, box = FALSE)

## End(Not run)

# create connectivity matrix using habitat suitability data for each feature,
# this could be useful if prioritisations should spatially clump
# together adjacent planning units that have suitable habitat
# for the same species (e.g. to maintain functional connectivity)

## let's use the raster data for this example, and we can generate the
## connectivity matrix that we would use in the prioritization by
## (1) generating a connectivity matrix for each feature separately, and
## and then (2) then summing the values together
cm_sum <- lapply(as.list(cd), connectivity_matrix, x = r) # make matrices
cm_sum <- Reduce("+", cm_sum) # sum matrices together

## plot data and matrix
## Not run: 
par(mfrow = c(1, 2))
plot(r, main = "planning units (raster)", axes = FALSE, box = FALSE)
plot(clamp(raster(as.matrix(cm_sum)), lower = 1e-5, useValues = FALSE),
     main = "connectivity", axes = FALSE, box = FALSE)

## End(Not run)

## we could take this example one step further, and use weights to indicate
## relative importance of maintaining functional connectivity
## for each feature (i.e. use the weighted sum instead of the sum)

## let's pretend that the first feature is 20 times more important
## than all the other species
weights <- c(20, 1, 1, 1, 1)

## calculate connectivity matrix using weighted sum
cm_wsum <- lapply(as.list(cd), connectivity_matrix, x = r) # make matrices
cm_wsum <- Map("*", cm_wsum, weights) # multiply by weights
cm_wsum <- Reduce("+", cm_wsum) # sum matrices together

## plot data and matrix
## Not run: 
par(mfrow = c(1, 2))
plot(r, main = "planning units (raster)", axes = FALSE, box = FALSE)
plot(clamp(raster(as.matrix(cm_wsum)), lower = 1e-5, useValues = FALSE),
     main = "connectivity", axes = FALSE, box = FALSE)

## End(Not run)

## since the statistical distribution of the connectivity values
## for each feature (e.g. the mean and standard deviation of the
## connectivity values) are different, it might make sense -- depending
## on the goal of the conservation planning exercise and the underlying
## data -- to first normalize the conductance values before applying the
## weights and summing the data for feature together

## one approach would be to linearly rescale the values between 0.01 and 1
## note that we wouldn't want to rescale them between 0 and 1 since
## a value of zero means that there is no connectivity at all (and
## and not a relatively small amount of connectivity)
## Not run: 
### define helper function
library(scales) # load scales library for rescale
rescale_matrix <- function(x) {x@x <- rescale(x@x, c(0.01, 1)); x}

### calculate functional connectivity matrix using the weighted sum of
### connectivity values that have been normalized by linearly re-scaling
### values
cm_lwsum <- lapply(as.list(cd), connectivity_matrix, x = r) # make matrices
cm_lwsum <- lapply(cm_lwsum, rescale_matrix) # rescale matrices to [0.01, 1]
cm_lwsum <- Map("*", cm_lwsum, weights) # multiply by weights
cm_lwsum <- Reduce("+", cm_lwsum) # sum matrices together

## End(Not run)

## plot data and matrix
## Not run: 
par(mfrow = c(1, 2))
plot(r, main = "planning units (raster)", axes = FALSE, box = FALSE)
plot(clamp(raster(as.matrix(cm_lwsum)), lower = 1e-5, useValues = FALSE),
     main = "connectivity", axes = FALSE, box = FALSE)

## End(Not run)

## another approach for normalizing the data could be using z-scores
## note that after normalizing the data we would need to add a constant
## value so that none of the connectivity values are negative

### define helper functions
zscore <- function(x) {x@x <- (x@x - mean(x@x)) / sd(x@x); x}
min_non_zero_value <- function(x) min(x@x)
add_non_zero_value <- function(x, y) {x@x <- x@x + y; x}

### calculate functional connectivity matrix using the weighted sum of
### connectivity values that have been normalized using z-scores,
### and transformed to account for negative values
cm_zwsum <- lapply(as.list(cd), connectivity_matrix, x = r) # make matrices
cm_zwsum <- lapply(cm_zwsum, zscore) # normalize using z-scores
min_value <- min(sapply(cm_zwsum, min_non_zero_value)) # find min value
min_value <- abs(min_value) + 0.01 # prepare constant for adding to matrices
cm_zwsum <- lapply(cm_zwsum, add_non_zero_value, min_value) # add constant
cm_zwsum <- Map("*", cm_zwsum, weights) # multiply by weights
cm_zwsum <- Reduce("+", cm_zwsum) # sum matrices together

## plot data and matrix
## Not run: 
par(mfrow = c(1, 2))
plot(r, main = "planning units (raster)", axes = FALSE, box = FALSE)
plot(clamp(raster(as.matrix(cm_zwsum)), lower = 1e-5, useValues = FALSE),
     main = "connectivity", axes = FALSE, box = FALSE)

## End(Not run)