Comparison with FRAGSTATS

landscapemetrics re-implements landscape metrics as they are mostly described in the FRAGSTATS software (McGarigal et al. 2012). Therefore, we compared our results with the results of FRAGSTATS. 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. We list the metrics with deviations at the lowest level.

Unfortunatly, we do not have access to the source code of FRAGSTATS. Therefore, we are not able to finally explain the present differences between the results, nevertheless, we try to guess the most likely reasons.

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 after the 4th decimal place in FRAGSTATS. Contrastingly, we do not round the numbers. Naturally, this can lead to small deviations between the results.

There are quite a few metrics on class- and landscape-level that summarise 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 standard deviation and the coefficent of variation. landscapemetrics uses base R functions for that, so we should assume that the calculation of such indices is correct.

In the following, we are comparing the cv 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 actual 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
fragstats_calculated <- fragstats_patch_landscape %>%
dplyr::group_by(TYPE) %>%
dplyr::summarise(cv = cv(AREA)) %>%
purrr::set_names("class", "fragstats_calculated")

# Output of FRAGSTATS
fragstats_output <- fragstats_class_landscape %>%
dplyr::select(TYPE, AREA_CV) %>%
purrr::set_names("class", "fragstats_output")

# Output of landscapemetrics
landscapemetrics_output <- lsm_c_area_cv(landscape) %>%
dplyr::select(class, value) %>%
purrr::set_names("class", "landscapemetrics")

fragstats <- dplyr::full_join(x = fragstats_output, y = fragstats_calculated,
by = "class")

cv_full <- dplyr::full_join(x = fragstats, y = landscapemetrics_output,
by = "class")
knitr::kable(cbind(cv_full))
class fragstats_output fragstats_calculated landscapemetrics
1 228.5906 242.4570 242.4570
3 162.6136 187.7700 187.7700
2 159.0109 165.0134 165.0134

As for the cv, the results for the sd are similiar. The result calculated from all patch areas of FRAGSTATS is identical to the result of landscapemetrics, but not the actual result of FRAGSTATS.

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

# Output of FRAGSTATS
fragstats_output <- fragstats_class_landscape %>%
dplyr::select(TYPE, AREA_SD) %>%
purrr::set_names("class", "fragstats_output")

# Output of landscapemetrics
landscapemetrics_output <- lsm_c_area_sd(landscape) %>%
dplyr::select(class, value) %>%
purrr::set_names("class", "landscapemetrics")

fragstats <- dplyr::full_join(x = fragstats_output, y = fragstats_calculated,
by = "class")

sd_full <- dplyr::full_join(x = fragstats, y = landscapemetrics_output,
by = "class")
knitr::kable(cbind(cv_full))
class fragstats_output fragstats_calculated landscapemetrics
1 228.5906 242.4570 242.4570
3 162.6136 187.7700 187.7700
2 159.0109 165.0134 165.0134

Specific differences

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
fragstats_area <- fragstats_patch_landscape %>%
dplyr::select(PID, AREA) %>%
purrr::set_names("id", "fragstats_area")

landscapemetrics_area <- lsm_p_area(landscape) %>%
dplyr::select(id, value) %>%
purrr::set_names("id", "landscapemetrics_area")

# Calculate GYRATE
fragstats_gyrate <- fragstats_patch_landscape %>%
dplyr::select(PID, GYRATE) %>%
purrr::set_names("id", "fragstats_gyrate")

landscapemetrics_gyrate <- lsm_p_gyrate(landscape) %>%
dplyr::select(id, value) %>%
purrr::set_names("id", "landscapemetrics_gyrate")

fragstats <- dplyr::full_join(x = fragstats_area, y = fragstats_gyrate,
by = "id") %>%
dplyr::filter(fragstats_area == 0.0001)

landscapemetrics <- dplyr::full_join(x = landscapemetrics_area,
y = landscapemetrics_gyrate, by = "id") %>%
dplyr::filter(landscapemetrics_area == 0.0001)
knitr::kable(cbind(fragstats, landscapemetrics))
id fragstats_area fragstats_gyrate id landscapemetrics_area landscapemetrics_gyrate
1 1e-04 0.5 1 1e-04 0
19 1e-04 0.5 4 1e-04 0
23 1e-04 0.5 5 1e-04 0
25 1e-04 0.5 9 1e-04 0
10 1e-04 0.5 15 1e-04 0
20 1e-04 0.5 17 1e-04 0

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

Error propagation (for metrics based on GYRATE metric)
Class level
• GYRATE_CV (lsm_c_gyrate_cv)
• GYRATE_MN (lsm_c_gyrate_mn)
• GYRATE_SD (lsm_c_gyrate_sd)
Landscape level
• GYRATE_CV (lsm_l_gyrate_cv)
• GYRATE_MN (lsm_l_gyrate_mn)
• GYRATE_SD (lsm_l_gyrate_sd)

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
fragstats <- fragstats_patch_landscape %>%
dplyr::select(PID, AREA, PERIM, PARA) %>%
purrr::set_names("id", "area", "perim", "para") %>%
dplyr::mutate(para_calculated_ha = perim / area,
para_calculated_m = perim / (area * 10000)) %>%
dplyr::arrange(area)

# Output of landscapemetrics
area_landscapmetrics <- lsm_p_area(landscape) %>%
dplyr::select(id, value) %>%
purrr::set_names("id", "area")

perim_landscapmetrics <- lsm_p_perim(landscape) %>%
dplyr::select(id, value) %>%
purrr::set_names("id", "perim")

para_landscapemetrics <- lsm_p_para(landscape) %>%
dplyr::select(id, value) %>%
purrr::set_names("id", "para")

landscapemetrics <- dplyr::full_join(x = area_landscapmetrics,
y = perim_landscapmetrics,
by = "id") %>%
dplyr::full_join(para_landscapemetrics, by = "id") %>%
dplyr::mutate(para_calculated_ha = perim / area,
para_calculated_m = perim / (area * 10000)) %>%
dplyr::arrange(area)
knitr::kable(cbind(fragstats, landscapemetrics))
id area perim para para_calculated_ha para_calculated_m id area perim para para_calculated_ha para_calculated_m
1 0.0001 4 40000.000 40000.000 4.0000000 1 0.0001 4 4.0000000 40000.000 4.0000000
19 0.0001 4 40000.000 40000.000 4.0000000 4 0.0001 4 4.0000000 40000.000 4.0000000
23 0.0001 4 40000.000 40000.000 4.0000000 5 0.0001 4 4.0000000 40000.000 4.0000000
25 0.0001 4 40000.000 40000.000 4.0000000 9 0.0001 4 4.0000000 40000.000 4.0000000
10 0.0001 4 40000.000 40000.000 4.0000000 15 0.0001 4 4.0000000 40000.000 4.0000000
20 0.0001 4 40000.000 40000.000 4.0000000 17 0.0001 4 4.0000000 40000.000 4.0000000
9 0.0002 6 30000.000 30000.000 3.0000000 11 0.0002 6 3.0000000 30000.000 3.0000000
11 0.0002 6 30000.000 30000.000 3.0000000 12 0.0002 8 4.0000000 40000.000 4.0000000
15 0.0002 8 40000.000 40000.000 4.0000000 14 0.0002 6 3.0000000 30000.000 3.0000000
26 0.0003 8 26666.667 26666.667 2.6666667 7 0.0003 8 2.6666667 26666.667 2.6666667
27 0.0003 8 26666.667 26666.667 2.6666667 18 0.0003 8 2.6666667 26666.667 2.6666667
12 0.0003 8 26666.667 26666.667 2.6666667 19 0.0003 8 2.6666667 26666.667 2.6666667
13 0.0003 8 26666.667 26666.667 2.6666667 23 0.0003 8 2.6666667 26666.667 2.6666667
18 0.0003 8 26666.667 26666.667 2.6666667 26 0.0003 8 2.6666667 26666.667 2.6666667
14 0.0004 8 20000.000 20000.000 2.0000000 21 0.0004 8 2.0000000 20000.000 2.0000000
8 0.0005 12 24000.000 24000.000 2.4000000 2 0.0005 12 2.4000000 24000.000 2.4000000
22 0.0005 10 20000.000 20000.000 2.0000000 8 0.0005 10 2.0000000 20000.000 2.0000000
24 0.0007 12 17142.857 17142.857 1.7142857 22 0.0007 12 1.7142857 17142.857 1.7142857
2 0.0009 16 17777.778 17777.778 1.7777778 25 0.0009 16 1.7777778 17777.778 1.7777778
21 0.0010 16 16000.000 16000.000 1.6000000 27 0.0010 16 1.6000000 16000.000 1.6000000
16 0.0014 20 14285.714 14285.714 1.4285714 6 0.0014 20 1.4285714 14285.714 1.4285714
7 0.0024 32 13333.333 13333.333 1.3333333 16 0.0024 32 1.3333333 13333.333 1.3333333
3 0.0035 38 10857.143 10857.143 1.0857143 10 0.0035 38 1.0857143 10857.143 1.0857143
6 0.0057 40 7017.544 7017.544 0.7017544 20 0.0057 40 0.7017544 7017.544 0.7017544
17 0.0098 82 8367.347 8367.347 0.8367347 13 0.0098 82 0.8367347 8367.347 0.8367347
5 0.0148 130 8783.784 8783.784 0.8783784 3 0.0148 130 0.8783784 8783.784 0.8783784
4 0.0457 348 7614.880 7614.880 0.7614880 24 0.0457 348 0.7614880 7614.880 0.7614880
Error propagation (for metrics based on PARA metric)
Class level
• PARA_MN (lsm_c_para_mn)
• PARA_SD (lsm_c_para_sd)

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.

# Output of FRAGSTATS
fragstats <- fragstats_patch_landscape %>%
dplyr::select(PID, AREA, CIRCLE) %>%
purrr::set_names("id_fs", "area_fs", "circle_fs") %>%
dplyr::filter(area_fs == 0.0001)

# Output of landscapemetrics
area_landscapmetrics <- lsm_p_area(landscape) %>%
dplyr::select(id, value) %>%
purrr::set_names("id_lsm", "area_lsm")

circle_landscapmetrics <- lsm_p_circle(landscape) %>%
dplyr::select(id, value) %>%
purrr::set_names("id_lsm", "circle_lsm")

landscapemetrics <- dplyr::full_join(x = area_landscapmetrics, y = circle_landscapmetrics,
by = "id_lsm") %>%
dplyr::filter(area_lsm == 0.0001)
knitr::kable(cbind(fragstats, landscapemetrics))
id_fs area_fs circle_fs id_lsm area_lsm circle_lsm
1 1e-04 0 1 1e-04 0.3633802
19 1e-04 0 4 1e-04 0.3633802
23 1e-04 0 5 1e-04 0.3633802
25 1e-04 0 9 1e-04 0.3633802
10 1e-04 0 15 1e-04 0.3633802
20 1e-04 0 17 1e-04 0.3633802

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 raster
ras <- raster::raster(mat, xmn = 0, xmx = 13, ymn = 0, ymx = 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_sp = 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 diamter /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_sp = 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.

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

Patch metrics

To get all availabel 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.

library(SDMTools) # remotes::install_version("SDMTools", version = "1.1-221.2")
##
## Attaching package: 'SDMTools'
## The following object is masked from 'package:raster':
##
##     distance
# binarize every class in the landscape and calculate patch metrics
sdmtools_result <- lapply(raster::unique(landscape), FUN = function(x){
tmp_land <- landscape
raster::values(tmp_land)[raster::values(tmp_land) != x] <- NA
raster::values(tmp_land)[raster::values(tmp_land) == x] <- 1
ccl <- SDMTools::ConnCompLabel(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")
## Warning: Please use 'check_landscape()' to ensure the input data is valid.
knitr::kable(bind_rows(sdmtools_result))
Class patchID n.cell n.core.cell n.edges.perimeter n.edges.internal area core.area perimeter perim.area.ratio shape.index frac.dim.index core.area.index
1 1 1 0 4 0 1 0 4 4.0000000 1.000000 NaN 0.0000000
1 2 1 0 4 0 1 0 4 4.0000000 1.000000 NaN 0.0000000
1 3 148 39 130 462 148 39 130 0.8783784 2.600000 1.393273 0.2635135
1 4 5 0 12 8 5 0 12 2.4000000 1.200000 1.365212 0.0000000
1 5 1 0 4 0 1 0 4 4.0000000 1.000000 NaN 0.0000000
1 6 14 0 20 36 14 0 20 1.4285714 1.250000 1.219707 0.0000000
1 7 3 0 8 4 3 0 8 2.6666667 1.000000 1.261859 0.0000000
1 8 5 0 10 10 5 0 10 2.0000000 1.000000 1.138647 0.0000000
1 9 1 0 4 0 1 0 4 4.0000000 1.000000 NaN 0.0000000
2 1 2 0 6 2 2 0 6 3.0000000 1.000000 1.169925 0.0000000
2 2 35 5 38 102 35 5 38 1.0857143 1.583333 1.266425 0.1428571
2 3 98 31 82 310 98 31 82 0.8367347 2.050000 1.317534 0.3163265
2 4 2 0 8 0 2 0 8 4.0000000 1.333333 2.000000 0.0000000
2 5 1 0 4 0 1 0 4 4.0000000 1.000000 NaN 0.0000000
2 6 2 0 6 2 2 0 6 3.0000000 1.000000 1.169925 0.0000000
2 7 1 0 4 0 1 0 4 4.0000000 1.000000 NaN 0.0000000
2 8 24 0 32 64 24 0 32 1.3333333 1.600000 1.308626 0.0000000
2 9 3 0 8 4 3 0 8 2.6666667 1.000000 1.261859 0.0000000
2 10 3 0 8 4 3 0 8 2.6666667 1.000000 1.261859 0.0000000
2 11 57 22 40 188 57 22 40 0.7017544 1.250000 1.139033 0.3859649
2 12 4 0 8 8 4 0 8 2.0000000 1.000000 1.000000 0.0000000
2 13 7 0 12 16 7 0 12 1.7142857 1.000000 1.129150 0.0000000
2 14 3 0 8 4 3 0 8 2.6666667 1.000000 1.261859 0.0000000
3 1 457 158 348 1480 457 158 348 0.7614880 4.046512 1.458331 0.3457330
3 2 9 0 16 20 9 0 16 1.7777778 1.333333 1.261859 0.0000000
3 3 3 0 8 4 3 0 8 2.6666667 1.000000 1.261859 0.0000000
3 4 10 0 16 24 10 0 16 1.6000000 1.142857 1.204120 0.0000000

References

• McGarigal, K., SA Cushman, and E Ene. 2012. FRAGSTATS v4: Spatial Pattern Analysis Program for Categorical and Continuous Maps. Computer software program produced by the authors at the University of Massachusetts, Amherst. Available at the following website: http://www.umass.edu/landeco/research/fragstats/fragstats.html

• 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