Week Seven: Reachable Functional Habitat

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 functional approach that builds on the species-oriented aspects of the previous analyses we have carried out. 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 the only new package this week is RColorBrewer so you may just choose to load just that one.

Set the working directory and install packages

setwd('./Week7')

#install necessary packages

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

# load libraries for spatial data
#NB raster library for use with gdistance
library(raster)

Loading required package: sp

library(terra)

terra 1.9.1

Load in a classified land cover layer of the study area created from a Random Forest model. This layer contains five classes where 1 = woodland, 2 = grassland, 3 = shrub, 4 = water and 5 = urban.

#actual classified land cover layer
classif=rast("classif_RF.tif")

#plot

plot(classif)

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 class
woodland=classif

#set all values that are not 1 to zero
woodland[woodland!=1]=0

# use the patches() function to
woodland=patches(woodland,directions=8,zeroAsNA=T)


#polygonize patches for cost distance analysis
patchesBool=as.polygons(woodland,dissolve=TRUE)
patchesBool$Area=expanse(patchesBool)

#subset patches >= 0.5 ha
patchesBool=patchesBool[patchesBool$Area>=5000]

Next we assign costs to each class in the land cover map.

#set the costs for each land cover
landCosts=c(1,10,2,10,5)

#create a Boolean layer of landscape cost
booleanCost=classify(classif,cbind(c(1,2,3,4,5),landCosts))

#plot
plot(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 magnitude higher than the simple “urban” class.

#get mask of motorway (high cost)
m62=vect("M62.shp")

#rasterise m62 to combine with the cost analysis
m62Cost=rasterize(m62,field=1,classif,background=0)

Now set all cells in the booleanCost object that are not zero in the m62 layer to 100 to assign very high cost.

#set the motorway cost to 100
booleanCost[m62Cost!=0]=100

Now we are ready to start modelling so first we build a transition layer as in previous weeks.

#load gdistance for cost analysis
library(gdistance)

Loading required package: igraph

Attaching package: 'igraph'

The following objects are masked from 'package:terra':

    blocks, compare, union

The following object is masked from 'package:raster':

    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

#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 delineation
sitesBool=SpatialPoints(crds(centroids(patchesBool)))

#create cost distance matrix
costMatBool = as.matrix(costDistance(trBool,sitesBool))

#multiply by the resolution of the raster
costMatBool=costMatBool*res(classif)[1]


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

Now build the adjacency matrix and associated graph for our “Boolean” patches.

library(igraph)

# negative exponential of colonization kernel x distance get probability
weightedMatBool = 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 matrix
weightedGraphBool = igraph::graph_from_adjacency_matrix(weightedMatBool, mode="undirected", weighted=T)

Compute shortest movement proababilities between patches and combine with patch area to calculate the PC metric.

## calculate Connectivity metric (PC)

# calculate all shortest paths between nodes
pstar.matBool = distances(weightedGraphBool, weights= -log(igraph::E(weightedGraphBool)$weight))

# back-transform to probabilities of connectedness
pstar.matBool = exp(-pstar.matBool)

# get study area in m2
AL = ncell(classif)*res(classif)[1]^2

# get patch areas
areaBool = patchesBool$Area

# get product of all patch areas ij and multiply by probabilities above
pcMat = outer(areaBool,areaBool) * pstar.matBool

pcMatSum = sum(pcMat)

# divide by total area of the study squared to get the PC metric
pcMod = (sqrt(pcMatSum) / as.numeric(AL))*100

print(pcMod)

[1] 5.224294

Next we get the contribution of each patch to connectivity using the same function as we built in the last practical.

#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 an empty vector "diffCon" 

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



dataDPCBool=data.frame(site=1:nrow(patchesBool),dPC=dPC_Bool,
                       areaBool=areaBool)

#add the key result to the patches object.
patchesBool$dPC=dataDPCBool$dPC

Next we do some plotting using the RClolourBrewer library for help with colour ramps.

library(RColorBrewer)

#get colour palette to use
cols = brewer.pal(7, "Greens")

#build function for setting colour ramp
pal = colorRampPalette(cols)


#plot the patches according to their contribution to connectivity
plot(patchesBool, "dPC", type="continuous", col=pal(100))

Next we will work with the type 1 membership values to create a multivariate habitat model.

Load in layers containing type 1 membership values obtained from a Random Forest classifier

#woodland type-1 membership
woodMean=rast("woodMeanClipped.tif")

#grassland membership
grassMean=rast("grassMeanClipped.tif")

#shrub layer membership
shrubMean=rast("shrubMeanClipped.tif")

#water class membership
waterMean=rast("waterMeanClipped.tif")

#urban cover membership
urbanMean=rast("urbanMeanClipped.tif")

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 a focal generic species. These values are taken from the supplementary data of a paper by Emma Gardner and colleagues based on a specialist woodland bird focal generic species available here

#############################Woodland weights


woodWeightNesting=0.78125
woodWeightForaging = 0.826


##################################################################
######### Grassland weights


grassWeightNesting=0
grassWeightForaging=0.001


#################################################
############shrub layer weights

shrubWeightNesting=0.43333
shrubWeightForaging=0.792



#######################################
#############Water weights

waterWeightNesting=0
waterWeightForaging=0.001



#####################################
#########################Urban weights

urbanWeightNesting=0.21111
urbanWeightForaging=0.187

As in the case of cost, we will give special attention to the M62 - setting this to zero in terms of habitat suitability.

#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 nesting and foraging potential
urbanHab=urbanMean*m62Hab

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 suitability
nestHab=(woodMean*woodWeightNesting)+(grassMean*grassWeightNesting)+(shrubMean*shrubWeightNesting)+(waterMean*waterWeightNesting)+(urbanHab*urbanWeightNesting)

#take a look at this initial layer

plot(nestHab)

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 distances
  for(i in 1: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 increase 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 in 1: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
  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.

  type1Cost=woodCostWeighted+grassCostWeighted+urbanCostWeight+waterCostWeight+shrubCostWeight

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 total cost.

  costUrbanFocal=focalUrban*5
  costGrassFocal=focalGrass*10

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(10,10),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 potential of each cell by summing the type one membership values of all land cover classes weighted by the foraging values set up at the start of this practical.

#sum all foraging potential values
forageWeights=woodMean*woodWeightForaging+grassMean*grassWeightForaging+shrubMean*shrubWeightForaging+waterMean*waterWeightForaging+urbanHab*urbanWeightForaging

To assess the accessibility of these foraging resources around each cell 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).

#
  trFunc=transition(raster(type1CostEdge),transitionFunction = function(x) 1/mean(x),directions = 4)

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 (by passing habitat cell coordinates to the accCost() function).

      #get cumulative cost of moving from each patch out into the matrix
      cumulativeCost=accCost(trFunc,crds(binaryPatchMask))*res(classif)[1]
      
      #assign a crs to this layer
      crs(cumulativeCost)="EPSG:27700"
      
      # this sets the shape of neg. exp. kernel to compute movement probability according to a 250 metre froaging distance from Gardner et al. 
      distForage = -log(0.05) / 250
     
    
      #compute movement probabilities
      forageProb = exp(-distForage * cumulativeCost)
      
      #set cells with a value of 1 to zero (to only consider neighbouring cells when calculating positive neighbourhood effects)
      forageProb[forageProb==1]=0
      
      #inspect
      plot(forageProb)

To complete our calculation of positive neighbourhood effects (foraging potential) we need to multiply our layer of movement probabilities 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(forageProb)*forageWeights)
      plot(forage)

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

With the neighbourhood weights set we will now ensure that the weights matrix for positive effects sums to one so that positive 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 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 calculate sum of cell foraging     potential values in the neighbourhood
      focalForageSum=focal(extend(forage,10,fill=0),w=winForage,fun="sum")
      
      #crop to orginal extent
      focalForageSum=crop(focalForageSum,woodMean)
     

  #keep only habitat cells by applying the binary mask
  forageFinal=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+forageFinal

  #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

  plot(funcHab)

Next, in order to work with this layer to use in a connectivity analysis, we need to delineate patches as usual.

  #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,]

#plot the patches 
  
 plot(patchesFunc)

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
  
 #print result
  print(RFH)

[1] 2.488799

#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=patchesFunc$Area,
                   dRH=weighted_dRH,habAmount=patchesFunc$habAmount)
   

#add dRH values to the functional patches object
patchesFunc$dRH=dataDRH$dRH

Do some plotting using the RClolourBrewer libarary again.

#set colour palette
cols = brewer.pal(7, "Greens")

#make function to set colour ramp for plotting
pal = colorRampPalette(cols)

Finally plot the Boolean and Functional habitat patches according to their respective estimates of patch contribution to connectivity.

#plot dPC
plot(patchesBool, "dPC", type="continuous", col=pal(100))

#plot dRH
plot(patchesFunc, "dRH", type="continuous", col=pal(100))

Clearly the two approaches place different emphases on the influence of the surrounding matrix with changes in the relative contribution of each patch resulting primarily from a combination of patch shape and negative neighbourhood effects.

Things to try next

Now that you have reached the end of this week’s practical, you may wish to explore further options. An obvious limitation of the approach we have taken here is that, working with type 1 membership values, we are constrained to one realisation of each cell’s membership to each land cover type. It also means that we are stuck with one single, crisp boundary for our habitat patches, which may still be unrealistic. However, we know that the membership values (to woodland, grassland, urban etc.) also comes with some uncertainty. To see an example of how we might leverage this uncertainty to build many realisations of the study area landscape and avoid having to resort to a hard boundary for patch delineation, take a look at the following github repository:

shrodinger repo

This repo contains dummy data and code that will run a Monte Carlo analysis to obtain a distribution of functional habitat and connectivity values for two simple patches as reproduced in this paper:

Dennis, M., Huck, J., Holt, C., McHenry, E., Andersson, E., Sharma, S. and Haase, D., 2026. In search of Schrödinger’s patch: a functional approach to habitat delineation. Landscape Ecology - link: here