--- title: "Calculating spatial metrics for IUCN RLE assessments" author: "Aniko Toth, José R. Ferrer-Paris, Calvin K.F. Lee, Nicholas J. Murray" date: "2026-01-27" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Calculating spatial metrics for IUCN RLE assessments} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: markdown: wrap: 72 --- ## 1. Introduction In this vignette, we provide an example of using the tools within `redlistr` to assess the spatial criteria of the [IUCN Red List of Ecosystems](https://iucnrle.org/). These criteria assess the change in the extent of an ecosystem over time (Criterion A) and properties of the geographic distribution size (Criterion B). Both of these criteria require the use of maps of ecosystem distributions. We begin by loading the packages we need for data import and for calculation of the spatial metrics. Please ensure that these are already installed. ```{r Loading packages, message=FALSE} library(terra) library(redlistr) ``` ## 2. Importing data We first import the spatial data we want to analyse. In this vignette, we will be using two example distributions included with the package. They are mangrove distributions from the northern regions of Western Port Bay, Victoria, Australia. The first distribution is from Giri et al., ([2011](https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1466-8238.2010.00584.x)), and represents mangrove distribution in 2000. The second distribution is a classification map generated using the `XGBoost` package and data from Landsat 8, and represents mangrove distribution in 2017. Both of these geographic distributions are stored as rasters, and a raster value of 1 denotes presence. We use function `rast` from the `terra` package to read these rasters. This function requires a valid path to find the files with the raster data. In this case both files are included in the `redlistr` package, so we use the `system.file` function to find the path to those files. More information about formats and functions of the `terra` package can be found [here](https://rspatial.github.io/terra/). ```{r Loading our example distributions} mangrove.2000 <- rast(system.file("extdata", "example_distribution_2000.tif", package = "redlistr")) mangrove.2017 <- rast(system.file("extdata", "example_distribution_2017.tif", package = "redlistr")) ``` ### Alternative formats Polygon or any shapefile data may be read in using `sf::st_read()` and point data in table format can be read in using `read.csv()` followed by `sf::st_as_sf()` to convert points into a spatial object. See function documentation in the sf package for more detailed instructions. ### Plotting out data We plot the data to view the two mangrove distributions to get a general idea of how the distribution might have changed over time. ```{r Plotting the two rasters, fig.show='hold', fig.width=7, fig.height=7, fig.alt='Overlay of two mangrove distribution maps'} plot(mangrove.2000, col = "grey30", legend = FALSE, main = "Mangrove Distribution") plot(mangrove.2017, add = T, col = "darkorange", legend = FALSE) ``` We can see that there has been some change in mangrove cover, although the losses appear to be relatively minor. ## 3. Basic information for the data We now use some of features of `redlistr` to get some basic information on the ecosystem distribution. However, before proceeding, we check their resolution (grain size), extent, and spatial properties. Specifically, the coordinate reference system (CRS) of the datasets must be measured in metres, and must be consistent between layers you wish to compare. ```{r Checking CRS} is.lonlat(mangrove.2000) is.lonlat(mangrove.2017) # If TRUE, they must be transformed to a projected coordinate system crs(mangrove.2000) == crs(mangrove.2017) ``` Reprojection can be done by using the `terra::project` function. Be sure to choose a coordinate reference system that is appropriate to your study system. We start by calculating the area of the ecosystem. Area is calculated here using the pixel count method, where the number of pixels classified as the ecosystem of interest is multiplied by the area of each pixel, giving the area of the ecosystem distribution. The area calculated is provided in square kilometres. ```{r Calculate area of rasters, eval = F} a.2000 <- getArea(mangrove.2000) a.2000 a.2017 <- getArea(mangrove.2017) a.2017 ``` The results: at time 1 (2000), the ecosystem was `r a.2000$area` km2, and at time 2, (2017) the ecosystem was `r a.2017$area` km2. ### Rasters with multiple ecosystems or layers The `getArea` function is capable of processing rasters with multiple ecosystems and layers. For example you can also calculate areas for both layers at once by combining them into a raster stack: ```{r Calculate area of multiple rasters, eval = F} a <- getArea(c(mangrove.2000, mangrove.2017)) ``` ## 4. Assessing Criterion A The IUCN Red List of Ecosystems criterion A requires estimates of the magnitude of change over time. This is typically calculated over a 50 year time period (past, present or future) or against a historical baseline ([IUCN 2024](https://portals.iucn.org/library/node/51533)). The first step towards achieving this change estimate over a fixed time frame is to assess the amount of change observed in your data. ### Area change Area change can be calculated for one of two scenarios: multiple ecosystems over two time slices using `getAreaChange()`, or one or more ecosystems over multiple time slices using `getAreaTrend()`. The `getAreaChange` function takes two inputs in the same format. It calculates the difference in area between two datasets, which can contain multiple ecosystems. It can accept data from `data.frame` such as those output by the `getArea()` function, `SpatRaster`, `SpatVector` or `sf` objects with POLYGON geometries. Make sure to provide area measure in $km^2$ if you are using a data.frame as the input. Note that when providing data frame and polygon inputs, the column name with unique ecosystem labels must be included in the function call. If the column name is different in the two inputs, both names must be provided in the same order as the input datasets. ```{r Using getAreaChange, eval = F} # input rasters area.change <- getAreaChange(mangrove.2000, mangrove.2017) # input data frames # value is the column that identifies the ecosystem in both tables. area.change <- getAreaChange(a.2000, a.2017, value) ``` The `getAreaTrend` can be used to work with multiple maps representing different ecosystem types and time frames. ### Rate of change In the Red List of Ecosystems, two methods are suggested to determine the rate of decline of an ecosystem, each of which assumes a different functional form of the decline ([IUCN 2024](https://portals.iucn.org/library/node/51533)). The proportional rate of decline (PRD) is a fraction of the previous year's remaining area, while the absolute rate of decline (ARD) is a constant fraction of the area of the ecosystem at the beginning of the decline ([IUCN 2024](https://portals.iucn.org/library/node/51533)). These rates of decline allow the use of two or more data points to extrapolate to the full 50 year timeframe required in an assessment. The annual rate of change (ARC) uses a compound interest law to determine the instantaneous rate of change ([Puyravaud 2003](https://www.sciencedirect.com/science/article/pii/S0378112702003353)). To calculate the rate of change using these methods and obtain forecasted estimates of change in area using extrapolation, we use the `declineForecast()` function. ```{r Using getDeclineStats} decline <- declineForecast(mangrove.2000, # ecosystem at time 1 mangrove.2017, # ecosystem at time 2 t1 = 2000, # year of time 1 year_diff = 17, # difference between time 1 and time 2 in years forecast_year = 2050) # year of desired forecast, defaults to t1+50 decline ``` Each method represent a different shape of decline. For further information about the choice of each of these methods to extrapolate refer to the IUCN Red List of Ecosystems guidelines ([IUCN 2024](https://portals.iucn.org/library/node/51533)). The results produced shows three estimated areas and estimated percent declines under the three estimation methods. It is important to note that this relatively simple exercise is founded on assumptions that should be fully understood before submitting your ecosystem assessment to the IUCN Red List of Ecosystems Committee for Scientific Standards. Furthermore, the guidelines suggest using area estimates from more than two time points to estimate change, and providing a measure of uncertainty if possible. Please see the guidelines ([IUCN 2024](https://portals.iucn.org/library/node/51533)) for more information. ## 4. Assessing Criterion B (distribution size) Criterion B utilizes measures of the geographic distribution of an ecosystem type to identify ecosystems that are at risk from catastrophic disturbances. This is done using two standardized metrics: the extent of occurrence (EOO) and the area of occupancy (AOO) [(Gaston and Fuller, (2009)](https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/j.1365-2664.2008.01596.x), [Keith et al., (2013)](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0062111)). It must be emphasised that EOO and AOO are not used to estimate the mapped area of an ecosystem like the methods we used in Criterion A; they are simply spatial metrics that allow us to standardise an estimate of risk due to spatially explicit catastrophes [(Murray et al., (2017)](https://onlinelibrary.wiley.com/doi/abs/10.1111/ddi.12533). Thus, it is critical that these measures are used consistently across all assessments, and the use of non-standard measures invalidates comparison against the thresholds. Please refer to the guidelines ([IUCN 2024](https://portals.iucn.org/library/node/51533)) for more information on AOO and EOO. ### Subcriterion B1 (calculating EOO) For subcriterion B1, we will need to calculate the extent of occurrence (EOO) of our data. We begin by creating the minimum convex polygon enclosing all occurrences of our focal ecosystem. ```{r Make EOO, eval = F} EOO.polygon <- getEOO(mangrove.2017) ``` There are some handy summary functions for the EOO object that is returned. We can also extract the area of the EOO polygon. ```{r View EOO, fig.width=7, fig.height=7, eval = F} #summarise and plot summary(EOO.polygon) plot(EOO.polygon) #area EOO.area <- getAreaEOO(EOO.polygon) EOO.area ``` This result can then be combined with additional information required under B1(a-c) to assess the ecosystem under subcriteria B1. ### Subcriterion B2 (calculating AOO) For subcriterion B2, we will need to calculate the number of 10x10 km grid cells occupied by our distribution. ```{r Creating AOO grid, eval = F} AOO.grid <- getAOO(mangrove.2017, cell_size = 10000) ``` Once more, the output object has some handy summary functions defined: ```{r view and summarise AOO, eval = F, fig.width=7, fig.height=7} # summary summary(AOO.grid) # handy plotting function plot(AOO.grid) ``` Both the `getEOO()` and `getAOO()` functions work on points and polygons too. #### Random grid search The `getAOO()` function uses a random grid search ("jitter") to find the optimal grid with the fewest cells. By default, the random grid search is omitted for ecosystems that are contained in more than 150 grid cells, because it is unlikely to yield results near the IUCN Criteria thresholds. However this option can be adjusted using the `jitter` argument (see below). We can run some diagnostics to ensure we have the best grid. The jitter-plot (jplot) shows how the minimum number of cells decreased as the random search iterations increased. The plot should fall to the minimum number early and stay flat. You can also plot a histogram of the AOO values generated with the random grid search iterations, and this should look like a bell curve in most cases, unless the number of grid cells is very low. ```{r Getting number of AOO grids, eval = FALSE, fig.alt='histogram of values generated with the jitter function'} # jplot jplot(AOO.grid) # histogram of values generated with each rep hist(AOO.grid) ``` If these plots look unconvincing, you can increase the number of iterations with the argument `n_jitter` directly in the `getAOO()` function. ```{r increase iterations, eval = F} # increase the number of iterations, if desired AOO.grid_n150 <- getAOO(mangrove.2017, cell_size = 10000, n_jitter = 150) # You can rerun diagnostics here to check for improvement jplot(AOO.grid_n150) hist(AOO.grid_n150) ``` You can also change the `jitter` argument to set the AOO threshold under which the random grid search is performed. For a fast approximation, turn off jitter (jitter = 0); to find a slower but more accurate solution, turn it on (jitter = 1) and to save time on larger datasets that are not near any threshold while maintaining accuracy for smaller datasets, set a custom threshold for the random grid search (defaults to jitter = 150). #### One percent rule As per the Red List of Ecosystems guidelines ([IUCN 2024](https://portals.iucn.org/library/node/51533)), the cells containing the smallest areas of the ecosystem, and collectively adding to less than 1% of the ecosystem's total area, could be excluded from the AOO calculation if their contribution to risk spreading is negligible. This prevents tiny, far-flung fragments of an ecosystem from inflating the AOO. This option is set per default but can be turned off using the argument, `bottom_1pct_rule = FALSE`, or the percentage changed using the `percent` argument. For assessments involving species, for instance, this option should be turned off. Here, we demonstrate the differences between including or not including the one percent rule and changing the percent parameter: ```{r One percent grid, eval = F} AOO.grid.keep.all <- getAOO(mangrove.2017, cell_size = 10000, bottom_1pct_rule = F, n_jitter = 150) AOO.grid.5.pct <- getAOO(mangrove.2017, cell_size = 10000, bottom_1pct_rule = T, n_jitter = 150, percent = 5) ``` Let's compare the results: ```{r View percent rule results, fig.width=14, fig.height=5, eval = TRUE, fig.alt='Testing options of AOO metric for 2017 mangrove distribution map', eval = F} par(mfrow = c(1,3)) plot(AOO.grid.keep.all, title = "Keep all") plot(AOO.grid_n150, title = "Remove bottom 1% (default)") plot(AOO.grid.5.pct, title = "Remove bottom 5%") par(mfrow = c(1,1)) ``` ## 5. Final words We have demonstrated in this vignette a typical workflow for assessing an example ecosystem under the RLE. The methods outlined here can easily be adapted for the Red List of Threatened Species as well, with slight adjustments in the parameters (specifically the cell_size parameter for `getAOO`). We hope that with this package, we help ensure a consistent implementation for both red lists criteria, minimising any misclassification as a result of using different software and methods.