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 at the centre 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).

This practical will explore how we can take measurements related to some of these concepts using basic commands in R and with the help of R packages (raster and SDMTools). As there was no lecture this week I strongly recommend you read chapter 3 of the core text for this module (Fetcher and Fortin, 2018) which you can access here: https://link.springer.com/chapter/10.1007/978-3-030-01989-1_3 (you should be able to download the content via Shibboleth using your university credentials).

Another useful source for reference is the Fragstats documentation on landscape metrics. Fragstats is a stand alone landscape ecology application developed at the University of Massachussetts and most subsequent packages tend to adopt the same definitions of individual metrics. Therefore the original Fragstats documentation continues to be a handy reference guide (and can be found here: http://www.umass.edu/landeco/research/fragstats/documents/Metrics/Metrics%20TOC.htm). A much more detailed document can be found here (http://www.umass.edu/landeco/research/fragstats/documents/fragstats.help.4.2.pdf). The latter is very long but useful as a reference for looking up individual metrics.

First let’s install the packages for doing the work:

setwd("P:/GEOG70922/Week7")


#install necessary packages

install.packages(c("raster","SDMTools","rasterVis","landscapemetrics"))

Now load all the libraries needed for the first part of the practical.

library(SDMTools)
library(rasterVis)
## Warning: package 'rasterVis' was built under R version 3.6.3
## Loading required package: raster
## Warning: package 'raster' was built under R version 3.6.2
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.6.2
## 
## Attaching package: 'raster'
## The following object is masked from 'package:SDMTools':
## 
##     distance
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer

Loading the data

The code used in this practical is largely borrowed from Chapter 3 of Fletcher and Fortin but adapted to work with some more familiar data. Copy the GM.tif from the Week7 folder in the GEOG70922 area of the T Drive to a folder on your P drive for this work and then load directly into R:

library(raster)# raster library for reading and working with raster objects.


#load some landcover for our analysis

GM<- raster("GM.tif")


plot(GM)# have a look at the data

In order to visualize our categorical raster better (beyond labels of 1-5) we’ll turn our raster into a factor (remember a factor in R is a way to describe categories). In the case of our GM raster we will convert the values to a factor, the original raster values will then each be “levels” in our categorical raster. By giving names to each level, we can visualize them when plotting. Along the way we’ll save ourselves some processing time for the practical work by aggregating our raster to a lower resolution (larger cell size).

GM<-aggregate(GM, fact=3, fun=modal) #increase cell size from 10 to 30 metres


GM<-as.factor(GM) # convert to categorical raster

levels(GM)#inspect levels
## [[1]]
##   ID
## 1  1
## 2  2
## 3  3
## 4  4
## 5  5
res(GM)#check resolution
## [1] 30 30
#create a holder for our raster leves (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","field layer","canopy")

#tell R that these are what we want to name each level

levels(GM) <- land_cover



#create a list of colours for plotting

land_col <- c("white","blue","green","yellow","dark green")
plot(GM, legend = T, col = land_col)

#we can plot the map also with the rasterVis package to visualize our labels

levelplot(GM, col.regions=land_col, xlab="", ylab="")

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

GM.cat <- unique(GM)# this holds a list of the unique values in the GM raster (i.e. 1,2,3,4,5).

gmF <- c(0,0,0,0,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.

reclass.mat <- cbind(GM.cat,gmF)# create a matrix for reclassification 
reclass.mat#inspect
##      GM.cat gmF
## [1,]      1   0
## [2,]      2   0
## [3,]      3   0
## [4,]      4   0
## [5,]      5   1
#forest binary layer from reclassification matrix
gmForest <- reclassify(GM,reclass.mat) #use the recassify function passing the raster as the first argument and the reclassification matrix as the second.  
plot(gmForest)#plot it!

More trees than you thought in Greater Manchester?

Now, with this new canopy layer we can carry out some basic analysis using some common landscape metrics using the SDMTools 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 is typically seen as one pixel width) and shape (the complexity of patch shape can influence ecological processes such as the severity of edge effects). See the fragstat documentation for more details on how these metrics area calculated.

Currently our raster can be considered an arrary of cells of value 0 and 1. So an initial task is to convert these individual pixels into groups of pixels that are connected (i.e. that form patches). This is not a trivial matter (try creating such a network in ArcGIS or python). Thankfully the SDMTools package in R comes with the clump() function which carries this out for us in a neat way. For this function we pass the raster as the first argument, followed by the number of surrounding cells we wish to consider when establishing connectedness. In the below code we state ‘8’ to use the queen neighbourhood around each cell.

############################################
#patch-level analysis
############################################

#8-neighbour rule
gmF8 <- clump(gmForest, directions=8)
## Loading required namespace: igraph
plot(gmF8)#plot the result

cellStats(gmF8, max)#number of patches identified
## [1] 28257
#Now, calculate patch-level metrics:

?PatchStat #call the function info to the 'help' panel
## starting httpd help server ...
##  done
# now use the PatchStat() function to compute a range of metrics. Note with large data sets the function works best when the raster is firt converted to a matrix of values.

gmFstat <- PatchStat(as.matrix(gmF8), cellsize=30)

#list the names of the metrics we have just computed.
names(gmFstat)
##  [1] "patchID"           "n.cell"            "n.core.cell"      
##  [4] "n.edges.perimeter" "n.edges.internal"  "area"             
##  [7] "core.area"         "perimeter"         "perim.area.ratio" 
## [10] "shape.index"       "frac.dim.index"    "core.area.index"

The number of patches in our data are numerous and individual values are not particularly meaningful for us so we will calculate the mean of each metric as a summary:

#get the mean of each column (excluding column 1 - ID)

gmFstat.mean <- colMeans(gmFstat[,2:ncol(gmFstat)])
gmFstat.mean
##            n.cell       n.core.cell n.edges.perimeter  n.edges.internal 
##      9.500301e+00      1.143150e+00      1.644470e+01      2.155650e+01 
##              area         core.area         perimeter  perim.area.ratio 
##      8.550271e+03      1.028835e+03      4.933411e+02      1.140507e-01 
##       shape.index    frac.dim.index   core.area.index 
##      1.246140e+00      1.042122e+00      2.807206e-03
#plot area frequencies on a log scale - this exaggerates changes at the lower end of area values which is useful as we are generally more interested in changes from say, 100 to 200 than from 1000 to 1100 m2. 
hist(log(gmFstat$area))

#correlation matrix
round(cor(gmFstat[,6:12]),2)
##                   area core.area perimeter perim.area.ratio shape.index
## area              1.00      0.91      0.98            -0.18        0.74
## core.area         0.91      1.00      0.82            -0.13        0.57
## perimeter         0.98      0.82      1.00            -0.19        0.78
## perim.area.ratio -0.18     -0.13     -0.19             1.00       -0.46
## shape.index       0.74      0.57      0.78            -0.46        1.00
## frac.dim.index    0.26      0.18      0.29            -0.64        0.75
## core.area.index   0.39      0.44      0.35            -0.35        0.44
##                  frac.dim.index core.area.index
## area                       0.26            0.39
## core.area                  0.18            0.44
## perimeter                  0.29            0.35
## perim.area.ratio          -0.64           -0.35
## shape.index                0.75            0.44
## frac.dim.index             1.00            0.28
## core.area.index            0.28            1.00
#plot correlation across for key metrics

pairs(gmFstat[,c("area","core.area", "perim.area.ratio", "shape.index")])

As we can see, most of our patch-level metrics are concerned with the area and perimeter of individual patches. These factors are important from an ecological and biogeographical perspective as they can sharply infuence species richness within habitats as a function of disturbance and (meta-)population dynamics.

We can also see that perimeter/area ratio and shape index area weekly correlated and that the latter correlates well with patch area, unlike the perimeter/area ratio metric. This is due to the fact the shape index accounts for area in the way it is calculated. In essence Shape Index seeks to measure the “compactness” of any given patch (where the smallest patch, here a square comprising a single cell, is the most compact shape possible).

Next we’ll go “up” one level to consider class-level characteristics.

############################################
#3.3.3.2 Class-level quantification
############################################

# check the documentation
?ClassStat

#calculation based on forest layer
GMFcstat <- ClassStat(as.matrix(gmForest), cellsize=res(gmForest)[[1]])

#let's just return a selection of the results 
GMFcstat[,c(2:4,19,37,38)]
##   n.patches total.area prop.landscape mean.shape.index effective.mesh.size
## 1      1876 1037401200      0.8110994         1.214419         827833210.0
## 2     28257  241605000      0.1889006         1.246140            255162.2
##   patch.cohesion.index
## 1             9.982225
## 2             9.693244
#calculation based on GM layer (all land-cover types)
GM.cstat <- ClassStat(as.matrix(GM), cellsize=res(GM)[[1]])

head(GM.cstat)
##   class n.patches total.area prop.landscape patch.density total.edge
## 1     1      7555  423204300      0.3308853  5.906930e-06   12562680
## 2     2      3192   21621600      0.0169050  2.495688e-06    1146720
## 3     3     23057  173291400      0.1354891  1.802728e-05    8573040
## 4     4     28725  419283900      0.3278201  2.245884e-05   14454180
## 5     5     28257  241605000      0.1889006  2.209293e-05   13940340
##   edge.density landscape.shape.index largest.patch.index mean.patch.area
## 1 0.0098222198             152.60787        0.2370549103       56016.453
## 2 0.0008965711              61.65161        0.0006987456        6773.684
## 3 0.0067028917             162.73804        0.0019097640        7515.783
## 4 0.0113011024             176.35652        0.0953622430       14596.480
## 5 0.0108993530             224.04918        0.0061008305        8550.271
##   sd.patch.area min.patch.area max.patch.area perimeter.area.frac.dim
## 1    3496644.49            900      303194700              0.05936911
## 2      34858.19            900         893700              0.10606340
## 3      37421.87            900        2442600              0.09893643
## 4     743940.79            900      121968900              0.06894420
## 5     107129.80            900        7803000              0.11539305
##   mean.perim.area.ratio sd.perim.area.ratio min.perim.area.ratio
## 1             0.1054213          0.03075563          0.012009238
## 2             0.1109722          0.02820588          0.007519302
## 3             0.1090785          0.03184678          0.012030075
## 4             0.1131639          0.02730220          0.009919028
## 5             0.1140507          0.02518630          0.015942029
##   max.perim.area.ratio mean.shape.index sd.shape.index min.shape.index
## 1            0.1333333         1.245472      1.4754591               1
## 2            0.1333333         1.215524      0.4790930               1
## 3            0.1333333         1.165902      0.4036288               1
## 4            0.1333333         1.203364      0.6076525               1
## 5            0.1333333         1.246140      0.7416017               1
##   max.shape.index mean.frac.dim.index sd.frac.dim.index min.frac.dim.index
## 1      112.610680            1.040071        0.05616369                  1
## 2        9.692308            1.043862        0.05913331                  1
## 3       10.085714            1.032838        0.04937338                  1
## 4       30.010855            1.036967        0.05647489                  1
## 5       27.417112            1.042122        0.06246819                  1
##   max.frac.dim.index total.core.area prop.landscape.core mean.patch.core.area
## 1           1.483780       131274900         0.102638205            17375.897
## 2           1.358529         5346900         0.004180511             1675.094
## 3           1.316182        32471100         0.025387758             1408.297
## 4           1.373185       162022500         0.126678432             5640.470
## 5           1.429040        29071800         0.022729991             1028.835
##   sd.patch.core.area min.patch.core.area max.patch.core.area
## 1         1193310.35                   0           103401900
## 2           20529.73                   0              687600
## 3           14615.94                   0              913500
## 4          511540.55                   0            84982500
## 5           24377.77                   0             2009700
##   prop.like.adjacencies aggregation.index lanscape.division.index
## 1             0.6358111          77.85007               0.9435263
## 2             0.4308517          60.61419               0.9999975
## 3             0.4587481          63.03987               0.9999795
## 4             0.5891288          74.25374               0.9902783
## 5             0.3959251          56.83561               0.9998005
##   splitting.index effective.mesh.size patch.cohesion.index
## 1        17.70736        72230188.607             9.969156
## 2    406543.16903            3146.053             9.300349
## 3     48700.77820           26262.541             9.238968
## 4       102.86236        12434152.146             9.931698
## 5      5012.52175          255162.225             9.693244
#check against PatchStat calculations:

#mean patch size for canopy layer ('5' in original raster)
GM.cstat[GM.cstat$class==5, "mean.patch.area"]
## [1] 8550.271
gmFstat.mean["area"]#mean patch size
##     area 
## 8550.271
#correlation matrix
cor(GM.cstat[,c(2:4,19,37,38)])#subset of metrics
##                       n.patches total.area prop.landscape mean.shape.index
## n.patches             1.0000000  0.3556487      0.3556487       -0.2378791
## total.area            0.3556487  1.0000000      1.0000000        0.2936711
## prop.landscape        0.3556487  1.0000000      1.0000000        0.2936711
## mean.shape.index     -0.2378791  0.2936711      0.2936711        1.0000000
## effective.mesh.size  -0.4206692  0.6596270      0.6596270        0.4880187
## patch.cohesion.index  0.1974428  0.9154469      0.9154469        0.6077973
##                      effective.mesh.size patch.cohesion.index
## n.patches                     -0.4206692            0.1974428
## total.area                     0.6596270            0.9154469
## prop.landscape                 0.6596270            0.9154469
## mean.shape.index               0.4880187            0.6077973
## effective.mesh.size            1.0000000            0.6629632
## patch.cohesion.index           0.6629632            1.0000000
#plot subset of metrics
pairs(GM.cstat[,c(2:4,19,37,38)])

As you can see, many of the metrics produced here are aggregations of those computed for the patch-level but summarized for each category in our raster object. Others relate to concepts such as connectivity (e.g. effective mesh size and patch cohesion) and density (patch density) which take into account the spatial distribution of patches rather than simple measures of size and shape. Connectivity, for example, is not meaningful when applied at the patch-level and so is considered a class or landscape-level metric. Look at the suggested reading for this practical for more detail on these metrics. Briefly, effective mesh size is the probability that two randomly chosen points in a landscape will occur in the same patch multiplied by the total area of the study region. A nice description by the original author can be found here: https://tinyurl.com/sguyvvl and the original paper (for those interested) here: https://link.springer.com/article/10.1023/A:1008129329289

Patch cohesion is another measure of connectivity that reflects the “clumpiness” of habitat landcover patches and is influenced more by physical connectedness and the size and shape of individual patches (more here: https://esajournals.onlinelibrary.wiley.com/doi/abs/10.2307/2265590).

So far the metrics that we have calculated are informative but are essentially singular numerical summaries of processes occuring across space. We do have a few options available however for projecting our results on a map by using some basic commands in R. Here we will consider core area defined as the patch minus the edge. We can simulate this effect by generating internal buffers for each patch. In other words, we can convert our patches to polygons and buffer inside their extent to estimate core area.

#-----------------------------------#
#core area and edge distances
#-----------------------------------#



forest.poly <- rasterToPolygons(gmF8, dissolve=T) #convert to vector format and dissole patches into single polygons
## Loading required namespace: rgeos
#plot
plot(forest.poly)

library(rgeos)#for buffer analysis
## Warning: package 'rgeos' was built under R version 3.6.2
## rgeos version: 0.5-2, (SVN revision 621)
##  GEOS runtime version: 3.6.1-CAPI-1.10.1 
##  Linking to sp version: 1.3-2 
##  Polygon checking: TRUE
#manual calculation of core area
core.poly <- gBuffer(forest.poly, width = -10, byid = T) #buffer the inside each polygon
core.area <- gArea(core.poly, byid=T)#get the core area

hist(log(core.area)) #plot frequencies of core areas

#plot
plot(forest.poly, col="green",border=NA)
plot(core.poly, col="blue",border=NA,add=TRUE)

Another category of metric that we are able to compute using the raster package relates to measures of isolation (an important measure from a landscape connectivity perspective).

The distance() function can be used to caculate distance (i.e. isolation) of each canopy cell to NA cell in our raster layer. We therefore need to set all non-canopy cells to NA before using this function.

The distance() function is computationally quite intensive to use so we will crop our study area to something more manageable:

extent.new<-extent(375000,380000,405000,410000)


treeCrop<-crop(gmForest,extent.new)

plot(treeCrop)

#calculate distance to edge with raster package
forestNA <- treeCrop
forestNA[forestNA == 0] <- NA
forest.dist <- raster::distance(forestNA)

#plot distance to forest edge
plot(forest.dist)

There is one ‘level’ left at which we can characterise our landcover data - the landscape level. Landscape level metrics are based on concepts that involve relationships between patches and/or classes in the landscape. In some cases metrics can be seen as both class and landcape-level appropriate. For example, measures of connectivity such as effective mesh size and patch cohesion (see below) can be calculated at the class-level (e.g. for each land-cover type in our example) at the landscape (study area) scale and so can be considered class-landcape metrics. Others, such as richness (or diversity more generally), can only be defined as landscape metrics as diversity metrics consider multiple classes and are calculated at the landscape-scale (however we choose to define the landscape). Below we will consider a landscape scale analysis in two ways. Firstly we will calculate the effective mesh size of all green-blue land-cover types in the landscape by reclassifying our raster and using the ClassStat() function on this reclassified layer. Essentially we are running a class-level analysis but re-classifying our data beforehand to consider all green-blue landcover as a single feature in the landscape.

Secondly we calculate measures of diversity by computing the richness (number of landcovers) in the landscape and then Shannon’s Diversity index based on the relative distribution of each of these landcover types.

############################################
#3.3.3.3 landscape-level quantification
############################################

#some summary metrics derived from class-level metrics



#create a reclassification matrix
gmG <- c(0,1,1,1,1)

reclass.mat <- cbind(GM.cat,gmG)
reclass.mat
##      GM.cat gmG
## [1,]      1   0
## [2,]      2   1
## [3,]      3   1
## [4,]      4   1
## [5,]      5   1
#forest binary layer from reclassification matrix
gmGreen <- reclassify(GM,reclass.mat)
plot(gmGreen)

#perform class stats on that aggregated green-blue landcover layer
GMGcstat <- ClassStat(as.matrix(gmGreen), cellsize=res(gmGreen)[[1]])




GMGcstat$n.patches
## [1]  7555 12142
GMGcstat$effective.mesh.size
## [1]  72230189 474537082
GMGcstat$patch.cohesion.index
## [1] 9.969156 9.980230
#----------------------------------#
#diversity-related metrics
#----------------------------------#






#richness
richness <- length(unique(values(GM)))
richness
## [1] 6
#ah, the NAs have been included but we don't want those!

#richness function to calulate richness that removes NAs
richness <- function(x) length(unique(na.omit(x)))

richness(GM)
## [1] 5

We now have a simple measure of richness oflandcover type but this is really only repeating information that we already knew (that there are 5 landcovers in our dataset). 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, and performs well in this regard, is the Shannon’s diversity index. We will look at how this is calculated below (but do see the suggested reading for other indices worth exploring).

The basic approach used by the Shannon index (SHDI) is to take each class and calculate the proportion of the landscape attributed to it, multiply this number by its logarithm then sum the result for all values. 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). Take a look at the entry for SHDI in the Fragstats documentation in the link provided at the start of this practical. The mathematical notation for the index is:

H’ = - ∑pi ln(pi)

Where:

pi is the proportion of each class and ln(pi) is the natural logarithm of the proportion of each class.

# Shannon's diversity index ,D

table(values(GM))
## 
##      1      2      3      4      5 
## 470227  24024 192546 465871 268450
C <- table(values(GM))


P <- C / sum(C)

D <- -sum(P * log(P))

D
## [1] 1.386177

Again, using the SDMTools package we have simply calculated a single metric that attempts to reflect what is happening at the landscape scale. It would be nice if we could acheive a more localized version of events that tells us, for a given point in space, the diversity of land-cover in its neighbourhood. If you consider the goal of most of our analysis on this course up to here, you can surely appreciate why this would be useful.

Unfortunately we cannot acheive this with a ready-made function from the SDMTools package but we can use the focal() function in the raster package to implement our own. We will use two main functions to do this. The first is the focalWeight() function and the second is the focal() function itself. focalWeight() allows us to describe an area around each cell in the our raster layer - a neighbourhood of n x n cells (if a rectangular shape is chosen) or a buffer of radius n (if a circular shape i chosen). The focalWeight() function then assigns values to all cells in this neighbourhood such that their sum equals 1. Type ?focalWeight into the console and press return to get more help on this function.

To get used to working with the focal() function we will begin by caculating the proportion of forest cover within a 100m buffer around each cell (before moving on to the more complex function for calculating SHDI).

Note that R reads the buffer as a matrix that below we re-scale to a binary 1/0 format. The rationale for doing this is that, when we apply our function, R will treat everything outside of our buffer as zero and multiply everything inside the buffer by 1. The default setting for the focal() function is to sum all the values and assign this to the focal cell. In this way the code written below will perform a calculation that returns values for each cell that reflect the total number of cells assigned to canopy cover in a 100 m radius.

#focal buffer matrix for moving windows
buffer.radius <- 100
fw.100m <- focalWeight(gmForest, buffer.radius, 'circle')#buffer in CRS units

#inspect the focalWeight object
fw.100m
##            [,1]       [,2]       [,3]       [,4]       [,5]       [,6]
## [1,] 0.00000000 0.00000000 0.02702703 0.02702703 0.02702703 0.00000000
## [2,] 0.00000000 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703
## [3,] 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703
## [4,] 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703
## [5,] 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703
## [6,] 0.00000000 0.02702703 0.02702703 0.02702703 0.02702703 0.02702703
## [7,] 0.00000000 0.00000000 0.02702703 0.02702703 0.02702703 0.00000000
##            [,7]
## [1,] 0.00000000
## [2,] 0.00000000
## [3,] 0.02702703
## [4,] 0.02702703
## [5,] 0.02702703
## [6,] 0.00000000
## [7,] 0.00000000
#re-scale weight matrix to 0,1
fw.100m <- ifelse(fw.100m > 0, 1, 0)

#inspect the re-scaled object
fw.100m
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]    0    0    1    1    1    0    0
## [2,]    0    1    1    1    1    1    0
## [3,]    1    1    1    1    1    1    1
## [4,]    1    1    1    1    1    1    1
## [5,]    1    1    1    1    1    1    1
## [6,]    0    1    1    1    1    1    0
## [7,]    0    0    1    1    1    0    0
#forest cover moving window; number of cells
forest.100m <- focal(gmForest, w=fw.100m)



plot(forest.100m)

Now, to transform this raster to one that gives us a nice standardized score (let’s take proportion here) we simply divide the forest.100m object by the total number of cells in the buffer:

#proportion of forest
forest.prop.100m <-forest.100m/sum(fw.100m) #proportion: divide by buffer size
plot(forest.prop.100m)

We can also compute a buffer based on Gaussian kernel (to assign more weight to nearby cells in our subsequent calculations). The first number inside the function is the sigma (denoting the kernel width - in the case of a Gaussian kernel this is actually the standard deviation) and the second is the radius of the buffer:

#weight matrix for a Gaussian kernel
focalGauss<-focalWeight(gmForest, c(3,100), type = "Gaus")

focalGauss

We now have a sense of canopy density at a more localised scale. Think how this might be incorporated into a species-based analysis such as the one we carried out on scale responses in Week 2.

Now let’s look at calculating SHDI 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 or temperature in urban areas).

Remember that the basic premise is to start from a position of considering the proportion of each landcover in the study area. Conceptually then, what we need to do is similar to the above calculation for canopy cover. We just need to calculate this proportion as above but for each landcover class. We then create a function that returns the logarithm of each cell for a given landcover class 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 calulation (by adding each layer to this blank raster in turn).

#diversity moving window function
diversity <- function(landcover, radius) {    # this function will take two arguments in order to run: a landcover raster and a radius (a number)
  
  n=unique(landcover) # create an object that signifies the number of unique layers in the landcover raster (i.e. the number 5 in our example)
  
  #Create focal weights matrix
  fw.i <- focalWeight(landcover, radius, type="circle")
  fw.i<- ifelse(fw.i > 0, 1, 0) # rescale as before
  
  #Compute focal means of indicators
  D <- landcover # blank raster to hold all (sum) values
  values(D) <- 0
  log.i <- function(x) ifelse(x==0, 0, x*log(x))# this is the function that will carry out log(p)*p
  
  #for each landcover, create a moving window map and make calculations
  for (i in 1:length(n)) {
    focal.i <- focal(landcover == i, fw.i,na.rm=FALSE) # this line sums cells of each landcover class
    focalProp.i<-  focal.i/sum(fw.i,na.rm=FALSE) #this line calculate proportion as in the example on canopy above
    D <- D + calc(focalProp.i, log.i)#take log(p)*p of each proportion value and add to the blank raster
    #plot(D)
  }
  
  D <- calc(D, fun=function(x){x * -1})#take negative
  return(D)
}

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

You’ve reached the end of this practical session. If you want to practice further you could revisit some of the landscape datasets used in previous weeks and apply the functions covered above. You could also revisit the analysis from Week 2 based on your new knoweldge on creating weighted kernel-based neighbourhoods to see how using a kernel window stacks up against the distance-based buffer.