Calculating spatial metrics for IUCN RLE assessments

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. 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.

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), 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.

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.

plot(mangrove.2000, col = "grey30", legend = FALSE, main = "Mangrove Distribution")
plot(mangrove.2017, add = T, col = "darkorange", legend = FALSE)

Overlay of two mangrove distribution maps

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.

is.lonlat(mangrove.2000)

## [1] FALSE

is.lonlat(mangrove.2017)

## [1] FALSE

# If TRUE, they must be transformed to a projected coordinate system

crs(mangrove.2000) == crs(mangrove.2017)

## [1] TRUE

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.

a.2000 <- getArea(mangrove.2000)
a.2000
a.2017 <- getArea(mangrove.2017)
a.2017

The results: at time 1 (2000), the ecosystem was 13.6590844 km2, and at time 2, (2017) the ecosystem was 12.3366736 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:

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). 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.

# 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. See vignette Scaling up Red List of Ecosystem assessments for instruction on how to work with multiple ecosystem maps.

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). 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). 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).

To calculate the rate of change using these methods and obtain forecasted estimates of change in area using extrapolation, we use the declineForecast() function.

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

## [[1]]
##   value  area.x  area.y area_diff percent_diff
## 1     1 13.6557 12.3336   -1.3221    -9.681671
## 
## [[2]]
##   value method        rate forecast.year forecast.area area.change pct.change
## 1     1    ARC -0.59899866          2050     10.121458   -3.534242  -25.88108
## 2     1    ARD  0.07777059          2050      9.767171   -3.888529  -28.47550
## 3     1    PRD  0.59720824          2050     10.121458   -3.534242  -25.88108

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).

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) 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), Keith et al., (2013)). 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). 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) 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.

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.

#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.

AOO.grid <- getAOO(mangrove.2017, cell_size = 10000)

Once more, the output object has some handy summary functions defined:

# 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.

# 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.

# 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 grid is jittered. 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), 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:

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:

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.