In this week’s practical we will return to a graph theoretic framework for modelling habitat connectivity but will focus more on the development of a function approach that builds on the species-oriented aspects of the previous analysis. In Week Five, this consisted of assigning cost values to the landscape based on some species specific assumptions. In this practical we will go a little deeper and consider the quantification of habitat itself from a species-oriented perspective (as opposed to the more common parameterisation of habitat as one or other land cover type).
First sort out the packages needed for this analysis. Note that there are no new packages this week so you may just choose to load all the necessary libraries instead.
The following objects are masked from 'package:terra':
intersect, union
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
library(raster)
Loading required package: sp
Attaching package: 'raster'
The following object is masked from 'package:dplyr':
select
library(gdistance)
Loading required package: igraph
Attaching package: 'igraph'
The following object is masked from 'package:raster':
union
The following objects are masked from 'package:dplyr':
as_data_frame, groups, union
The following objects are masked from 'package:terra':
blocks, compare, union
The following objects are masked from 'package:stats':
decompose, spectrum
The following object is masked from 'package:base':
union
Loading required package: Matrix
Attaching package: 'gdistance'
The following object is masked from 'package:igraph':
normalize
library(sf)
Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE
Set the weorking directory
setwd('./Week7')
Load in a classified land cover layer of the study area created from a Random Forest model.
#actual classified land cover layerclassif<-rast("classif_RF.tif")
From this classified raster we will isolate the woodland class as the designated habitat class and calculate the overall habitat connectivity score for the study area and the contribution of each patch in the same way as in the previous practical.
#isolate the woodland classwoodland=classif#set all values that are not 1 to zerowoodland[woodland!=1]=0# use the patches() function towoodland<-patches(woodland,directions=8,zeroAsNA=T)#polygonize patches for cost distance analysispatchesBool<-as.polygons(woodland,dissolve=TRUE)patchesBool$Area<-expanse(patchesBool)#subset patches >= 0.5 hapatchesBool=patchesBool[patchesBool$Area>=5000]
Next we assign costs to each class in the land cover map.
#set the coasts for each land coverlandCosts=c(1,10,2,10,5)#create a Boolean layer of landscape costbooleanCost=classify(classif,cbind(c(1,2,3,4,5),landCosts))#plotplot(booleanCost)
Given that the M62 motorway runs through the study area and we can expect this to have much greater cost than the rest of the urban class, we will set this feature an order of magintude higher than the simple “urban” class.
#get mask of motorway (high cost)m62<-vect("M62.shp")#rasterise m62 to combine with the cost analysism62Cost<-rasterize(m62,field=1,classif,background=0)
We can also set this landscape feature to zero in terms of habitat suitability.
#set the motorway cost to 100booleanCost[m62Cost!=0]=100
Now we are ready to start modelling so first we build a transition layer as in previous weeks.
library(gdistance)#create a transition layer from the Boolean cost layer.trBool=transition(raster(booleanCost),transitionFunction =function(x) 1/mean(x),directions =4)#get patch centroids for least cost path delineationsitesBool<-SpatialPoints(crds(centroids(patchesBool)))#create cost distance matrix#init list for least cost paths costMatBool <-as.matrix(costDistance(trBool,sitesBool))costMatBool=costMatBool*res(classif)[1]costMatBool=costMatBool*9# set alpha to determine colonization probability of the species (max dispersal capacity as 10km after Dennis et al. 2024)alpha=log(0.05) /10000# adjacency graphgraph.Aprob_bool <- igraph::graph_from_adjacency_matrix(costMatBool, mode="undirected", weighted=T)# negative exponential of colonization kernel x distance get probabilityweightedMatBool <-exp(alpha * costMatBool)# set diag to zero (so we we don't count movement from patches to themselves)diag(weightedMatBool) <-0# create a graph from this weighted adjacency matrixweightedGraphBool <- igraph::graph_from_adjacency_matrix(weightedMatBool, mode="undirected", weighted=T)## calculate Connectivity metric (RFH)# calculate all shortest paths between nodespstar.mat <-distances(weightedGraphBool, weights=-log(igraph::E(weightedGraphBool)$weight))# back-transform to probabilities of connectednesspstar.mat <-exp(-pstar.mat)# get study area in m2AL <-ncell(classif)*res(classif)[1]^2# get area vectorareaBool <- patchesBool$Area# get product of all patch areas ij and multiply by probabilities abovepcMat <-outer(areaBool,areaBool) * pstar.matpcMatSum <-sum(pcMat)# divide by total area of the study squared to get the PC metricpcMod <- (sqrt(pcMatSum) /as.numeric(AL))*100#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 diff_connectivity <-function(prob.matrix,area,baseline,landarea){#each row is a patch N <-nrow(prob.matrix) #empty vector to store results for each patch diffCon <-rep(NA, N) #iterate over all patchesfor (i in1: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 graph.i <-graph_from_adjacency_matrix(prob.matrix.i, mode="undirected", weighted=TRUE)#create all shortest paths between patches pstar.mat.i <-distances(graph.i, weights=-log(E(graph.i)$weight))# back-transform to probabilities of connectedness pstar.mat.i <-exp(-pstar.mat.i)#multiply all areas and movement probabilities together productMat <-outer(area.i, area.i)*pstar.mat.i#final calculation con.i <-sqrt(sum(productMat))/landarea*100 diffCon[i] <- (baseline-con.i)/baseline*100 }return(diffCon)}dPC_Bool <-diff_connectivity(prob.matrix=weightedMatBool,area=areaBool,landarea=AL,baseline=pcMod)spNodesBool=vect(sitesBool)dataDPCBool=data.frame(site=1:nrow(patchesBool),coords=crds(spNodesBool),dPC=dPC_Bool,areaBool=areaBool)
Now that we have our type one layer inputs, we need to assign nesting and foraging weights to each cover type that reflect their suitability for the focal generic bird species.
#rasterise the motorway again but to use as habitat (setting habitat suitability to zero)m62Hab<-rasterize(m62,field=0,classif,background=1)#remove motorway from urban foraging potentialurbanMean<-urbanMean*m62Hab
Now we create a stacked raster with all layers combined.
By multiplying all layers by their respective nesting suitability values we can achieve a baseline multivariate habitat estimate for each cell in the landscape.
#sum all layers weighted by nesting suitabilitynestHab<-(woodMean*woodWeightNesting)+(grassMean*grassWeightNesting)+(shrubMean*shrubWeightNesting)+(waterMean*waterWeightNesting)+(urbanMean*urbanWeightNesting)
Next we create a binary layer delimiting all cells greater than or equal to 0.5 as habitat patches. This is so a) we can subset our data a little to exclude very small patches and b) we can limit the spatial extent of our patches to this original nesting “baseline” throughout all steps.
#init binary mask of habitat patches with values >0.5 to create polygons and compute areas binaryPatchMask= nestHab binaryPatchMask[binaryPatchMask>=0.5]=1 binaryPatchMask[binaryPatchMask!=1]=NA#combine adjacent habitat cells nestPatches<-patches(binaryPatchMask,directions=8)#polygonize patches for cost distance analysis nestPatches<-as.polygons(nestPatches,dissolve=TRUE) nestPatches$Area<-expanse(nestPatches)#subset patches >= 0.5 ha nestPatches<-nestPatches[nestPatches$Area>=5000]
Now we have some initial patches to work with, we will consider neighbourhood effects. As seen in the lecture, we can increase or decrease the habitat amount (~habitat quality) according to the influence of the surrounding matrix, where land uses with negative neighbourhood effects (“edge effects” in common parlance) reduce quality and land covers that offer complementary resources (i.e. foraging potential) have a positive influence on habitat quality.
We the patches tentatively defined above, the next step is to create everyone’s favourite bit of spatial ecology - the neighbourhood weights matrix.
The first of these will have the goal of parameterising negative edge effects which we will consider to stem from the urban and grassland layers (following the approach taken in the reading for this week and making the simplifying assumption that all grassland in the study area is improved agricultural grassland).
#simulate neighbourhood effects.#urban edge effect radius (in metres - from Eycott et al. 2011) radiusUrban<-76#neighbourhood weights matrix#get number of pixels required in each direction to cover 76 metres nPixUrban<-round(radiusUrban/res(classif)[1])#add one to ensure matrix is has odd number of columns and rows (i.e. so the focal cell is bang in the middle) nPixUrban<-(nPixUrban*2)+1#now create the matrix matDistUrban<-matrix(1:nPixUrban^2,nrow=nPixUrban,ncol=nPixUrban)#get x and y position of focal cell x<-ceiling(ncol(matDistUrban)/2) y<-ceiling(nrow(matDistUrban)/2)#get focal cell from x and y focalCellUrban<-matDistUrban[x,y]#identify focal cell and its index in the matrix (cell number) indFocalUrban<-which(matDistUrban==focalCellUrban,arr.ind = T)#creat a list to hold distance values distListUrban<-list()#compute distancesfor(i in1:nPixUrban^2){ ind.i=which(matDistUrban==i,arr.ind=T) diffX<-abs(ind.i[1,1]-indFocalUrban[1,1])*res(classif)[1] diffY<-abs(ind.i[1,2]-indFocalUrban[1,2])*res(classif)[1] dist.i<-sqrt(diffX^2+diffY^2) distListUrban[[i]]=dist.i }#add distance values matDistUrban[]<-unlist(distListUrban)#set cells outside search radius to zero matDistUrban[matDistUrban>radiusUrban]=NA#calculate edge effect decay rate based on Eycott et al. #get constant for kernel shape based on max 76 metre edge effect (for urban). We set 1/1000000 just to get a very small value (but >0) fo edge effects at 76 metres. alphaUrban=log(1/1000000)/76#compute cell edge effect values matKernelUrban<-exp(alphaUrban * matDistUrban)#set focal cell to zero - because we are only interested in the effects of neighbouring cells matKernelUrban[x,y]<-0
With the neighbourhood weights set up, we can pass these to the focal() function to compute the edge effect values applied to each cell in the landscape. We will use the extend() function to increas the size of our raster at the edges by 10 pixels and giving those pixel a value of zero - otherwise the focal function will remove cells near the edge as they would be involve in calculations with NA cells (because the 76 metre neighbourhood would extend beyond the study area boundary).
#sum edge effects from all surrounding cells focalUrban<-focal(extend(urbanMean,c(10,10),fill=0),matKernelUrban,fun=sum)#crop to orginal extent focalUrban<-crop(focalUrban,ext(classif))plot(focalUrban)
Now just repeat the same process but for grassland for which we apply a 30 metre edge effect according to the same literature as above.
############repeat for grassland cover#set up radius and equivalent number of pixels radiusGrass<-30 nPixGrass<-round(radiusGrass/res(classif)[1]) nPixGrass<-(nPixGrass*2)+1#create distance/neighbourhood matrix distMatGrass<-matrix(1:nPixGrass^2,nrow=nPixGrass,ncol=nPixGrass)#get row and column number of the focal cell x<-ceiling(ncol(distMatGrass)/2) y<-ceiling(nrow(distMatGrass)/2)#identify focal cell and its index in the matrix (cell number) focalCellGrass<-distMatGrass[x,y] indFocalGrass<-which(distMatGrass==focalCellGrass,arr.ind = T)#init list for distances distListGrass<-list()for(i in1:nPixGrass^2){ ind.i=which(distMatGrass==i,arr.ind=T) diffX<-abs(ind.i[1,1]-indFocalGrass[1,1])*res(classif)[1] diffY<-abs(ind.i[1,2]-indFocalGrass[1,2])*res(classif)[1] dist.i<-sqrt(diffX^2+diffY^2) distListGrass[[i]]<-dist.i }#add distances to the matrix distMatGrass[]=unlist(distListGrass) distMatGrass[distMatGrass>radiusGrass]=NA#plot(rast(distMatG))#res(wood.i)[1]^2#get constant for kernel shape based on max 30 metre edge effect (for urban). We set 1/1000000 just to get a very small value (but >0) fo edge effects at 30 metres. alphaGrass=log(1/1000000)/30#compute cell edge effect values matKernelGrass<-exp(alphaGrass * distMatGrass)#set focal cell to zero matKernelGrass[x,y]<-0#compute focal layer focalGrass<-focal(extend(grassMean,c(10,10),fill=0),matKernelGrass,fun=sum)#crop to original extent focalGrass<-crop(focalGrass,ext(classif))plot(focalGrass)
With the negative neighbourhood effects done for habitat suitability, we can now focus on setting our cost values based on type 1 membership. We do this in the same way as for computing nesting suitability values - weighting each cell as a function of its land cover type and the cost assosiated with that type.
############################Create cost layers#assign cost weights to other cover types woodCostWeighted<-woodMean grassCostWeighted<-grassMean*10 urbanCostWeight<-urbanMean*5 waterCostWeight<-waterMean*10 shrubCostWeight<-shrubMean*2
Now add these together for a total weighted cost per cell.
To this initial cost calculation we can add in neighbourhood effects as for habitat delineation. Here we just use the neighbourhood effect surfaces created above, weight these by the corresponding cost values and add to the mix.
We should also consider the motorway again though, so create a separate layer for this (assuming the same edge effect distance for urban but applying our much higher cost value that we set up earlier).
#get contribution of motorway to functional cost focalM62Cost<-focal(extend(m62Cost,c(30,30),fill=0),matKernelUrban,fun="sum")#crop focalM62Cost<-crop(focalM62Cost,woodMean)#weight according to cost focalM62Cost<-focalM62Cost*100#add in neighbourhood effects on cost and the motorway cover type1CostEdge<-type1Cost+costUrbanFocal+costGrassFocal+focalM62Cost
Foraging potential
The next step in the functional habitat delineation process is to assign foraging potential values for all cells and add these to our original habitat suitability based on nesting potential.
First, create a layer to present the total foraging potetial of each cell by summing the type one menbership values of all land cover class weighted by the foraging values set up at the start of this practical.
#sum all foraging potential valuesforageWeights=woodMean*woodWeightForaging+grassMean*grassWeightForaging+shrubMean*shrubWeightForaging+waterMean*waterWeightForaging+urbanMean*urbanWeightForaging
The next step is to assess the accessibility of these foraging resources around each cell. This requires a little more thought than the computation of the negative neighbourhood effects carried out earlier. This is because we have to consider organismal movement as part of the process. In other words, whereas the negative neighbourhood effects can be considered as a function of Euclidean distance, in the case of positive neighbourhood effects, the ability to access such resources depends on the cost of moving through the immediate landscape. So, we have to build our neighbourhood weights matrix on cost distances rather than straight-line distance. That means we need to build a transition layer based on our newly created functional cost layer.
So we create a new transition layer based on the more functional derivation created in the previous step (where costs are now a weighted some of type 1 membership values).
Now, to get the cost of moving out from each of our habitat patches (which are currently contained in the binaryPatchMask object created earlier) towards neighbouring resources, we compute the cumulative cost of moving out from each habitat cell.
#get cumulative cost of moving from each patch out into the matrix cumulativeCost=accCost(trFunc,crds(binaryPatchMask))*9#assign a crs to this layercrs(cumulativeCost)="EPSG:27700"# this sets shape of neg. exp. kernel to compute movement probability distForage =-log(0.05) /250# this creates the raster kernel#compute movement probabilities forageRast <-exp(-distForage * cumulativeCost)#set cells with a value of 1 to zero (to only consider neighbouring cells when calculating positive neighbourhood effects) forageRast[forageRast==1]<-0
To complete our calculation of positive neighbourhood effects (foraging potential) we need to multiply our layer of movement probability by the foraging potential layer, then pass this to the focal() function to sum all values around each cell.
#multiply the movement probability layer by the foraging weights layer forage<-(rast(forageRast)*forageWeights)#create a neighbourhood weights matrix to pass to the focal() function winForage=matrix(1,ncol=round(250/res(classif)[1])+1,nrow=round(250/res(classif)[1])+1)
In the next step we will ensure that the weights matrix for positive neighbourhood effects sums to one so that positive neighbourhood effects are returned as a probability (range from 0-1). This avoids having very large numbers that would give unrealistically high habitat values (ultimately we want the final habitat suitability score for each cell to be between 0 and 1). Also, we can’t assume that all resources are equally and continuously accessible, so truncating positive neighbourhood effects as 0-1 gives a more conservative estimate and gives us an idea of the potential foraging suitability.
#make sure the weights matrix sums to one winForage=winForage/sum(winForage)#pass all this to the focal() function to caluclate sum of all foraging potential in the neighbourhood focalForageSum<-focal(extend(forage,10,fill=0),w=winForage,fun="sum")#crop toorginal extent focalForageSum<-crop(focalForageSum,woodMean)#keep only habitat cells by applying the binary mask forage<-focalForageSum*binaryPatchMask
Finally, get the habitat suitability scores for each cell by add positive and removing negative neighbourhood effects.
#get final habitat values by remove negative and adding positive neighbourhood effects funcHab<-nestHab-focalUrban-focalGrass+forage
Next, in order to work with this layer to use in a connectivity analysis, we need to delineate patches as usual.
#any values that have crept >1 need to be set to max of 1 funcHab[funcHab>1]<-1#use a 0.5 threshold to delineate habitat cells funcHab[funcHab<0.5]<-0#get patches of adjacent cells based on the final functional habitat amount patchesFunc<-patches(funcHab,direction=8,zeroAsNA=TRUE)#convert to polygons patchesFunc<-as.polygons(patchesFunc,dissolve=TRUE)#subset patches >=0.5 ha patchesFunc$Area<-expanse(patchesFunc) patchesFunc<-patchesFunc[patchesFunc$Area>=5000,]
As per the graph theoretic approach to conectivity, we need to compute patch area for the equation. In the case of functional habitat however, we want to use the “amount” of habitat (the sum of all suitability scores within the patch) rather than the area of the contiguous area of the patch. We do this just by using the extract function to sum all cell values within the patch.
#get habitat amount based on functional habitat values per cell habAmount<-extract(funcHab,patchesFunc,method="simple",fun="sum")#"amount" = cell suitability values multiplied by the cell area habAmount<-habAmount[2]*(res(classif)[2]^2)#add to the patch polygon object for use later on. patchesFunc$habAmount<-habAmount
With this set up, we can proceed to computing the cost distance matrix and habitat connectivity metric.
#get patch centroids for least cost path delineation sitesFunc<-SpatialPoints(crds(centroids(patchesFunc)))#create cost distance matrix costMatFunc <-as.matrix(gdistance::costDistance(trFunc,sitesFunc)) costMatFunc=costMatFunc*res(classif)[1]# set alpha to determine colonization probability of the species (max dispersal capacity as 10km after Dennis et al. 2024) alphaDisp=log(0.05) /10000# negative exponential of colonization kernel x distance get probability weightedMatFunc <-exp(alphaDisp * costMatFunc)# set diag to zero (so we we don't count movement from patches to themselves)diag(weightedMatFunc) <-0# create a graph from this weighted adjacency matrix weightedGraphFunc <- igraph::graph_from_adjacency_matrix(weightedMatFunc, mode="undirected", weighted=T)## calculate Connectivity metric (RFH)# calculate all shortest paths between nodes pstar.matFunc <-distances(weightedGraphFunc, weights=-log(igraph::E(weightedGraphFunc)$weight))# back-transform to probabilities of connectedness pstar.matFunc <-exp(-pstar.matFunc)# get study area in m2 AL <-ncell(classif)*res(classif)[1]^2# get product of all patch areas ij and multiply by probabilities above rhMat <-outer(patchesFunc$habAmount,patchesFunc$habAmount) * pstar.matFunc rhMatSum <-sum(rhMat)# divide by total area of the study squared to get the PC metric RFH <- (sqrt(rhMatSum) /as.numeric(AL))*100#compute contribution of each patch to RFH weighted_dRH <-diff_connectivity(prob.matrix=weightedMatFunc,area=patchesFunc$habAmount,landarea=AL,baseline = RFH)dataDRH=data.frame(site=1:nrow(patchesFunc),areaBool=expanse(patchesFunc),dRH=weighted_dRH,habAmount=patchesFunc$habAmount)spNodesFunc=vect(sitesFunc)
Next we do some plotting using the RClolourBrewer library for help with colour ramps.
