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