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)
# this holds a list of the unique values in the GM raster (i.e. 1,2,3,4,5). 
land_cover = levels(GM)[[1]]

Land cover reclassification to obtain woodland patches

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

#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(land_cover$ID,gmW) 

#inspect
reclass.mat

       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 the LCM (from Eycott et al. 2011 - 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)

Subsetting the woodland class and converting to polygons for connectivity analysis

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

We will use a 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).

#Create patches from woodland raster
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 10 hectares and over.

#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)

Cost distance analysis

Our next job is to create a transition matrix to model the cost of moving between cells in the GM layer. We do this using the transition() function in the gdistance package.

We’ll choose a 4-neighbour rule this week. 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 all woodland patches by passing their coordinates to the costDistance() function.

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

Next we build a matrix from the cost values and multiply these cost distances by 30 (the resolution of the raster) in order to set the values in the same units (metres) as the dispersal distance that we will use to parameterise movement probability in the next step. The idea here is that if the cost of movement between two cells (using a four neighbour rule) is 1, then this is simply the actual distance between cell centres (30 metres). If it is 10, then the cost distance should be equivalent to 300 metres. Note, that this calculation becomes much more complicated if we use an 8-neighbour rule. Can you tell me why?

distMat<-as.matrix(cost_dist)

distMat=distMat*30

Next we can create a simple binary adjacency matrix and set all patch distances lower than the anticipated dispersal distance (10km) to 1 and everything else to 0.

#Create a binary adjacency matrix (max dispersal capacity as 10km after Dennis et al. 2024)
binaryMat <- matrix(0, nrow=nrow(distMat), ncol=ncol(distMat))
#patches closer than the dispersal distance get 1
binaryMat[distMat < 10000] <- 1
#everything else gets 0
diag(binaryMat) <- 0

Building graphs in the igraph package

Now we can build a graph (an object with nodes and edges) from this binary adjacency matrix with the igraph package.

library(igraph)
# Build a binary graph from this matrix making sure to set weighted = NULL so that the function knows we are not dealing with probabilities. 
binaryGraph <- graph_from_adjacency_matrix(binaryMat, mode="undirected", weighted=NULL)

#Degree centrality for binary graph
binaryDegree <- degree(binaryGraph) 

#Build a data frame of these values for plotting later. 
binaryDegreeData<-data.frame(site=1:nrow(woodCentres),coords=woodCentres,degree=binaryDegree)

We can also compute the betweenness centrality for this binary graph.

#betweenness centrality
binaryBetweenness <- betweenness(binaryGraph) 

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

It would be nice to plot these values on the study area so let’s create a points layer from the coordinates and plot them with size corresponding to their centrality. Note, you may have to raise the degree value by the power of some small number to get the relative sizes to something reasonable so that they fit on your map. I have chosen degree^0.5 below which works well on my machine but you may want to adapt this.

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

#set plotting scheme
par(mfrow=c(1,2))

#Plot the study area 
plot(GM,legend=F,main="Binary Degree",cex.main=2)

#add the patch nodes with degree values. "pch" sets the symbol type (we'll just go for a filled circle) and "cex" sites the size.
plot(spNodes,pch=21,cex=binaryDegreeData$degree^0.5,col="black", bg="white", lwd=1,add=T)

#Do the same for betweenness
plot(GM,legend=F,main="Binary Between",cex.main=2)

#add the patch nodes
plot(spNodes,pch=21,cex=binaryBetweenData$between^0.2,col="black", bg="white", lwd=1,add=T)

Building a weighted graph

Next we can build a weighted graph based on the 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.

#clear the plotting panel
dev.off()

null device 
          1 

plot(exp(alpha*0:10000))

Now we apply this approach to our actual between-patch distances and build a weighted graph (based on probability of movement).

  # 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.
  
#For weighted graphs, betweenness is calculated based on cost, not conductance or probability of movement, so we take the inverse of the weights from our weightedGraph object.
weightedBetweenness <- betweenness(weightedGraph, weights=1/E(weightedGraph)$weight) 

#collect as a data frame
weightedBetweenData<-data.frame(site=1:nrow(woodCentres),coords=woodCentres,between=weightedBetweenness)

#plot against the study area
plot(GM,legend=F,main="Betweenness",cex.main=1)

#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
plot(spNodes,pch=21,cex=weightedBetweenData$between^0.1,col="black", bg="white", lwd=1,add=T)

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)

We can set the colours for plotting the different clusters manually to ensure they stand out. First create a data frame that collects patch the coordinates with their cluster memebership, then create a data frame with a set of colours we want to use for each cluster. Then to join these two together, we can use the merge() function to join both objects together by a common field (“cluster” in the example below). Then just plot.

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

#data frame of clusters and colours
clusterCols=data.frame(clusters=1:9,cols=c("green","darkgreen","yellow","pink","purple","brown","black","red","blue"))

#merge together
clusterData=merge(clusterData,clusterCols,by="clusters")

#plot GM
plot(GM,legend=FALSE)

#add the patch nodes by creating a points object with the vect() function, setting pch to 21 for a filled point with an outline then set fill to the cluster colours we just made
plot(vect(clusterData[,3:4]),pch=21,cex=2,col="black",bg=clusterData$cols,add=T)

Caluclating probability of connectivity for the whole landscape

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 these by the values (probabilities) from the first step.

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”). This formulation has the advantage of rendering the resulting connectivity score as a percentage, which allows us to describe more intuitively the percentage area of the landscape that is connected.

  # 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[!is.na(GM)])*res(GM)[1]^2
  
  # get areas of each patch
  area <- polyWood$area
  
  #get proportion of landscape that is covered by these woodland patches (not a part of the PC calculation but just for info)
  pc_patches=sum(area)/AL*100
  
  # get product of all patch areas ij and multiply by   probabilities above
  PCmat <- outer(area,area) * pstar.mat
  #and sum them
  pcMatSum <- sum(PCmat)
  # divide the above by total area of the study squared to get the PC metric
  pcMod <- sqrt(pcMatSum) / as.numeric(AL)*100
  
  
#print the result with some explanatory text and rounding to two decimal places.  

#percent cover by woodland patches
print(paste("percentage cover by woodland patches =",round(pc_patches,2)))

[1] "percentage cover by woodland patches = 7.31"

#percent connected
print(paste("percentage of the landscape that is connected =",round(pcMod,2)))

[1] "percentage of the landscape that is connected = 1.12"

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){

  #each row is a patch
  N <- nrow(prob.matrix) 
  #empty vector to store results for each patch
  dPC <- rep(NA, N) 
  
  
  #iterate over all patches
  for (i in 1:N) {
    #remove patch.i from the probability matrix
    prob.matrix.i <- prob.matrix[-i,-i]
    #remove patch.i from the patch areas
    area.i <-area[-i]
    #create adjacency matrix
    pc.graph.i <- graph_from_adjacency_matrix(prob.matrix.i, mode="undirected", weighted=TRUE)
    #create all shortest paths between patches
    pstar.mat.i <- distances(pc.graph.i, weights= -log(E(pc.graph.i)$weight))
    # back-transform to probabilities of connectedness
    pstar.mat.i <- exp(-pstar.mat.i)
    #multiply all areas and movement probabilities together
    PCmat.i <- outer(area.i, area.i)*pstar.mat.i
    #final calculation
    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