In landscapeecology/landscapemetrics: Landscape Metrics for Categorical Map Patterns

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(landscapemetrics)
library(terra)
library(bench)
library(dplyr)
library(purrr)

# internal data needs to be read
landscape <- terra::rast(landscapemetrics::landscape)

Visualisization functions

Visualizing landscapes

We can use the plot function from terra to have look at our landscape.

# Plot landscape
plot(landscape) 

This is how we typically inspect our landscape, but which also makes it quite hard to relate to the landscape metrics we are interested in. This why we show in the following how to dissect this landscape visually into the components, that drive the calculation of landscape metrics.

Visualizing patches

To visualize patches in a landscape and encode each patch with an ID that can be used to compare a landscape metric with the actual landscape you can use the auxiliary visualisation function show_patches():

# Plot landscape + landscape with labeled patches
show_patches(landscape) 

You can also plot all patches of each class grouped.

# show patches of all classes
show_patches(landscape, class = "all", labels = FALSE)

To show only the core area, there is the visualization function show_cores. The arguments are similar to show_patches():

# show core area of class 1 and 3
show_cores(landscape, class = c(1, 2), labels = FALSE)

Lastly, you can also "fill" the colours of each patch according to its value of a certain patch level metric, e.g., the patch area, using show_lsm(). You can chose if the label should be the patch id or the actual value of the landscape metric (label_lsm = TRUE/FALSE).

# fill patch according to area
show_lsm(landscape, what = "lsm_p_area", class = "global", label_lsm = TRUE)

To get the result as a SpatRaster, there is spatialize_lsm().

spatialize_lsm(landscape, what = "lsm_p_area")

Show correlation

Selecting meaningful landscape metrics for your field of research is difficult, as many landscape metrics are very abstract and the common approach is often simply to calculate as many as possible.

To select at the least that ones for your landscape and research question that are not highly correlated, you can use the function show_correlation() to get a first insight into the correlation of the metrics you calculated:

metrics <- calculate_lsm(landscape, what = "patch")
show_correlation(metrics, method = "pearson")

Building blocks

Get patches

landscapemetrics makes internally heavy use of an connected labeling algorithm and exports an re-implementation of this algorithm (get_patches()). This function return a list, where each list entry includes all patches of the corresponding class. The patches are labeled from 1...to n.

# get a list of all patches for each class
get_patches(landscape)

Get adjacencies

Adjacencies are a central part for landscape metrics, so calculating them quick and in a flexible way is key for, e.g., developing new metrics. Hence, landscapemetrics exports a function that can calculate adjacencies in any number if directions when provided with a binary matrix (NA / 1 - NA are cells that would be left out for looking at adjacencies).

# calculate full adjacency matrix
get_adjacencies(landscape, neighbourhood = 4)

# count diagonal neighbour adjacencies
diagonal_matrix <- matrix(c(1,  NA,  1,
                            NA,  0, NA,
                            1,  NA,  1), 3, 3, byrow = TRUE)
get_adjacencies(landscape, diagonal_matrix)

# equivalent with the terra package:
adj_terra <- function(x){
    adjacencies <- terra::adjacent(x, 1:terra::ncell(x), "rook", pairs=TRUE)
    table(terra::values(x, mat = FALSE)[adjacencies[,1]],
          terra::values(x, mat = FALSE)[adjacencies[,2]])
}

# compare the two implementations
bench::mark(
    get_adjacencies(landscape, neighbourhood = 4),
    adj_terra(landscape),
    iterations = 250, 
    check = FALSE
)

adj_terra(landscape) == get_adjacencies(landscape, 4)[[1]]

Get nearest neighbour

landscapemetrics implements a memory efficient and quite fast way to calculate the nearest neighbour between classes in a raster (or matrix).

# run connected labeling for raster
patches <- get_patches(landscape, class = 1)

# calculate the minimum distance between patches in a landscape
min_dist <- get_nearestneighbour(patches$layer_1$class_1)

# create a function that would do the same with the raster package
nn_terra <- function(patches) {

    np_class <- terra::values(patches[[1]][[1]]) |>
        na.omit() |>
        unique() |>
        length()

    points_class <- terra::as.data.frame(patches[[1]][[1]], xy = TRUE)

    minimum_distance <- seq_len(np_class) |>
        purrr::map_dbl(function(patch_ij) {

            patch_focal <- dplyr::filter(points_class, lyr.1 == patch_ij) |> 
                dplyr::select(x, y) |> 
                as.matrix(ncol = 2)

            patch_others <- dplyr::filter(points_class, lyr.1 != patch_ij) |> 
                dplyr::select(x, y) |> 
                as.matrix(ncol = 2)

            minimum_distance <- terra::distance(patch_focal, patch_others,
                                                lonlat = FALSE) |>
                min()
        })

    data.frame(id = unique(sort(points_class$lyr.1)), distance = minimum_distance)

}


# compare the two implementations
bench::mark(
    get_nearestneighbour(patches$layer_1$class_1)[, 2:3],
    nn_terra(patches$layer_1$class_1),
    iterations = 250, check = FALSE
)

# check if results are identical
get_nearestneighbour(patches$layer_1$class_1)[, 2:3] == nn_terra(patches$layer_1$class_1)

Get circumscribing circle

To get the smallest circumscribing circle that includes all cells of the patch, simply run get_circumscribingcircle(). The result returns the diameter for each circle that includes all cells of each patch. This includes not only the cell centers but the whole cells using the cells corners.

# get all patches of class 1
class_1 <- get_patches(landscape, class = 1)

# get smallest circumscribing circle for each patch
circle <- get_circumscribingcircle(class_1$layer_1$class_1)
circle

landscapeecology/landscapemetrics documentation built on April 7, 2024, 11:11 p.m.