---
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 ESA Landcover 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, we will also load the `{dplyr}` 
and `{tidyr}`packages. Then, we will read an 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(tidyr)

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 behavior 
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}
# copying package internal resource to a temporary location
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:2015,
  outdir = file.path(outdir, "res"),
  tmpdir = outdir,
  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 extent does 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 respective directories nested within `outdir`. 
Any intermediate files during the calculation of an indicator will be written 
to `tmpdir`. Thus, please ensure that you have write access to both directories 
and that there is sufficient free disk space to support the analysis of your 
portfolio. 

Once you request a specific resource for your portfolio, only those files will be 
downloaded that are missing to match its spatio-temporal extent. This behavior 
is beneficial, e.g. in case you share the `outdir` between different projects to 
ensure that only resources matching your current portfolio are returned.

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 
errors or warnings.

# 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 `population_count` indicator.
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 [population_count](https://mapme-initiative.github.io/mapme.biodiversity/reference/population_count.html) documentation.

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

By inspecting the help page we learned that this indicator requires the 
`worldpop` resource and it requires to specify two extra arguments: the population
statistic to calculate and the eninge to be used for the calculation (learn
more about engines [here](https://mapme-initiative.github.io/mapme.biodiversity/articles/terminology.html)).

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 `worldpop` 
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
[worldpop](https://mapme-initiative.github.io/mapme.biodiversity/reference/worldpop.html)
resource documentation.

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

We can now make the `worldpop` 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_esalandcover}
sample_portfolio <- get_resources(x = sample_portfolio, resources = "worldpop")
```

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("worldpop", "gfw_treecover"),
  vers_treecover = "GFC-2021-v1.9"
)
```

# 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 `population_count` indicator which is based on the `worldpop` 
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 = "population_count",
                                    stats_popcount = "sum", engine = "zonal")
```

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, population_count))
```

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 
6 rows and two columns. Let's have a closer look at one of these objects.

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

For each asset, the result is a tibble in long format indicating the population
sum per year (make sure to read the detailed indicator documentation via `?population_count`).
Let's quickly visualize the results for a single asset:

```{r plot_landcover, echo = FALSE}
data <- sample_portfolio %>%
  filter(assetid == 10) %>%
  st_drop_geometry() %>%
  unnest(population_count)

pop <- data$popcount_sum
names(pop) <- data$year
barplot(pop, main = "Population totals over time", 
        xlab = "Year", ylab = "Persons", 
        col = "steelblue")
```

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(population_count) %>%
  filter(assetid == 3)
```

# Enabling parallel computing

{mapme.biodiversity} follows the parallelization paradigm of the
{[future](https://cran.r-project.org/package=future)} 
package. That means that you as a user are in the control if and how you would 
like to set up parallel processing. Currently, `{mapme.biodiversity}` supports 
parallel processing on the asset level of the `calc_indicators()` function only. We
also currently assume that parallel processing is done on the cores of a single 
machine. In future developments, we would like to support distributed processing.
If you are working on a distributed use-cases, please contact the developers,
e.g. via the [discussion board](https://github.com/mapme-initiative/mapme.biodiversity/discussions) or mail.

To process e.g. 6 assets in parallel and report a progress bar you will have to set 
up the following in your code:

```{r parallel, eval = FALSE}
library(future)
library(progressr)

plan(multisession, workers = 6) # set up parallel plan with 6 concurrent threads

with_progress({
  portfolio <- calc_indicators(
    sample_portfolio,
    indicators = "population_count",
    stats_popcount = "sum", 
    engine = "zonal"
  )
})

plan(sequential) # close child processes
```

Note, that the above code uses `future::multisession()` as the parallel backend.
This backend will resolve the calculation in multiple background R sessions.
You should use that backend if you are operating on Windows, using R Studio or
otherwise are not sure about which backend to use. In case you are operating on 
a system that allows process forking and are *not* using R Studio, consider
using `future::multicore()` for more efficient parallel processing.

# 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)
```

