GEOG71922: Week Five: Landscape connectivity

In today’s practical we will look at different conceptualizations of connectivity in ecology and explore a graph-theoretic approach to modelling connectivity.

First you will need to set your working directory and load in the data for todays practical (this consists of one land cover layer of Greater Manchester that you will find in the “Landscape” folder on the course Dropbox).

setwd("./GEOG71922/WeekFive")

Install the packages for today’s work (we will call the libraries on first use):

#Install packages

install.packages(c("terra","gdistance","raster","igraph")

#library for reading and working with raster objects.

library(terra)

Load in the data we need for the analysis.

#load some landcover for our analysis
GM= rast("GM_LC.tif")

#and plot
plot(GM)

Carry out some initial processing to represent the different land cover classes. Note the following few lines of codes simply repeat some of the steps taken in the Week Four practical so you only need to repeat them if necessary.

# convert to categorical raster
GM=as.factor(GM) 

#inspect levels
levels(GM)

#create an object to hold our raster levels (values)

land_cover = levels(GM)[[1]]

#create a vector of strings (words) to assign to each level

Reclassify the land cover data to isolate woodland patches (which we will tentatively consider to be “habitat” for this practical).

# this holds a list of the unique values in the GM raster (i.e. 1,2,3,4,5). These values are the first column in the GM levels object (which is a list) so we need [[1]][1]
GM_cat=levels(GM)[[1]][1]

#create a new list where all values 1-4 are set to zero and 5 (canopy) is set to one. So, basically we create a raster where trees are 1 and all other values NA.
gmW = c(NA,NA,NA,NA,1)

# create a matrix for reclassification using the cbind() function
reclass.mat = cbind(GM_cat,gmW) 

#inspect
reclass.mat

  ID gmW
1  1  NA
2  2  NA
3  3  NA
4  4  NA
5  5   1

#use the classify function passing the raster as the first argument and the reclassification matrix as the second.  
gmWood = classify(GM,reclass.mat) 

#plot it
plot(gmWood,col="darkgreen")

The next step is to create the cost layer (as we did in Week Four).

########Landscape Resistance and Least Cost Path

#create cost values for LCM (from Eycott et al. 2011)


gm_levels=levels(as.factor(GM))


#costs from Eycott et al.as adapted in Dennis et al.,2026 
costs=c(5,10,10,2,1)


#combine with land cover classes (taken from the land_cover object made earlier)
costMat=cbind(land_cover$ID,costs)

#reclassify the raster
costGM=classify(GM,costMat)

#inspect the result
plot(costGM)

Now we will set up our woodland patches that will form the network for use in the connectivity analysis.

We will use an 4-neighbour rule to delineate patches of woodland and then polygonise them so that we can more easily filter out smaller patches (we will just focus on patches >= 10 Ha).

patchWood=patches(gmWood, directions=4)

#Convert to polygons using the as.polygons() function
polyWood=as.polygons(patchWood,dissolve=T)

We will only select from the larger patches so let’s subset patches of woodland over 10 hectares.

#Calculate area so we can select only large patches.
polyWood$area=expanse(polyWood)

#Subset to woodland patches > 10 hectares
polyWood=polyWood[polyWood$area>100000,]

#get the coordinates of the patch centroids using the crds() function to get coordinates and centroids() to get patch centroids.
woodCentres=crds(centroids(polyWood))

plot(polyWood)

Our next job is to create a transition matrix to model the cost of moving between cells in the GM object. We do this using the transition() function in the gdistance package. Here we have to convert to a raster object using the raster package because gdistance still works with this older rSpatial package.

Note also we set the transitionFunction argument to 1/mean(x) because the transition function actually calculates conductance (ease of movement) whereas we want the inverse of this (cost of movement). We’ll choose an 4-neighbour rule. We do this so that later we can multiply the cost layer by the resolution of the raster to get an idea of total cost distance that is commensurate with values used for dispersal in our connectivity analysis. If we allow movement in 8 directions, this is more complicated (possible, but requires some non-trivial extra coding so we will make the simplifying assumption that cost distance is based on movement in four directions).

#load libraries
library(gdistance)
library(raster)

#Create transition layer.
land_cond = transition(raster(costGM), transitionFunction=function(x) 1/mean(x), directions=4)

Now we can compute the least cost path between the all woodland patches by passing their coordinates to the costDistance() function.

  #compute cost distance between all patches
 
cost_dist <- costDistance(land_cond, woodCentres)

distMat<-as.matrix(cost_dist)

distMat=distMat*30

#Create a binary adjacency matrix (max dispersal capacity as 10km after Dennis et al. 2024)
binaryMat <- matrix(0, nrow=nrow(distMat), ncol=ncol(distMat))
binaryMat[distMat < 10000] <- 1
diag(binaryMat) <- 0
# first build a binary graph from the A.mean matrix
binaryGraph <- graph_from_adjacency_matrix(binaryMat, mode="undirected", weighted=NULL)

#betweenness centrality for weighted graph
binaryDegree <- degree(binaryGraph) #betweenness is calculated based on shortest paths (i.e. lower values are better) so we take the inverse of the weights from our graph.Aprob object.

binaryDegreeData<-data.frame(site=1:nrow(woodCentres),coords=woodCentres,degree=binaryDegree)

We can also compute the betweenness centrality for the weighted graph.

#betweenness is calculated based on shortest paths (i.e. lower values are better) so we take the inverse of the weights from our binaryGraph object.
binaryBetweenness <- betweenness(binaryGraph, weights=1/E(binaryGraph)$weight) 

binaryBetweenData=data.frame(site=1:nrow(woodCentres),coords=woodCentres,between=binaryBetweenness)

# create vector object from the patch coordinates
spNodes<-vect(woodCentres)

par(mfrow=c(1,2))

plot(GM,legend=F,main="binary Degree",cex.main=2)
plot(spNodes,pch=21,cex=binaryDegreeData$degree*0.2,col="black", bg="white", lwd=1,add=T)#add the patch nodes


plot(GM,legend=F,main="Binary Between",cex.main=2)
plot(spNodes,pch=21,cex=binaryBetweenData$between^0.2,col="black", bg="white", lwd=1,add=T)#add the patch nodes

Next we can build a weighted graph based on use of a negative exponential kernel. Here we can set the probability of dispersal at the chosen threshold to be 0.05 (this is the standard approach in studies on functional connectivity). To do this, we build alpha (the constant in the Probability of Connectivity metric) by taking the log of 0.05 (our probability) and dividing this by the chosen dispersal distance.

  # set alpha to determine colonization probability of the species (max dispersal capacity as 10km after Dennis et al. 2024)
  alpha= log(0.05)/10000

We can plot this against a sensible value range to confirm the shape of the kernel created.

plot(exp(alpha*0:10000))

Now we apply this approach to our actual between-patch distances and build a weighted graph.

  # negative exponential of colonization kernel x distance get probability
  weightedMat <- exp(alpha * distMat)
 
  # set diag to zero (so we we don't count movement from patches to themselves)
  diag(weightedMat) <- 0
  
  # create a graph from this weighted adjacency matrix
  weightedGraph <- igraph::graph_from_adjacency_matrix(weightedMat, mode="undirected", weighted=T)
  

#betweenness centrality for weighted graph
#Note betweenness is calculated based on shortest paths (i.e. lower values are better) so we take the inverse of the weights from our graph.Aprob object.
weightedBetweenness <- betweenness(weightedGraph, weights=1/E(weightedGraph)$weight) 

betweenData<-data.frame(site=1:nrow(woodCentres),coords=woodCentres,between=weightedBetweenness)






plot(GM,legend=F,main="Betweenness",cex.main=1)
plot(spNodes,pch=21,cex=betweenData$between^0.1,col="black", bg="white", lwd=1,add=T)#add the patch nodes, set pch to 21 for a filled point with outline then set fill to white and the outline to black and lwd (line width) = 1

Another useful thing we can do is to see which patches form clusters of well-connected patches using the cluster_louvain() function.

binaryModularity <- cluster_louvain(weightedGraph)


clusterData<-data.frame(site=1:nrow(woodCentres),coords=woodCentres,clusters=binaryModularity$membership)

plot(GM,legend=FALSE)
plot(spNodes,pch=21,cex=2,col="black",bg=binaryModularity$membership, lwd=1,add=T)#add the patch nodes, set pch to 21 for a filled point with outline then set fill to white and the outline to black and lwd (line width) = 1

Next, we will calculate the probability of connectivity metric (PC) which is a metric that symbolizes the overall connectivity of entire networks. The PC metric is worked out first by calculating probability of movement between all nodes in the network. Once we have this we get the product of (i.e. multiply) the areas of all pairs of patches in the network and multiply the values (probabilities) from the first step by these area products. Dividing this by the area of the landscape squared gives us the PC metric. Note that here we are applying the modification suggested in Dennis et al. (2024) (in which we take the square root of the values before dividing by the study area) so we’ll call this pcMod (for “modified”).

  # calculate all shortest paths between nodes
  pstar.mat <- distances(weightedGraph, weights= -log(igraph::E(weightedGraph)$weight))
  
  # back-transform to probabilities of connectedness
  pstar.mat <- exp(-pstar.mat)
  
  # get study area in m2
  AL <- ncell(GM)*res(GM)[1]^2
  
  # get area vector
  area <- polyWood$area
 
  
  # sum areas from the vector
  areaSum <- sum(area)
  
  # get product of all patch areas ij and multiply by probabilities above
  PCmat <- outer(area,area) * pstar.mat
 
  pcMatSum <- sum(PCmat)
  
  # divide by total area of the study squared to get the PC metric
  pcMod <- (sqrt(pcMatSum) / as.numeric(AL))*100

The final task we will carry out is to determine the importance of each patch using a leave-one-out approach. Here we build a function that will iterate over each patch and remove that patch from the graph before calculating pcMod (which we will call pcMod.i). The percentage difference between pcMod and pcMod.i (where we pcMod.i must be < pcMod) gives us the relative importance of each patch from the point of view of connectivity.

#function for calculating dPC - contribution of each patch to PC. Here we simply remove each patch using a for loop and repeat the above function, adding the difference between PC above and each new calculation (minus patch "i") to am empty vector "dPC 



prob.connectivity <- function(prob.matrix,area,landarea){

  #dPC
  N <- nrow(prob.matrix) #each row is a patch
  dPC <- rep(NA, N) #empty vector to store results for each patch
  
  

  for (i in 1:N) {
    prob.matrix.i <- prob.matrix[-i,-i]
    area.i <-area[-i]
    pc.graph.i <- graph_from_adjacency_matrix(prob.matrix.i, mode="undirected", weighted=TRUE)
    pstar.mat.i <- distances(pc.graph.i, weights= -log(E(pc.graph.i)$weight))
    pstar.mat.i <- exp(-pstar.mat.i)
    PCmat.i <- outer(area.i, area.i)*pstar.mat.i
    
    pcMod.i <- sqrt(sum(PCmat.i))/landarea*100
    
    dPC[i] <- (pcMod-pcMod.i)/pcMod*100
  }

    return(dPC)
}

Aprob.PC <- prob.connectivity(prob.matrix=weightedMat,area=area, landarea=AL)

dataDPC=data.frame(site=1:nrow(woodCentres),coords=woodCentres,dPC=Aprob.PC)


plot(GM,legend=FALSE)
plot(spNodes,pch=21,cex=dataDPC$dPC, lwd=2,add=T,col="white",bg="blue")#add the patch nodes, set pch to 21 for a filled point with outline then set fill to white and the outline to black and lwd (line width) = 1