---
title: "Calculating spatial metrics for IUCN red list assessments"
author: "Calvin K.F. Lee, Nicholas J. Murray"
date: "2017-07-07"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Calculating spatial metrics for IUCN red list assessments}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
## 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 the assessment. Please ensure
that these are already installed.
```{r Loading packages, message=FALSE}
library(sf)
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/10.1111/j.1466-8238.2010.00584.x/abstract)),
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.

```{r Loading our example distributions}
mangrove.2000 <- raster(system.file("extdata", "example_distribution_2000.tif", 
                                    package = "redlistr"))
mangrove.2017 <- raster(system.file("extdata", "example_distribution_2017.tif", 
                                    package = "redlistr"))
```

### Importing data as rasters
To import your own raster data, you can use the `raster::raster` function, which
can handle many different types of georeferenced spatial data. More information
on the `raster` package can be found 
[here](https://CRAN.R-project.org/package=raster/raster.pdf).

### Importing data as shapefile or .KML file
To import data as a shapefile (.shp), we can use the `sf::st_read` function.
```{r Importing shapefile, eval=FALSE}
library(sf)
my.shapefile <- st_read('./path/to/folder/', 'shapefile.shp')
my.KML.file <- st_read('./path/to/folder/kmlfile.kml')
```

Once these files are imported, they can then be converted into raster format via
the `raster::raster` or `raster::rasterFromXYZ` functions.

### Plotting out data
We plot the data to view the two 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}
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.

At this stage it is important to check that the data are georectified (in the
right location on Earth). This can be acheived by simply plotting the map and
checking the coordinates. Otherwise,  check the distribution maps against
satellite images, perhaps by using packages such as `ggmap`, `plotgooglemaps`,
`googleVis` etc. It is also important to check that the data are suitable for
the task at hand.

## 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 it is necessary to know
the data by investigating their resolution (grain size), extent and spatial
properties. Specifically, the coordinate reference system (CRS) of the datasets
must be measured in metres, and is consistent to compare area.

```{r Checking CRS}
isLonLat(mangrove.2000) 
isLonLat(mangrove.2000)
# If TRUE, they must be reprojected to a projected coordinate system

crs(mangrove.2000)@projargs == crs(mangrove.2017)@projargs
```

Reprojection can be done by using the `raster::projectRaster` function.

We start by calculating the area of the dataset. 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}
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` km2, and at time
2, (2017) the ecosystem was `r a.2017` km2.

### Creating binary ecosystem raster
An additional parameter in `getArea` is to specify which class to count if your 
raster has more than one value (multiple ecosystem data stored within a single
file). However, for further functions, it may be wise to convert your raster
into a binary format, containing only information for your target ecosystem.

First, if you do not know the value which represents your target ecosystem, you 
can plot the ecosystem out, and use the `raster::click` function. Once you know
the value which represents your target ecosystem, you can create a new raster
object including only that value. Here's an example with a dummy raster.
```{r Binary object from multiclass, eval=FALSE}
# Create dummy raster for example
r <- raster(nrows=10, ncols=10)
r.multiple <- r
values(r.multiple) <- rep(c(1:10), 10)

# If the target ecosystem is represented by value = 5
r.bin <- r.multiple == 5 # Has values of 1 and 0
# Convert 0s to NAs
values(r.bin)[values(r.bin) != 1] <- NA
```

## 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 ([Bland et al., 2016](https://portals.iucn.org/library/sites/library/files/documents/2016-010.pdf)).
The first step towards acheiving 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 using the `getAreaLoss` function. It simply
calculates the difference in area between two datasets. It can accept data from
numbers, `RasterLayers` or `SpatialPolygons`. Make sure to use area measure in
km2 if you are using numbers as the input.

```{r Using getAreaChange}
area.lost <- getAreaLoss(a.2000, a.2017)
# getAreaLoss(mangrove.2000, mangrove.2017) generates identical results

```

### 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 ([Bland et al., 2016](https://portals.iucn.org/library/sites/library/files/documents/2016-010.pdf)).
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 ([Bland et al., 2016](https://portals.iucn.org/library/sites/library/files/documents/2016-010.pdf)).
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 estimate the rate of changing using these methods, we use the
`getDeclineStats` function. 

```{r Using getDeclineStats}
decline.stats <- getDeclineStats(a.2000, a.2017, 2000, 2017, 
                                 methods = c('ARD', 'PRD', 'ARC'))
decline.stats
```

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 ([Bland et al., 2016](https://portals.iucn.org/library/sites/library/files/documents/2016-010.pdf)).

Now, it is possible to extrapolate, using only two estimates of an ecosystems'
area, to the full 50 year period required for a Red List of Ecosystems
assessment. 

```{r Estimating future area}
extrapolated.area <- futureAreaEstimate(a.2000, year.t1 = 2000, 
                                        ARD = decline.stats$ARD, 
                                        PRD = decline.stats$PRD, 
                                        ARC = decline.stats$ARC, 
                                        nYears = 50)
extrapolated.area
```

50 years from our first estimate is the year
`r extrapolated.area$forecast.year`.

As we included all three methods of calculating rate of decline currently
included in the package, the results produced shows three estimated areas under
the various decline scenarios. 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 ([Bland et al., 2016](https://portals.iucn.org/library/sites/library/files/documents/2016-010.pdf))
for more information.

If we were to use the Proportional Rate of Decline (PRD) for our example
assessment, our results here will be suitiable for criterion A2b (Any 50 year
period), and the percent loss of area is:
```{r Percent loss}
predicted.percent.loss <- (extrapolated.area$A.PRD.t3 - a.2000)/a.2000 * 100
predicted.percent.loss
```


## 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://onlinelibrary.wiley.com/doi/10.1111/j.1365-2664.2008.01596.x/full), 
[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/10.1111/ddi.12533/abstract). 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 ([Bland et al., 2016](https://portals.iucn.org/library/sites/library/files/documents/2016-010.pdf))
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, fig.width=7, fig.height=7}
EOO.polygon <- makeEOO(mangrove.2017)
plot(EOO.polygon)
plot(mangrove.2017, add = T, col = "green", legend = FALSE)
```

We can then calculate the area of the EOO polygon.

```{r Calculating EOO area} 
EOO.area <- getAreaEOO(EOO.polygon)
EOO.area
```

The calculated EOO area for this subset of mangroves in Victoria is `r EOO.area`
km2. 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. We begin by creating the appopriate grid cells.

```{r Creating AOO grid, fig.width=7, fig.height=7}
AOO.grid <- makeAOOGrid(mangrove.2017, grid.size = 10000,
                            min.percent.rule = F)
plot(AOO.grid)
plot(mangrove.2017, add = T, col = "green", legend = FALSE)
```

Finally, we can use the created grid to calculate the AOO for our mangroves.

```{r Getting number of AOO grids}
n.AOO <- length(AOO.grid)

# the getAOO function can also be used to directly get the AOO
# n.AOO <- getAOO(mangrove.2017, grid.size = 10000, 
#                 min.percent.rule = T, percent = 0.1)
n.AOO
```


#### Grid uncertainty functions
Although we have a number for AOO `r n.AOO`, this may not be the actual minimum
AOO for this ecosystem, because the placement of the AOO grid will influence this
number. To take this into account, we have included a few additional functions.
We will be using `gridUncertainty` in this vignette.

`gridUncertainty` simply moves the AOO grid systematically (with a
small random movement around fixed points), and continues searching  for a minimum 
AOO until additional shifts no longer produce improved results. 
```{r gridUncertainty}
gU.results <- gridUncertainty(mangrove.2017, 10000,
                              n.AOO.improvement = 5, 
                              min.percent.rule = F)
# If it takes your computer very long to run this command, consider reducing 
# n.AOO.improvement

gU.results$min.AOO.df
```
The result here shows that the minimum AOO does indeed change in response to 
where the AOO grid is placed, and that our original AOO value of `r n.AOO` is
not the true minimum AOO. The true minimum AOO (and the one which should be used
for the assessment) is `r gU.results$min.AOO.grid$AOO.number`.

It is also possible to plot out the grid which was used to produce the minimum 
AOO result, which can be used to report our final results.

```{r Plotting out minimum AOO grid, fig.width=7, fig.height=7}
plot(gU.results$min.AOO.grid$out.grid)
plot(mangrove.2017, add = T, col = "green", legend = FALSE)
```

#### Exporting results
The results of the above plot from `gridUncertainty` is a shapefile, which can
be exported with the `raster::shapefile` function.

#### One percent rule
In addition to the size of the grids used (which will be different for species 
assessments), there is also an option in the `makeAOOGrid` function to specify 
whether a minimum of percent of the grid cell area must be occupied before they 
are included as an AOO grid. Typically, this is only used in very special cases
for ecosystems with highly skewed distribution of patch sizes ([Bland et al., 2016](https://portals.iucn.org/library/sites/library/files/documents/2016-010.pdf)).

Here, we demonstrate the differences between including, or not including the one
percent rule:

```{r One percent grid, fig.width=7, fig.height=7}
AOO.grid.one.percent <- makeAOOGrid(mangrove.2017, grid.size = 10000, 
                                    min.percent.rule = T, percent = 1)
plot(AOO.grid.one.percent)
plot(mangrove.2017, add = T, col = "green", legend = FALSE)
```

There is an additional parameter - `percent` - which adjusts the threshold for
the AOO grids. Here, we set it to 0.1% to demonstrate its functionalities.

```{r AOO Grid 0.1percent, fig.width=7, fig.height=7}
AOO.grid.min.percent <- makeAOOGrid(mangrove.2017, grid.size = 10000,
                                    min.percent.rule = T, percent = 0.1)
par(mfrow = c(2,2))
plot(AOO.grid, main = 'AOO grid without one percent rule')
plot(mangrove.2017, add = T, col = "green", legend = FALSE)
plot(AOO.grid.one.percent, main = 'AOO grid with one percent rule')
plot(mangrove.2017, add = T, col = "green", legend = FALSE)
plot(AOO.grid.min.percent, main = 'AOO grid with one percent rule at 0.1%')
plot(mangrove.2017, add = T, col = "green", legend = FALSE)
```

## 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 grid.size parameter for `getAOO` or
`gridUncertainty`). 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.

