Practical 1: Landscape Metrics

library(landscapemetrics)
library(landscapetools)
library(raster)
library(tidyverse)
library(janitor)

## Import a raster and plot it. 
esett <- raster("data/esett.asc")
psett <- raster("data/psett.asc")
show_landscape(esett, discrete = TRUE)

Task 1: Landscape composition

Task 1.1 Calculate the total area of each class and use this to calculate the proportion of each class type in each landscape. Summarise in Table 1.1a and 1.1b (2 pts). You may want to switch to the Visual Editor to fill in these tables.

esett_metrics <-calculate_lsm(esett, #This is our input raster
                              what = c("lsm_c_ca")) #class area for each cover type

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

esett_metrics #This will print out the table

## # A tibble: 3 × 6
##   layer level class    id metric value
##   <int> <chr> <int> <int> <chr>  <dbl>
## 1     1 class     1    NA ca        72
## 2     1 class     2    NA ca        23
## 3     1 class     3    NA ca         5

psett_metrics <-calculate_lsm(psett, #This is our input raster
                              what = c("lsm_c_ca")) #class area for each cover type

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

psett_metrics #This will print out the table

## # A tibble: 3 × 6
##   layer level class    id metric value
##   <int> <chr> <int> <int> <chr>  <dbl>
## 1     1 class     1    NA ca        33
## 2     1 class     2    NA ca        47
## 3     1 class     3    NA ca        20

Table 1.1a – Early Settlement

Cover type	Cells	Pi
Forested	72	.72
Agricultural	23	.23
Urban	5	.5
—————	———	——-
Total	100

Table 1.1b – Post Settlement

Cover type	Cells	Pi
Forested	33	.33
Agricultural	47	.47
Urban	20	.2
—————	———	——-
Total	100

Task 1.2 Calculate Dominance: S is the number of land cover types Pi is the proportion of the ith land cover type. Values near 1 indicate the landscape is dominated by one land cover type; a value of zero means all cover types are equally abundant. Do this using the equation below. Summarise in Table 2 (2 pts).

Hint: use the log() and sum() functions to recreate this function in R.

(log(3)+sum((0.72*log(0.72))+(0.23*log(0.23))+(0.05*log(0.05))))/log(3)

## [1] 0.3406819

(log(3)+sum((0.33*log(0.33))+(0.47*log(0.47))+(0.20*log(0.20))))/log(3)

## [1] 0.05097834

Task 1.3 Using the landscapemetrics package in R calculate the Shannon’s Evenness Index: Values near 0 indicate the landscape is dominated by one land cover type; a value of 1 means all cover types are equally abundant. It is a measure of dominance. Summarise in Table 2 (2 pts).

#Task 1.3 code here
esett_metrics <-calculate_lsm(esett, #This is our input raster
                              what = c("lsm_l_shei"))#Shannon's evenness index

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

esett_metrics #This will print out the table

## # A tibble: 1 × 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA shei   0.659

psett_metrics <-calculate_lsm(psett, #This is our input raster
                              what = c("lsm_l_shei"))#Shannon's evenness index

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

psett_metrics #This will print out the table

## # A tibble: 1 × 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA shei   0.949

Table 1.3

Metric	Pre-Settlement	Post-settlement
Dominance	0.3406819	0.05097834
Shannon’s Evenness Index	0.6593181	0.9490217

Question 1.4 What is the relationship between Dominance and Shannon’s Evenness Index? Do you need to analyse both (1 pts)?

Shannon’s Evenness and Dominance are inversely related: as dominance increases, evenness decreases

Question 1.5 How would you interpret/describe the changes in this landscape between the two time periods (2 pts)?

After settlement, the landscape is more evenly distributed and less dominated by one class.

Task 1.6 Using Table 1.6 below calculate Shannon’s Evenness Index for each landscape. You will use these to answer the next questions in your notebook (2 pts).

#Task 1.6 code here
1-(log(4)+sum((0.1*log(0.1))+(0.8*log(0.8))+(0.1*log(0.1))))/log(4)

## [1] 0.460964

1-(log(4)+sum((0.8*log(0.8))+(0.1*log(0.1))+(0.1*log(0.1))))/log(4)

## [1] 0.460964

1-(log(4)+sum((0.65*log(0.65))+(0.2*log(0.2))+(0.15*log(0.15))))/log(4)

## [1] 0.639449

1-(log(4)+sum((0.65*log(0.65))+(0.2*log(0.2))+(0.15*log(0.15))))/log(4)

## [1] 0.639449

Table 1.6

Landscape	Pforest	Pagricultural	Purban	Shannon’s Evenness Index
W	0.1	0.8	0.1	0.460964
X	0.8	0.1	0.1	0.460964
Y	0.65	0.2	0.15	0.639449
Z	0.15	0.2	0.65	0.639449

Question 1.7 Which of these hypothetical landscapes might be considered “similar” when only comparing Shannon Evenness Index (SHEI)? Considering your interpretation of the data in Table 3, what other types of information and/or metrics would be necessary to distinguish these landscapes? (4 pts)

W and X are similar, and Y and Z are similar, because they have the same SHEI values.

To tell them apart, you need information on which class is dominant and other metrics like patch density or edge density

Task 2: Landscape structure

Task 2.1 Using the landscapemetrics package in R calculate the Number of Patches and Mean Patch Size for each landscape and summarise in Table 2.1 (4 pts).

#Task 2.1 code here
esett_metrics <-calculate_lsm(esett, #This is our input raster
                              what = c("lsm_c_np", #number of patches
                                        "lsm_c_area_mn"), #Mean patch area
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

esett_metrics

## # A tibble: 6 × 6
##   layer level class    id metric  value
##   <int> <chr> <int> <int> <chr>   <dbl>
## 1     1 class     1    NA area_mn 72   
## 2     1 class     2    NA area_mn  4.6 
## 3     1 class     3    NA area_mn  1.67
## 4     1 class     1    NA np       1   
## 5     1 class     2    NA np       5   
## 6     1 class     3    NA np       3

psett_metrics <-calculate_lsm(psett, #This is our input raster
                              what = c("lsm_c_np", #number of patches
                                        "lsm_c_area_mn"), #Mean patch area
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

psett_metrics

## # A tibble: 6 × 6
##   layer level class    id metric  value
##   <int> <chr> <int> <int> <chr>   <dbl>
## 1     1 class     1    NA area_mn  3   
## 2     1 class     2    NA area_mn  5.22
## 3     1 class     3    NA area_mn  4   
## 4     1 class     1    NA np      11   
## 5     1 class     2    NA np       9   
## 6     1 class     3    NA np       5

Table 2.1

Class type	# patches - Early	# patches - Post	Mean patch size - Early	Mean patch size - Post
Forested	72	3	1	11
Agricultural	4.6	5.2222	5	9
Urban	1.667	4	3	5

Question 2.2 What characteristics of a landscape do the number of patches and the average patch size represent (2 pts)?

The number of patches shows how fragmented the landscape is, while mean patch size shows how large and continuous the patches are.

Task 2.3 You will also need to calculate the edge/area ratio (i.e. Number of Edges of class I /Total Area of class I) for each class for both the early and post settlement time periods. You can do this manually based on the outputs from the code above. Summarise these in Table 2.3 (4 pts).

#Task 2.3 code here
esett_metrics <-calculate_lsm(esett, #This is our input raster
                              what = c("lsm_c_ed",# Edge density
                                       "lsm_c_ca"), #Class area 
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

esett_metrics

## # A tibble: 6 × 6
##   layer level class    id metric value
##   <int> <chr> <int> <int> <chr>  <dbl>
## 1     1 class     1    NA ca        72
## 2     1 class     2    NA ca        23
## 3     1 class     3    NA ca         5
## 4     1 class     1    NA ed        29
## 5     1 class     2    NA ed        43
## 6     1 class     3    NA ed        16

psett_metrics <-calculate_lsm(psett, #This is our input raster
                              what = c("lsm_c_ed",# Edge density
                                       "lsm_c_ca"), #Class area 
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

psett_metrics

## # A tibble: 6 × 6
##   layer level class    id metric value
##   <int> <chr> <int> <int> <chr>  <dbl>
## 1     1 class     1    NA ca        33
## 2     1 class     2    NA ca        47
## 3     1 class     3    NA ca        20
## 4     1 class     1    NA ed        58
## 5     1 class     2    NA ed        86
## 6     1 class     3    NA ed        40

Table 2.3

Class type	Number of Edges- Early	Number of Edges - Post	Edge:Area Ratio- Early	Edge:Area Ratio - Post
Forested	29	58	0.40	1.76
Agricultural	43	86	1.87	1.83
Urban	16	40	3.20	2.00

Question 2.4 Two landscapes are the same size and both contain the same amount of a given cover type. Landscape A has four patches of that cover type, and Landscape B has 17 patches of the same cover type. Which of the landscapes will have the higher edge to area ratio of that cover type (1 pt)?

Landscape B will have the higher edge-to-area ratio because more patches mean more edges for the same total area.

Task 2.5 Use landscapemetrics package in R to calculate contagion and summarise in Table 2.5 (2 pts).

#Task 2.5 code here
esett_metrics <-calculate_lsm(esett, #This is our input raster
                              what = c("lsm_l_contag", "lsm_c_np"), #contag.with NP, can't run a single metric
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

esett_metrics

## # A tibble: 4 × 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 class         1    NA np       1  
## 2     1 class         2    NA np       5  
## 3     1 class         3    NA np       3  
## 4     1 landscape    NA    NA contag  41.8

#alternative code to run one metric only:

esett_metrics<-lsm_l_contag(esett)
esett_metrics

## # A tibble: 1 × 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA contag  41.8

psett_metrics <-calculate_lsm(psett, #This is our input raster
                              what = c("lsm_l_contag", "lsm_c_np"), #contag.with NP, can't run a single metric
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

psett_metrics

## # A tibble: 4 × 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 class         1    NA np     11   
## 2     1 class         2    NA np      9   
## 3     1 class         3    NA np      5   
## 4     1 landscape    NA    NA contag  8.78

#alternative code to run one metric only:

psett_metrics<-lsm_l_contag(psett)
psett_metrics

## # A tibble: 1 × 6
##   layer level     class    id metric value
##   <int> <chr>     <int> <int> <chr>  <dbl>
## 1     1 landscape    NA    NA contag  8.78

Table 2.5

Cover Type	Pre-Settlement	Post-Settlement
Contagion	41.82543	8.780757

Question 2.6 What does the contagion value tell you about changes in landscape structure over time (2 pts)?

Contagion decreased from pre to post, showing the landscape became more fragmented and less connected over time, with patches more dispersed and mixed.

Task 3: Neighbourhood rules

Task 3.1 Adapt the provided code to calculate the suite of landscape metrics for both neighbourhood rules. Once you have created four objects with each set of metrics, you can use the provided code to bind them into a table. Put all of your code here: (4 pts):

#Paste task 3.1 code here
esett_4 <-calculate_lsm(esett, #This is our input raster
                              what = c("lsm_c_ca", #Class area
                                       "lsm_c_np",#number of patches
                                       "lsm_c_area_mn", #Mean patch area
                                       "lsm_c_ed",# Edge density
                                       "lsm_l_shdi",#Shannon's diversity index
                                       "lsm_l_contag"), #contagion
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

esett_4$layer <- "esett4" #This gives the layer a name associated with the map and neighborhood rule so we can bind them!
psett_4 <-calculate_lsm(psett, #This is our input raster
                              what = c("lsm_c_ca", #Class area
                                       "lsm_c_np",#number of patches
                                       "lsm_c_area_mn", #Mean patch area
                                       "lsm_c_ed",# Edge density
                                       "lsm_l_shdi",#Shannon's diversity index
                                       "lsm_l_contag"), #contagion
                              directions = 4) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

psett_4$layer <- "psett4" #This gives the layer a name associated with the map and neighborhood rule so we can bind them!

esett_8 <-calculate_lsm(esett, #This is our input raster
                              what = c("lsm_c_ca", #Class area
                                       "lsm_c_np",#number of patches
                                       "lsm_c_area_mn", #Mean patch area
                                       "lsm_c_ed",# Edge density
                                       "lsm_l_shdi",#Shannon's diversity index
                                       "lsm_l_contag"), #contagion
                              directions = 8) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

esett_8$layer <- "esett8" #This gives the layer a name associated with the map and neighborhood rule so we can bind them!

psett_8 <-calculate_lsm(psett, #This is our input raster
                              what = c("lsm_c_ca", #Class area
                                       "lsm_c_np",#number of patches
                                       "lsm_c_area_mn", #Mean patch area
                                       "lsm_c_ed",# Edge density
                                       "lsm_l_shdi",#Shannon's diversity index
                                       "lsm_l_contag"), #contagion
                              directions = 8) #This is the neighborhood rule we are using

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

psett_8$layer <- "psett8" #This gives the layer a name associated with the map and neighborhood rule so we can bind them!

allmetrics <- list(esett_4, esett_8, psett_4, psett_8) #This creates a list of the data tables
metric_table <- do.call(rbind, allmetrics) #This binds the tables together
metric_table <- tidyr::spread(metric_table, layer, value) #This puts the table into wide format

metric_table[order(metric_table$level, metric_table$metric),] #This prints the table and orders it

## # A tibble: 14 × 8
##    level     class    id metric  esett4 esett8 psett4 psett8
##    <chr>     <int> <int> <chr>    <dbl>  <dbl>  <dbl>  <dbl>
##  1 class         1    NA area_mn 72     72       3      4.12
##  2 class         2    NA area_mn  4.6    7.67    5.22  15.7 
##  3 class         3    NA area_mn  1.67   1.67    4      4   
##  4 class         1    NA ca      72     72      33     33   
##  5 class         2    NA ca      23     23      47     47   
##  6 class         3    NA ca       5      5      20     20   
##  7 class         1    NA ed      29     29      58     58   
##  8 class         2    NA ed      43     43      86     86   
##  9 class         3    NA ed      16     16      40     40   
## 10 class         1    NA np       1      1      11      8   
## 11 class         2    NA np       5      3       9      3   
## 12 class         3    NA np       3      3       5      5   
## 13 landscape    NA    NA contag  41.8   41.8     8.78   8.78
## 14 landscape    NA    NA shdi     0.724  0.724   1.04   1.04

[Table 3.1

This will be your metrics table. You can print it in an intuitive order using this code:

#Paste table 3.1 code here
llmetrics <- list(esett_4, esett_8, psett_4, psett_8) #This creates a list of the data tables
metric_table <- do.call(rbind, allmetrics) #This binds the tables together
metric_table <- tidyr::spread(metric_table, layer, value) #This puts the table into wide format

metric_table[order(metric_table$level, metric_table$metric),] #This prints the table and orders it

## # A tibble: 14 × 8
##    level     class    id metric  esett4 esett8 psett4 psett8
##    <chr>     <int> <int> <chr>    <dbl>  <dbl>  <dbl>  <dbl>
##  1 class         1    NA area_mn 72     72       3      4.12
##  2 class         2    NA area_mn  4.6    7.67    5.22  15.7 
##  3 class         3    NA area_mn  1.67   1.67    4      4   
##  4 class         1    NA ca      72     72      33     33   
##  5 class         2    NA ca      23     23      47     47   
##  6 class         3    NA ca       5      5      20     20   
##  7 class         1    NA ed      29     29      58     58   
##  8 class         2    NA ed      43     43      86     86   
##  9 class         3    NA ed      16     16      40     40   
## 10 class         1    NA np       1      1      11      8   
## 11 class         2    NA np       5      3       9      3   
## 12 class         3    NA np       3      3       5      5   
## 13 landscape    NA    NA contag  41.8   41.8     8.78   8.78
## 14 landscape    NA    NA shdi     0.724  0.724   1.04   1.04

Question 3.2 What metrics changed between the two neighbourhood rules and what does this mean (4 pts)?

The metrics that changed between 4 vs 8 directions are mean patch size and number of patches.

Using 8 directions groups more neighbouring cells together, so patches become larger and fewer in number.

This means the landscape appears less fragmented under the 8-direction rule compared to the 4-direction rule.

Task 4: Real maps

Task 4.1 Create maps (show_landscape code) in R for each landscapes (2 pts).

## Import a raster and plot it. 
pre_fire <- raster("data/pre_fire.grd")
show_landscape(pre_fire, discrete = TRUE)

Your figure caption here

post_fire.grd <- raster("data/post_fire.grd")
show_landscape(post_fire.grd, discrete = TRUE)

Your figure caption here

## Import a raster and plot it. 

library(raster)
library(landscapemetrics)

land1 <- raster("data/pre_fire.grd")
land2 <- raster("data/post_fire.grd")

land1_mets <- calculate_lsm(
  land1,
  what = c(
    "lsm_l_contag",
    "lsm_l_pd",
    "lsm_l_ed",
    "lsm_l_lsi",
    "lsm_l_lpi",
    "lsm_l_pr"
  ),
  directions = 8
)

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

land2_mets <- calculate_lsm(
  land2,
  what = c(
    "lsm_l_contag",
    "lsm_l_pd",
    "lsm_l_ed",
    "lsm_l_lsi",
    "lsm_l_lpi",
    "lsm_l_pr"
  ),
  directions = 8
)

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

land1_mets$layer <- "pre"
land2_mets$layer <- "post"

Question 4.2 What is the extent and grain for the landscape (2 pts)?

The landscape extent is defined by its bounding coordinates where the total spatial area covered by the raster.

The grain is the resolution of the raster, represents the smallest unit of the landscape.

Task 5: Landscape metrics

Question 5.1 What do each of the metrics mean (include interpretations!)? (6 pts)

Contagion
Measures how clumped or connected patches are. High values mean patches are aggregated and connected; low values mean the landscape is fragmented and dispersed.

Patch Density
Measures the number of patches per area. Higher values indicate more fragmentation; lower values indicate a more continuous landscape.

Edge Density
Measures the amount of edge per area. Higher values mean more boundaries between patches; lower values mean fewer edges.

Landscape Shape Index
Measures how complex patch shapes are. Higher values indicate more irregular and complex shapes; lower values indicate simpler, more compact shapes.

Largest Patch Index
Shows the percentage of the landscape occupied by the largest patch. Higher values mean one patch dominates; lower values mean no single patch dominates.

Patch Richness
Measures the number of different patch types. Higher values indicate greater diversity; lower values indicate fewer types.

Task 6: Plotting Landscape metrics

Task 6.1 For each metric, graph using a bar graph (metric on the Y-axis, landscape on the X-axis); your plots will have two, one for each landscape (4 pts).

land1_mets$layer <- "pre"
land2_mets$layer <- "post"
# Task 6.1
library(ggplot2)

# Combine metric tables
allmetrics <- list(land1_mets, land2_mets)
metric_table <- do.call(rbind, allmetrics)

# Make sure layer is a factor in the right order
metric_table$layer <- factor(metric_table$layer, levels = c("pre", "post"))

# Plot
p <- ggplot(metric_table, aes(x = layer, y = value, fill = layer)) +
  geom_bar(stat = "identity") +
  facet_wrap(~metric, scales = "free") +
  scale_fill_manual(values = c("pre" = "blue", "post" = "red")) +
  xlab("Landscape") +
  ylab("Value") +
  ggtitle("Landscape level metrics for the pre- and post-fire landscapes") +
  theme_minimal()

print(p)

Question 6.2 Which of the landscapes is the most fragmented, and which appears least fragmented? How did you determine this? (3 pts)

The post-fire landscape is the most fragmented, while the pre-fire landscape is the least fragmented.

Higher patch density and edge density in the post-fire landscape. Lower contagion in the post-fire landscape – less connected, more dispersed patches

Question 6.3 : How would the correlation among landscape metrics influence your choice of what to report in an analysis that describes landscape pattern? (3 pts)

Many landscape metrics are correlated, they show similar information – choose a few key metrics to avoid repetition and keep the analysis clear.

Task 7 - Landscape Change Over Time

Task 7.1 Generate bar graphs for these class metrics (4 pts):

Patch Density
Edge Density
Landscape Shape Index
Largest Patch Index

Hint You can use the code from previous questions to calculate your metrics and put together your table.

#Task 7 code here
library(raster)
library(landscapemetrics)
library(dplyr)
library(tidyr)
library(ggplot2)

# Load rasters
land1 <- raster("data/pre_fire.grd")
land2 <- raster("data/post_fire.grd")

# Calculate class-level metrics
land1_mets <- calculate_lsm(
  land1,
  what = c(
    "lsm_c_pd",   # Patch Density
    "lsm_c_ed",   # Edge Density
    "lsm_c_lsi",  # Landscape Shape Index
    "lsm_c_lpi"   # Largest Patch Index
  ),
  directions = 8
)

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

land2_mets <- calculate_lsm(
  land2,
  what = c(
    "lsm_c_pd",
    "lsm_c_ed",
    "lsm_c_lsi",
    "lsm_c_lpi"
  ),
  directions = 8
)

## Warning: Please use 'check_landscape()' to ensure the input data is valid.

# Add layer labels
land1_mets$layer <- "pre"
land2_mets$layer <- "post"

# Combine tables
metric_table3 <- bind_rows(land1_mets, land2_mets)

# Order layer
metric_table3$layer <- factor(metric_table3$layer, levels = c("pre", "post"))

# Bar plot
p <- ggplot(
  data = metric_table3,
  mapping = aes(x = factor(class), y = value, fill = layer)
)

p +
  geom_bar(position = "dodge", stat = "identity") +
  facet_wrap(~metric, scales = "free") +
  xlab("Class") +
  ylab("Value") +
  scale_fill_manual(
    values = c("pre" = "blue", "post" = "red"),
    labels = c("Pre-fire", "Post-fire"),
    name = "Landscape"
  ) +
  ggtitle("Class-level metrics for pre- and post-fire landscapes") +
  theme_minimal()

Question 7.2 : Based on the metrics how would you describe the change in the structure of each land cover class over time? (8 pts)

Class 1–3 – These show higher PD and ED, meaning more patches and more edges, indicating increased fragmentation. LSI increases, showing more irregular patch shapes. This suggests these classes became more patchy and dispersed after the fire.
Class 4–5 – these also show increases in PD and LSI, indicating greater fragmentation and more complex shapes. Changes in LPI suggest reduced dominance of large continuous patches.
Class 6 – This class maintains a high LPI, meaning it still dominates the landscape, although slight increases in ED and PD suggest some fragmentation.
Class 7 – Shows little change across metrics, indicating it remains relatively stable and limited in extent.
The fire led to more fragmented, smaller, and more irregular patches across most classes, with reduced dominance of large patches except for class 6.

Copyright information:

You need to knit your handout including this copyright information, do not remove either image. Removing this information will result in a copyright infringement. We cannot grade submitted assignments with missing copyright information.