In landscapeecology/landscapemetrics: Landscape Metrics for Categorical Map Patterns

knitr::opts_chunk$set(message = FALSE, warning = FALSE, include = FALSE)

library(landscapemetrics)

library(dplyr)
library(purrr)
library(raster)
library(SDMTools) # remotes::install_version("SDMTools", version = "1.1-221.2")
library(stringr)
library(terra)
library(tidyr)

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

fs_patch_landscape <- landscapemetrics:::internal_data$fs_data$patch |> 
    dplyr::filter(LID == "landscape") |> 
    dplyr::select(-LID) |> 
    dplyr::mutate(TYPE = stringr::str_remove_all(TYPE, pattern = " "),
                  TYPE = as.integer(stringr::str_remove(TYPE, pattern = "cls_")))

fs_class_landscape <- landscapemetrics:::internal_data$fs_data$class |> 
    dplyr::filter(LID == "landscape") |> 
    dplyr::select(-LID) |> 
    dplyr::mutate(TYPE = stringr::str_remove_all(TYPE, pattern = " "),
                  TYPE = as.integer(stringr::str_remove(TYPE, pattern = "cls_")))

Comparison with FRAGSTATS

landscapemetrics re-implements landscape metrics as they are mostly described in the FRAGSTATS software (McGarigal et al. 2023). Therefore, we compare our results with the results of FRAGSTATS. On patch-level we use a correlation approach testing values calculated by landscapemetrics and FRAGSTATS have at least a correlation of p>=0.975 because we are not able to exactly match patches. On class- and landscape-level we can match individual values and test if no deviation is larger than 5%. In the process, we recognized a few differences between the results. Some metrics in FRAGSTATS are interdependent across scales. Thus, if there is a deviation at the patch level, it propagates through the class- and landscape-level. In this document, we list the metrics with deviations at the lowest level only. Unfortunately, we do not have access to the source code of FRAGSTATS. Therefore, we are not able to finally explain all present differences of results.

General differences

Firstly, the patch ID is ordered in a different way, due to technical reasons (how connected patches are specified). Therefore, one has to pay attention comparing the results on patch level for FRAGSTATS and landscapemetrics.

All double precision floating point numbers are rounded in FRAGSTATS. Contrastingly, we do not round the numbers. Naturally, this can lead to small deviations between the results. However, these small differences should have no influence on any interpretation of the results.

There are quite a few metrics on class- and landscape-level that summarize patch level metrics (e.g., the mean, standard deviation (sd) or coefficient of variation (cv) of all values belonging to class i). While the results are identical for single patches and the mean of all patches, there are some slight differences between lanscapemetrics and FRAGSTATS for the sd and the cv. landscapemetrics uses base R functions and standard statistical definitions to summarize the patch level metrics. Again, difference are only very minor and should not affect interpretation of the results. In the following, we are comparing the cv and sd for the patch area.

We are including the cv calculated from all patch areas and the actual output of FRAGSTATS as well as the output of landscapemetrics. Interestingly, the cv calculated from all patches of FRAGSTATS is identical to the cv of landscapemetrics, but the already summarized result of FRAGSTATS is slightly different.

# function to calculate coefficient of variation
cv <- function(x) {
    (sd(x) / mean(x)) * 100
}

# CV calculated from patch values of FRAGSTATS 
fs_calculated <- fs_patch_landscape |>
    dplyr::group_by(TYPE) |>
    dplyr::summarise(cv = cv(AREA)) |>
    purrr::set_names("class", "fs_calc") 

# Output of FRAGSTATS
fs_output <- fs_class_landscape |> 
    dplyr::select(TYPE, AREA_CV) |>
    purrr::set_names("class", "fs_output")

# Output of landscapemetrics
lsm_calc <- lsm_p_area(landscape) |>
    dplyr::group_by(class) |>
    dplyr::summarise(value = cv(value)) |> 
    purrr::set_names("class", "lsm_calc")

lsm_output <- lsm_c_area_cv(landscape) |>
    dplyr::select(class, value) |>
    purrr::set_names("class", "lsm_fun")

fragstats <- dplyr::full_join(x = fs_output, y = fs_calculated, 
                              by = "class") 

lsm <- dplyr::full_join(x = lsm_calc, y = lsm_output, 
                        by = "class") 

cv_full <- dplyr::full_join(x = fragstats, y = lsm,
                            by = "class")

Click for details

knitr::kable(cv_full)

As for the cv, there are also minor differences for the sd. The result calculated from all patch areas of FRAGSTATS is identical to the result of landscapemetrics, but not the summarized result of FRAGSTATS.

# SD calculated from patch values of FRAGSTATS
fs_calculated <- fs_patch_landscape |>
  dplyr::group_by(TYPE) |>
  dplyr::summarise(sd = sd(AREA)) |>
  purrr::set_names("class", "fs_calculated") 

# Output of FRAGSTATS
fs_output <- fs_class_landscape |>
  dplyr::select(TYPE, AREA_SD) |>
  purrr::set_names("class", "fs_output")

# Output of landscapemetrics
lsm_calc <- lsm_p_area(landscape) |>
    dplyr::group_by(class) |> 
    dplyr::summarise(value = sd(value)) |> 
    purrr::set_names("class", "lsm_calculated")

lsm_output <- lsm_c_area_sd(landscape) |>
  dplyr::select(class, value) |>
  purrr::set_names("class", "lsm_output")

fragstats <- dplyr::full_join(x = fs_output, y = fs_calculated,
                              by = "class")

lsm <- dplyr::full_join(x = lsm_calc, y = lsm_output, 
                        by = "class") 

sd_full <- dplyr::full_join(x = fragstats, y = lsm,
                            by = "class")

Click for details

knitr::kable(sd_full)

Specific differences

All core-related metrics

Different definitions of the core area of a patch result in differences for all metrics that somehow related on the core area or number of disjunct core areas. landscapemetrics uses a cell-centered definition of the core area, namely "a cell is defined as core area if the cell has no neighbour with a different value than itself (rook's case)." Contrastingly, FRAGSTATS uses "[...] a method involving the use of a variably-sized masked placed on cells on the perimeter of a patch [...]". Especially for patches with rather complicated shapes this seem to result in different core-metric values.

GYRATE metric

According to FRAGSTATS the radius of gyration for a patch consisting of only a single cell should equal GYRATE = 0.

[...] GYRATE = 0 when the patch consists of a single cell [...]

However, for patches containing a single cell FRAGSTATS returns a value of GYRATE = 0.5. In the following table, patches with an area of area = 0.0001 consist of only one cell.

# Calculate patch area
# Calculate GYRATE
fs <- fs_patch_landscape |>
    dplyr::select(PID, AREA, GYRATE) |>
    purrr::set_names("fs_id", "fs_area", "fs_gyrate") |> 
    dplyr::filter(fs_area == 0.0001)

lsm <- calculate_lsm(landscape, what = c("lsm_p_area", "lsm_p_gyrate")) |>
    tidyr::pivot_wider(names_from = "metric", values_from = "value") |> 
    dplyr::select(id, area, gyrate) |> 
    purrr::set_names("lsm_id", "lsm_area", "lsm_gyrate") |> 
    dplyr::filter(lsm_area == 0.0001)

Click for details

knitr::kable(cbind(fs, lsm))

Additionally, we recognized small differences for all other patches as well. However, we could not find an explanation for this difference, yet.

PARA metric

The documentation of FRAGSTATS defines the perimeter-area ratio the following:

[...] PARA equals the ratio of the patch perimeter (m) to area (m2). [...]

Contrastingly, the output of FRAGSTATS gives the result as the ratio of the patch perimeter in meters to area in hectares.

We implemented PARA as documented in the FRAGSTATS manual using square meters. Nevertheless, the differences between the softwares are only based on different units, as shown by converting the FRAGSTATS output to meters per square meters.

# Output of FRAGSTATS
fs <- fs_patch_landscape |>
    dplyr::select(PID, AREA, PERIM, PARA) |>
    purrr::set_names("fs_id", "fs_area", "fs_perim", "fs_para") |>
    dplyr::arrange(fs_area)

lsm <- calculate_lsm(landscape, what = c("lsm_p_area", "lsm_p_perim", "lsm_p_para")) |> 
    dplyr::select(-layer, -level, -class) |> 
    tidyr::pivot_wider(names_from = "metric", values_from = "value") |> 
    dplyr::mutate(lsm_para_ha = perim / area) |> 
    dplyr::select(id, area, perim, para, lsm_para_ha) |> 
    purrr::set_names("lsm_id", "lsm_area", "lsm_perim", "lsm_para", "lsm_para_ha") |>
    dplyr::arrange(lsm_area)

Click for details

knitr::kable(cbind(fs, lsm))

SHAPE metric

We are following the definition provided in the latest version of the FRAGSTATS manual, in which the shape index is defined as the perimeter (m) divided by the square root of the area (m2) of patch ij, adjusted by a constant to adjust for a square standard. Yet, for some patches there are minor differences between landscapemetrics and FRAGSTATS.

CIRCLE metric

According to FRAGSTATS, for patches with only one cell CIRCLE = 0.

[...] CIRCLE = 0 for one cell patches. [...]

Nevertheless, because also patches with only one cell have a dimension in the raster context, we decided to calculate the CIRCLE metric for such patches.

# Calculate patch area
# Calculate GYRATE
fs <- fs_patch_landscape |>
    dplyr::select(PID, AREA, CIRCLE) |>
    purrr::set_names("fs_id", "fs_area", "fs_circle") |> 
    dplyr::filter(fs_area == 0.0001)

lsm <- calculate_lsm(landscape, what = c("lsm_p_area", "lsm_p_circle")) |>
    tidyr::pivot_wider(names_from = "metric", values_from = "value") |> 
    dplyr::select(id, area, circle) |> 
    purrr::set_names("lsm_id", "lsm_area", "lsm_circle") |> 
    dplyr::filter(lsm_area == 0.0001)

Click for details

knitr::kable(cbind(fs, lsm))

It seems like FRAGSTATS uses the largest distance between the corner points of the cells to calculate the diameter of the circumscribing circle. However, this does not necessarily result in the correct circumscribing circle. While this approach works well for more compact patch shapes, there are some cases in which the approach fails. One example are T-shaped patches. Contrastingly to FRAGSTATS, our algorithm calculates the true circumscribing circle also for such patches.

# create matrix
mat <- matrix(data = NA, nrow = 13, ncol = 13)

# create T-shaped patches
mat[4:9, 7] <- 1
mat[4, 4:10] <- 1

# convert to terra
ras <- terra::rast(mat, extent = c(0, 13, 0, 13))

# get circumscribing circle 
circle_lsm <- get_circumscribingcircle(ras)

# construct circle using diameter / 2
circle_lsm <- construct_buffer(coords = matrix(c(circle_lsm$center_x, 
                                                 circle_lsm$center_y),
                                               nrow = 1, ncol = 2),
                               shape = "circle", size = circle_lsm$value / 2,
                               return_vec  = FALSE)

# calculate max distance between corner points of cells
circle_max_dist <- dist(matrix(data = c(3, 10, 7, 4), byrow = TRUE, nrow = 2, ncol = 2))

# construct circle using diameter /2
circle_max_dist <- construct_buffer(coords = matrix(c(6.5, 7.5), nrow = 1, ncol = 2),
                                    shape = "circle", size = circle_max_dist / 2,
                                    return_vec = FALSE)

# plot results
plot(ras, col = "#D5B528", legend = FALSE)
polygon(circle_lsm, border = "#84A98E", lwd = 2.5) # green circle
polygon(circle_max_dist, border = "#922418", lwd = 2.5) # red circle

Comparison with SDMTools

SDMTools (VanDerWal et al. 2014) (still available, but apparently not longer maintained) offers landscape metrics on patch and class level. However, it does not return the same results as FRAGSTATS. The main reason for this are different standard defaults (e.g., SDMTools always considers the global landscape boundary) and that SDMTools returns results in map units and not in m^2/hectar, as FRAGSTATS/landscapemetrics. This also explains differences between our package and SDMTools.

Jemma Stachelek was so nice to remind us of these issues and provided the comparison.

Patch metrics

To get all available metrics on, e.g., patch level with SDMTools, you have to make a binary landscape for every class in your landscape, perform connected components labelling on it and then calculate the patch metrics.

# binarize every class in the landscape and calculate patch metrics
sdmtools_result <- lapply(terra::unique(landscape), FUN = function(x){
  tmp_land <- landscape
  terra::values(tmp_land)[terra::values(tmp_land) != x] <- NA
  terra::values(tmp_land)[terra::values(tmp_land) == x] <- 1
  ccl <- SDMTools::ConnCompLabel(raster::raster(tmp_land))
  PatchStat(ccl)
})

# combine to one data frame
sdmtools_result <- dplyr::bind_rows(sdmtools_result, .id = "Class")

landscapemetrics offers for such tasks the function get_patches() and for the metrics itself all of that is done internally. To get all metrics on patch level with landscapemetrics you could for example do:

patch_metrics <- calculate_lsm(landscape, what = "patch")

Click for details

knitr::kable(sdmtools_result)

References

• McGarigal K., SA Cushman, and E Ene. 2023. FRAGSTATS v4: Spatial Pattern Analysis Program for Categorical Maps. Computer software program produced by the authors; available at the following web site: https://www.fragstats.org
• Jeremy VanDerWal, Lorena Falconi, Stephanie Januchowski, Luke Shoo and Collin Storlie (2014). SDMTools: Species Distribution Modelling Tools: Tools for processing data associated with species distribution modelling exercises. R package version 1.1-221. https://CRAN.R-project.org/package=SDMTools

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