GEOG71922: Week Four: Landscape pattern and process

In this practical you will be guided through some spatial techniques for exploring basic concepts in landscape ecology. So far on the course we have been mainly focused on species ecology in that we have considered point data and species-environment responses as the basis of our analysis. There are some broader ideas and approaches derived from the discipline of landscape ecology that are very pertinent and useful from a spatial ecology point of view. Although distributuion data are useful and versatile, our analyses necessarily depend on the context of our study, that being the underlying landscape. So far we have only looked at landscape characteritics as a function of land-use-land-cover but there are many other ways in which landscape qualities influence ecological processes. Typical qualities that concern ecologists are those related to fragmentation/aggregation, connectivity/isolation and richness (diversity).

In this week’s practical we will look at different ways of characterising landscapes and take an initial look at connectivity based on the use of least cost paths.

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/WeekFour")

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

#Install packages

install.packages(c("terra","landscapemetrics","gdistance","ggplot2","patchwork","raster"))

#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

# 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

land_cover$landcover<- c("built","water", "ground veg","shrub layer","canopy")

#tell R that these are what we want to name each level by passing two column data frame of the classes (column 1) and their names (column 3)

levels(GM) <- land_cover[,c(1,3)]

#create a list of colours for plotting

land_col <- c("white","blue","green","yellow","darkgreen")

Plot the land covers with a new legend

plot(GM, legend = T, col = land_col)

For the first part of the practical we will work only with the canopy layer of the GM object. So the next step is to create a new object that describes only the canopy cover in the study area:

# 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 binary raster where trees are 1 and all other values 0.
gmW <- c(NA,NA,NA,NA,1)

# create a matrix for reclassification with the IS
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

#woodland binary layer from re-classification matrix

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

Now, with this new canopy layer we can carry out some basic analysis using some common landscape metrics using the landscapemetrics package. The first kind of metric we will use are called Patch-level metrics. Patch-level means that these metrics are meaningful at the level of individual patches (i.e. stands of trees in this case). Such metrics include mean patch size, core area (the area of a single patch minus the edge - where the “edge” defaults to as one pixel width) and shape (the complexity of patch shape can influence ecologically relevant variables such as the length of edge habitat in the landscape). See the link to the fragstat documentation on Canvas for more details on how these metrics area calculated.

We will calculate patch area, core area, shape index and contiguity index. The contiguity index measure the extent to which land cover cells form a single contigous patch where a value of 0 reflects a one-pixel patch and a value of 1 reflects a situation where all cells form a single patch.

############################################
#patch-level analysis for woodland
############################################

#load library for computing landscape metrics
library(landscapemetrics)

#Calculate patch-level metrics.

#patch area
area_p=lsm_p_area(gmWood)

#inspect the first 10 results
print(area_p,n=10)

# A tibble: 27,400 × 6
   layer level class    id metric value
   <int> <chr> <int> <int> <chr>  <dbl>
 1     1 patch     1     1 area    0.36
 2     1 patch     1     2 area    0.72
 3     1 patch     1     3 area    0.09
 4     1 patch     1     4 area    0.09
 5     1 patch     1     5 area    0.18
 6     1 patch     1     6 area    0.72
 7     1 patch     1     7 area    0.36
 8     1 patch     1     8 area    0.18
 9     1 patch     1     9 area   28.4 
10     1 patch     1    10 area    0.27
# ℹ 27,390 more rows

We can look at the distribution of patch area using the hist() function. Note - where results are returned as areas, the landscapemetrics package gives this in Hectares. For our study area, the vast majority of patches and < 1 Ha so we can subset values below this when producing the histogram.

hist(area_p[area_p$value<1,]$value)

Now compute core area (area minus the edge - set to one pixel width).

#Core Area Index
ca_p=lsm_p_cai(gmWood,edge_depth = 1)

Contiguity Index

#Contiguity
contig_p=lsm_p_contig(gmWood)

Shape Index

#Shape
shape_p=lsm_p_shape(gmWood)

Now let’s put these metrics together in a data frame and plot their associated values.

#build data frame of patch-level metrics
patchStats=data.frame(area=area_p$value,coreArea=ca_p$value,
          shapeComplexity=shape_p$value,contiguity=contig_p$value)


#use the pairs() function to plot all pairs of variables in a data frame
pairs(patchStats)

Plotting the results in this way shows that, apart from the clear relationships where shape complexity increases with patch size, these patch-level metric appear to exhibit low redundancy (the measure relatively distinct patch attributes).

Now we will move on to some class level metrics to see if these shed light on how attributes vary between land cover classes. We will compute the number of patches per class, the proportion cover by each class, the interspersion and juxtaposition index, the contiguity index and the clumpiness index

#class level metrics

#proportion cover
pland_c=lsm_c_pland(GM)

#mean patch area
area_c=lsm_c_area_mn(GM)

#number of patches
np_c=lsm_c_np(GM)

#interspersion and juxtaposition
iji=lsm_c_iji(GM)


#Clumpiness
clump=lsm_c_clumpy(GM)


#Euclidean Nearest Neighbour
enn_c=lsm_c_enn_mn(GM)

Again, we can collect these metrics in a data frame and then plot these class-level metrics as bar plots using ggplot, and arrange them as a grid using the patchwork package.

#Create data frame of class-level metrics
class_stats=data.frame(pland=pland_c$value,nPatches=np_c$value,interspersion=iji$value,clumpiness=clump$value,area=area_c$value,ENN=enn_c$value)

class_stats$cover=c("built","water","grass","shrub","canopy")

Build separate plots for each metric using ggplot(). Note we add “theme(axis.text.x = element_text(angle = 45))” to set the x-axis labels at an angle (because they overlapped when I first ran this)

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.5.2

#Pland
propPlot=ggplot(class_stats,aes(x=cover,y=pland))+geom_bar(stat="identity",fill="white",color="black")+
theme_minimal()+theme(axis.text.x = element_text(angle = 45))

#Mean patch area
areaPlot=ggplot(class_stats,aes(x=cover,y=area))+geom_bar(stat="identity",fill="white",color="black")+
  theme_minimal()+theme(axis.text.x = element_text(angle = 45))

#Number of patches
patchPlot=ggplot(class_stats,aes(x=cover,y=nPatches))+geom_bar(stat="identity",fill="white",color="black")+
theme_minimal()+theme(axis.text.x = element_text(angle = 45))

#interspersion and juxtaposition
ijiPlot=ggplot(class_stats,aes(x=cover,y=interspersion))+geom_bar(stat="identity",fill="white",color="black")+
  theme_minimal()+theme(axis.text.x = element_text(angle = 45))

#Clumpiness
clumpPlot=ggplot(class_stats,aes(x=cover,y=clumpiness))+geom_bar(stat="identity",fill="white",color="black")+
  theme_minimal()+theme(axis.text.x = element_text(angle = 45))

#Mean nearest neighbour
ennPlot=ggplot(class_stats,aes(x=cover,y=ENN))+geom_bar(stat="identity",fill="white",color="black")+
  theme_minimal()

Plot them together using the patchwork package

library(patchwork)

Attaching package: 'patchwork'

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

    area

propPlot+areaPlot+patchPlot+ijiPlot+clumpPlot+ennPlot

Next we will look at diversity with respect to land cover in Greater Manchester. An initial, very simple metric is just to count the number of land cover classes. This is of course arbitrary because I decided how many there would be when I first created the classification, but just for completeness:

#----------------------------------#
#diversity-related metrics
#----------------------------------#


#richness as just the number of unique values in the GM raster
richness <- nrow(unique(GM))

#print
print(richness)

[1] 5

Richness, though a useful description and intuitive is generally considered a poor measure of diversity as it does not take into consideration the proportion of each class in the dataset. One example of an index that does is the Shannon’s diversity index (or Shannon’s H). We will look at how this is calculated below (but see the suggested reading for other indices worth exploring).

The basic approach used by the Shannon index is to take each class and calculate the proportion of the landscape/sampling site/quadrat etc. attributed to it, then multiply this number by its logarithm and sum the result for all values.

Where:

is the proportion of each class and is the natural logarithm of the proportion of each class.

We can compute this easily in R as follows(note the minus sign on the last step to return positive values - this is just a “mathematical convenience” to make the result more intuitive).

#Shannon's diversity index 


#Get number of all values (classes) in the land cover map
C <- table(values(GM))

#get this as a proportion
P <- C / sum(C)

#Shannon Index is just the sum of these proportions multiplied by their logarithm
H <- -sum(P * log(P))

print(H)

[1] 1.380972

This number is the Shannon’s Index for the whole of Greater Manchester which isn’t particularly informative (unless you are very interested in comparing this value to that of another city region for some large-scale investigation). To get a better idea of how land cover diversity varies across Greater Manchester we need to use a “moving window” type of analysis, which we do using the focal() function that we explored last week.

So let’s look at calculating Shannon’s Index for a neighbourhood around each cell. Such a calculation might be useful when considering the effects of landcover diversity at a local scale (for example on species richness, resilience of vegetation to drought or flooding in urban areas).

Remember that the basic premise is to start from a position of considering the proportion of each landcoer. Conceptually then, what we need to do is similar to the above calculation for Shannon’s H. We just need to calculate this proportion for each landcover class within a set neighbourhood (a search radius around each cell) and then multiply by the logarithm of these proportions and add them together. Note that in the code below we create a “blank” raster layer (with values of 0) that we then use as the basis to sum the products of our calculations (by adding each layer to this blank raster in turn). Then, don’t forget to multiply by minus one so that results are positive.

There is one difference here compared to who we carried out neighbourhood analysis in Week 3. In that case we built our spatial weight matrix from scratch (because I wanted to teach you how to do this). We can build one with less code however by using the focalMat() function built into the terra package. This function just takes the raster you are working with and a radius to set the neighbourhood. The default is to return a spatial weights matrix that sums to one (just like we did in last week’s practical). As we are building a relatively complicated function below, we’ll spare ourselves having to build the spatial weighted matrix, and just use the focalMat() function.

#diversity moving window function
# this function will take two arguments in order to run: a landcover raster and a radius (a number)

diversity <- function(landcover, radius) {
  # create an object that signifies the number of unique layers in the landcover raster (i.e. the number 5 in our example)
  n=unique(values(landcover)) 
  #remove the NA class
  n=n[!is.na(n)]
  
  #Create focal weights matrix
  fw.i = focalMat(landcover, radius, type="rectangle")

  
  #Compute focal means of indicators
  D = landcover # blank raster to hold all (sum) values
  values(D) = 0
  
  #for each landcover, create a moving window map and make calculations
  for (i in n) {
   
    # raster for summing cells of landcover class i in the neighbourhood of each cell. use na.rm=T to ignore NA        cells
    focal.i <- focal(landcover == i, fw.i,na.rm=TRUE) 

    #take log(p)*p of each proportion value 
    logP=focal.i*log(focal.i)
     #ensure NAs are zero so we can sum the results
    logP[is.na(logP)]=0
    
    #and add to the blank raster
    D <- D + logP 
    
  }
  
  D <- D*-1#take negative
  
  D[D==0]=NA #add NAs back in
  return(D)
}

Now run this function for a radius of 100 metres and plot the result.

#calculate diversity moving window
diversity.100m <- diversity(landcover=GM, radius=100)
plot(diversity.100m)

Part Two: landscape resistance and least cost path

In Part Two of this practical we will look at the land cover within Greater Manchester from the perspective of movement cost also known as “landscape resistance”. We will 1. reclassify our land cover raster according to assumed movement cost for a generic focal woodland species. This is a hypothetical species that captures typical characteristics of a species group such as nesting requirements, foraging preferences and sensitivity to land cover (that can be used to parameterise thinsg such as movement cost or edge effects). We will be using values taken from Eycott et al. which were produced through a Delphi process (expert consultation). That paper can be found here: Eycott et al. 2011 and an example of these being used can be found here: Dennis et al. 2024.

First we creat two numerical vectors, one for the classes in the GM object and one for the corresponding cost values from Eycott et al.

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

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

lcm_classes=1:21

#costs from Eycott et al.
costs=c(1,3,rep(10,2),rep(4.35,3),2.5,rep(2.22,2),2.44,5,20,10,14.29,rep(20,4),rep(5,2))

#combine
costMat=cbind(lcm_classes,costs)

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

#inspect the result
plot(costGM)

The next thing we need to do is to create some destinations between which we will compute the least cost path. We will focus on the woodland class (so we can use the Eycott values for a focal generic woodland species). For this analysis we need two xy coordinates, so the most straightforward way to get this is to convert the woodland class in the GM object to polygons and then, choosing two patches of woodland, get the coordinates of the patch centroids and calculate the least cost path between the two.

The first step is to use the patches() function to assign all cells in the gmWood object we created earlier to individual patches. We will use an 4-neighbour rule for this to select more contiguous (compact) patches.

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 largest patches so let’s subset patches of woodland over 100 hectares.

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

#Subset to woodland patches > 100 hectares
polyWood=polyWood[polyWood$area>1000000,]

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

Now we have the coordinates, separate X and Y so we can pin point two patches at either end of the study area.

#Get the X and Y coordinates
polyWood$X=woodCentres[,1]
polyWood$Y=woodCentres[,2]


#Choose the patch further to the NW of the study area - this is called Bottling Wood.
Bottling=centroids(polyWood[which.max(polyWood$Y-polyWood$X),])

#Get the patch furthest to the East (this is in Etherow country park)

Etherow=centroids(polyWood[which.max(polyWood$X),])

#take a look at where these are.
plot(polyWood)

plot(Bottling,add=T,cex=1.5,col="red")

plot(Etherow,add=T,cex=1.5,col="red")

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) wheras we want the inverse of this (cost of movement i.e. distance).

library(gdistance)

Loading required package: raster

Loading required package: sp

Loading required package: igraph

Attaching package: 'igraph'

The following object is masked _by_ '.GlobalEnv':

    diversity

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

    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(raster)
#Create transition layer.
land_cond <- transition(raster(costGM), transitionFunction=function(x) 1/mean(x), 8)

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

lcp <- shortestPath(land_cond,origin=crds(Bottling),goal=crds(Etherow),output="SpatialLines")

We can plot this over our land cover object

#land cover
plot(GM)

#Bottling Wood
plot(Bottling,add=T,cex=1.5,col="red")

#Etherow patch
plot(Etherow,add=T,cex=1.5,col="red")

#Least cost path
plot(lcp,add=T,col="red",lwd=2)

The least cost path is informative but not particularly realistic so we might instead choose to create a broader least cost corridor instead. This is easily done by creating a layer of the cumulative cost moving out from each point and then subsetting the locations (cells) below a certain percentile of the summed cost values.

#get cumulative costs from each centroid
Bottling.cost <- accCost(land_cond, crds(Bottling))

Etherow.cost <- accCost(land_cond, crds(Etherow))


#get least-cost corridor - use the overlay() functio to add both rasters together. 
leastcost_corridor <- overlay(Bottling.cost,Etherow.cost, fun=function(x, y){return(x + y)})

#plot
plot(leastcost_corridor, legend=F, axes=F)
#Bottling
plot(Bottling, add=T)
#Etherow
plot(Etherow,add=T)

#get lower quantile (I have chosen 5 here which seemed to work well but this number isn't special!)
quantile5 <- quantile(leastcost_corridor, probs=0.05, na.rm=TRUE)

Now we have our values to use for subsetting we use the ifelse() function to do the subsetting. The first argument is a condition that is test, for us this testing whether values are greater than the 5th percentile of cost values, and the following two arguments decide what happens if the statememt is true. We will set values to NA if this statement is true (because the cost is higher than our threshold) and to 1 is the statement is false.

#make new truncated layer
leastcost_corridorVals=ifelse(values(leastcost_corridor>quantile5),NA,1)

leastcost_corridor[]=leastcost_corridorVals

#plot
plot(GM)
plot(leastcost_corridor,add=T,col="blue")

Now we have our least cost corridor. Note that this corridor is still based on a single interpretation of landscape cost and assumes complete homomgeneity within the corridor. In later weeks we will look at ways to better capture environmental heterogeneity in the way we delineate habitat and related metrics such as movement cost and edge effects. For now, you have reached the end of this week’s practical. Well done.