---
title: "Quickstart"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Quickstart}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# Introduction

In the following we will demonstrate an idealized workflow based on a subset
of the CHIRPS dataset that is delivered together with this package. You can follow
along the code snippets below to reproduce the results. Please note that
to reduce the time it takes to process this vignette, we will not download
any resources from the internet. In a real use case, thus processing time
might substantially increase because resources have to be downloaded and real
portfolios might be larger than the one created in this example. 

This vignette assumes that you have already followed the steps in [Installation](https://mapme-initiative.github.io/mapme.biodiversity/articles/installation.html) and have familiarized yourself with the 
terminology used in the package. If you are unfamiliar with the terminology used here, please head over to the [Terminology](https://mapme-initiative.github.io/mapme.biodiversity/articles/terminology.html)
article to learn about the most important concepts. 

The idealized workflow for using `{mapme.biodiversity}` consists of the following
steps:

- initialize your portfolio using an sf-object containing only geometries of type POLYGON
- decide which indicator(s) you wish to calculate and download the required resource(s)
- conduct your indicator calculation, which adds a nested list column to your portfolio object
- continue your analysis in R or decide to export your results to a GeoPackage to use it with
other geospatial software

# Getting started

First, we will load the `{mapme.biodiversity}` and the `{sf}` package for handling
spatial vector data. For tabular data handling and the creation of a plot, we will
also load the `{dplyr}` and `{ggplot2}` packages.
Then, we will read a package internal GeoPackage which includes the geometry of a 
protected area in the Dominican Republic from the WDPA database.

```{r setup, message=FALSE}
library(mapme.biodiversity)
library(sf)
library(dplyr)
library(ggplot2)

aoi_path <- system.file("extdata", "sierra_de_neiba_478140.gpkg", package = "mapme.biodiversity")
(aoi <- read_sf(aoi_path))
```

The sf-object contains a single object of geometry type `MULTIPOLYGON`. The 
`{mapme.biodiversity}` package, however, only supports geometries of type `POLYGON`,
thus we need to cast the geometry before we advance. The resulting sf object 
also contains some metadata, that will be retained throughout the complete
workflow. Because some of the cast geometries represent artefacts of the digitization 
process, in this example we will subset to only the largest polygon.

```{r cast}
(aoi <- st_cast(aoi, to = "POLYGON")[1, ])
```

In the following, we will simulate a portfolio consisting of several polygons (assets,
in the jargon of this package). To this end, we create smaller polygons within 
the original extent of the main polygon. This way, we can showcase the behaviour 
of the `{mapme.biodiversity}` package for portfolios that contain multiple assets. 
We will only select single assets with geometry type `POLYGON` that lie within 
the original boundary of the protected area.

```{r chunking}
aoi_gridded <- st_make_grid(
  x = st_bbox(aoi),
  n = c(10, 10),
  square = FALSE
) %>%
  st_intersection(aoi) %>%
  st_as_sf() %>%
  mutate(geom_type = st_geometry_type(x)) %>%
  filter(geom_type == "POLYGON") %>%
  select(-geom_type, geom = x) %>%
  st_as_sf()

metanames <- names(st_drop_geometry(aoi))
aoi_gridded[metanames] <- st_drop_geometry(aoi)
```

# Initialization of a portfolio

Now, we are ready to initiate a portfolio object containing multiple assets.
We use the <code>init_portfolio()</code> function and set some attributes that
are important for the subsequent processing. The function will add a unique
identifier column called 'assetid' that is used to uniquely identify each
asset in the portfolio.

```{r init_portfolio, dpi = 50, out.width="120%", out.height="120%"}
# copying package internal resource to tempdir
outdir <- file.path(tempdir(), "mapme.biodiversity")
dir.create(outdir)
resource_dir <- system.file("res", package = "mapme.biodiversity")
file.copy(resource_dir, outdir, recursive = TRUE)

sample_portfolio <- init_portfolio(
  x = aoi_gridded,
  years = 2010,
  outdir = file.path(outdir, "res"),
  tmpdir = outdir,
  add_resources = FALSE,
  cores = 1,
  verbose = TRUE
)
plot(sample_portfolio["assetid"])
```

The first argument, `x`, is the sf-object that we want to turn into a portfolio.
The argument `years` allows us to restrict our analysis to certain years only.
Certain resources with a temporal dimension are only processed for the portfolio's 
temporal extent. All resource and indicator functions will inform the 
user if the portfolio's temporal extents do not intersect. The `outdir` and `tmpdir` 
arguments point towards directories on the local file system of your machine. 
If these directories do not exist, the package attempts to create them.
The `outdir` cannot be equal to the `tmpdir` argument. All downloaded resources
will be written to their respective directories within `outdir`. 
We will set the `{terra}`'s package temporal directory to directory within `tmpdir`
and also any intermediate files during the calculation of an indicator will be written 
there. Thus, please ensure that you have write access to these directories and that
there is sufficient free disk space to support the analysis of your portfolio. 

In case your share a common `outdir` across different portfolios, with the 
`add_resources` arguments controls the behaviour of the portfolio initialization. 
If set to `TRUE`, the function will automatically search for all available 
resources and add these as attributes to the portfolio, without further checking
if these resources actually match its spatio-temporal extent. Only use this if 
you are sure that for your current portfolio, all resources have been previously
downloaded. If set to `FALSE`, no pre-existing resources will be added to the
portfolio. Once you request a specific portfolio for your resources, only
those files will be downloaded that are missing in `outdir`. This behaviour is
beneficial if you share the `outdir` between different projects to ensure that 
only matching resources are returned.

The argument `cores` specifies the number of CPU cores available to the indicator 
routines for parallel processing. Parallel processing is deactivated on Windows 
by default. If you use another operating system, you may set to a value greater 
1 if you wish to activate parallel processing. 

Finally, the `verbose` logical controls whether or not the package
will print informative messages during the calculations. Note, 
that even if set to `FALSE`, the package will inform users about any potential 
error in specifying of arguments in the form of warnings and errors.

# Getting the right resources

You can check which indicators are available via the <code>available_indicators()</code>
function. For this, we inspect the names of the returned object:

```{r query_indicator}
names(available_indicators())
```

Say, we are interested in the precipitation indicator based on the CHIRPS dataset.
We can learn more about this indicator and its required resources by using
either of the commands below or, if you are viewing the online version, head
over to the [precipitation_chirps](https://mapme-initiative.github.io/mapme.biodiversity/reference/precipitation_chirps.html) documentation.

```{r helppage_indicator, eval = FALSE}
?precipitation_chirps
help(precipitation_chirps)
```

By inspecting the help page we learned that this indicator requires the `chirps` 
resource and it accepts three arguments, namely the `scales_spi`, `spi_prev_years`
and `engine` arguments. `scales_spi` determines the time-scale for the calculation of the Standardized-Precipitation-Index (SPI), while `spi_prev_years` indicates how many 
previous years are included in the fitting process.
The `engine` argument controls how zonal statistics for each asset in the portfolio
are extracted.

With that information at hand, we can start to retrieve the required resource. 
We can learn about all available resources using the <code>available_resources()</code> 
function:

```{r query_resources}
names(available_resources())
```

For the purpose of this vignette, we are going to download the `chirps` resource.
We can get more detailed information about a given resource, by using either of 
the commands below to open up the help page. If you are viewing the online 
version of this documentation, you can simply head over to the [chirps](https://mapme-initiative.github.io/mapme.biodiversity/reference/chirps.html)
resource documentation.

```{r helppage_resource, eval = FALSE}
?chirps
help(chirps)
```

We can now make the `chirps` resource available for our portfolio. We will
use a common interface that is used for all resources, called <code>get_resources()</code>.
We have to specify our portfolio object and the names of the resource(s) we wish
to download. Additional arguments for the specific resource can be specified.
The output of the function is the portfolio object with its attributes appended
for the new resource, thus we simply can overwrite the `sample_portfolio` variable.

```{r get_chirps}
sample_portfolio <- get_resources(x = sample_portfolio, resources = "chirps")
```

In case you want to download more than one resource, you can use the same interface
and the resources will be made available sequentially. Required arguments 
for a resource are simply added as usual:

```{r get_multi_resources, eval = FALSE}
sample_portfolio <- get_resources(
  x = sample_portfolio,
  resources = c("chirps", "treecover2000"),
  vers_treecover = "GFC-2020-v1.8"
)
```

# Calculate specific indicators

The next step consists of calculating specific indicators. Note 
that each indicator requires one or more resources that were made available 
via the <code>get_resources()</code> function explained above. Here, we are going
to calculate the `precipitation_chirps` indicator which is based on the `chirps` resource.
Since the resource has been made available in the previous step, we can continue
requesting the calculation of our desired indicator. Note the command below
would issue an error in case a required resource has not been made
available via <code>get_resources()</code> beforehand.

```{r calc_indicator}
sample_portfolio <- calc_indicators(sample_portfolio,
  indicators = "precipitation_chirps",
  scales_spi = 3,
  spi_prev_years = 8,
  engine = "extract"
)
```

The function call will inform us that we have not specified the two required 
arguments but that the default values have been inserted. If you get such a 
message for other indicators and you do not know what they mean, head over to 
the indicator help-page to learn about the available arguments.

Now let's take a look at the results. We will select only some of the metadata
and the output indicator column to get a clearer picture of what has happened.

```{r select_cols}
(sample_portfolio <- sample_portfolio %>% select(assetid, WDPAID, precipitation_chirps))
```

We obtained a new listed column in our sf object that is called like the requested
indicator. For each asset in our portfolio, this column
contains a tibble with 12 rows and 4 columns. Let's have a closer look at these 
objects.

```{r investigate_indicator}
sample_portfolio$precipitation_chirps[10]
```

For each asset, the result is a tibble in long format for 12 months in the year 
2010. The precipitation indicator with the default setting delivered the absolute
rainfall sum as well as the rainfall anomaly compared to the 1981-2010 climate 
normal period (make sure to read the detailed indicator documentation via `?precipitation_chirps`).
Let's quickly visualize the results for a single asset:

```{r plot_precipitation, echo = FALSE, warning=FALSE, dpi = 50}
sample_portfolio %>%
  filter(assetid == 10) %>%
  st_drop_geometry() %>%
  tidyr::unnest(precipitation_chirps) %>%
  mutate(sign = ifelse(anomaly < 0, "lower than average", "higher than average")) %>%
  ggplot() +
  geom_bar(aes(x = dates, y = anomaly, fill = sign), stat = "identity") +
  scale_fill_manual(values = c("darkblue", "brown2")) +
  labs(
    title = "Monthly precipitation anomaly for 2010 in regard to the 1981-2010 climate normal", x = "", y = "Precipitation anomaly [mm]",
    fill = "Anomaly"
  ) +
  theme_classic() +
  theme(legend.position = "bottom")

sample_portfolio %>%
  filter(assetid == 10) %>%
  st_drop_geometry() %>%
  tidyr::unnest(precipitation_chirps) %>%
  mutate(
    spi_3 = ifelse(is.na(spi_3), 0, spi_3),
    sign = ifelse(spi_3 < 0, "lower than average", "higher than average")
  ) %>%
  ggplot() +
  geom_bar(aes(x = dates, y = spi_3, fill = sign), stat = "identity") +
  scale_fill_manual(values = c("darkblue", "brown2")) +
  labs(
    title = "Monthly SPI value for 2010 (timescale = 3)", x = "Month", y = "SPI",
    fill = "SPI"
  ) +
  theme_classic() +
  theme(legend.position = "bottom")

sample_portfolio %>%
  filter(assetid == 10) %>%
  st_drop_geometry() %>%
  tidyr::unnest(precipitation_chirps) %>%
  ggplot() +
  geom_bar(aes(x = dates, y = absolute), stat = "identity", fill = "darkblue") +
  labs(title = "Sum of monthly precipitation 2010", x = "", y = "Precipitation sum [mm]") +
  theme_classic()
```

If you wish to conduct your statistical analysis in R, you can use `{tidyr}` functionality
to unnest one or multiple columns. Especially for large portfolios, it is usually
a good idea to keep the geometry information in a separated variable
to keep the size of the data object relatively small.

```{r unnest}
geometries <- select(sample_portfolio, assetid)
sample_portfolio %>%
  st_drop_geometry() %>%
  tidyr::unnest(precipitation_chirps) %>%
  filter(assetid == 3)
```

# Exporting an portfolio object

You can use the <code>write_portfolio()</code> function to save a processed 
portfolio object to disk as a GeoPackage. This allows sharing your data with 
others who might not be using R, but any other geospatial software. Simply 
point towards a non-existing file on your local disk to write the portfolio. 
The function will create an individual table for all processed indicators. Via 
the <code>read_portfolio()</code> function, a portfolio which has been written 
to disk in such a way can be read back into R. However, users should note that 
the portfolio-wide arguments that were set during the portfolio initialization 
are **not**  reconstructed (e.g. the temporal extent, `outdir` and `tmpdir`, etc.).
Thus if you wish to continue to use `{mapme.biodiversity}` functionality on 
such a portfolio object, make sure to re-run <code>init_portfolio()</code> 
on it.

```{r portfolio_io}
tmp_output <- tempfile(fileext = ".gpkg")
write_portfolio(
  x = sample_portfolio,
  dsn = tmp_output
)
(portfolio_from_disk <- read_portfolio(tmp_output))
```

```{r delete_tmp, echo=FALSE, include=FALSE}
file.remove(tmp_output)
```

