Quickstart
Introduction
In the following we will demonstrate an idealized workflow based on a
subset of the Global Forest Watch (GFW) data set 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
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
article to learn about the most important concepts.
The idealized workflow for using {mapme.biodiversity}
consists of the following steps:
- prepare your sf-object containing only geometries of type
'POLYGON' or 'MULTIPOLYGON'
- decide which indicator(s) you wish to calculate and make the
required resource(s) available
- 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
spatial data format 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 part of the geometry of a protected area in the Dominican
Republic from the WDPA database.
library(mapme.biodiversity)
library(sf)
library(dplyr)
library(tidyr)
aoi_path <- system.file("extdata", "gfw_sample.gpkg", package = "mapme.biodiversity")
aoi <- st_read(aoi_path, quiet = TRUE)
aoi
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -71.73773 ymin: 18.63179 xmax: -71.69 ymax: 18.68691
#> Geodetic CRS: WGS 84
#> geom
#> 1 POLYGON ((-71.73417 18.6435...
Setting standard option
We use the mapme_options() function and set some
arguments, such as the output directory, that are important to govern
the subsequent processing. For this, we create a temporary directory.
Internally, to save time on downloading when building this vignette, we
copied already existing files to that output location (code not shown
here).
outdir <- file.path(tempdir(), "mapme-resources")
dir.create(outdir, showWarnings = FALSE)
mapme_options(
outdir = outdir,
verbose = TRUE
)
The outdir argument points towards a directory on the
local file system of your machine. All downloaded resources will be
written to respective directories nested within outdir.
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.
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
available_indicators() function:
available_indicators()
#> # A tibble: 25 × 3
#> name description resources
#> <chr> <chr> <list>
#> 1 biome Areal statistics of biomes from TEOW <tibble>
#> 2 burned_area Monthly burned area detected by MODIS satell… <tibble>
#> 3 deforestation_drivers Areal statistics of deforestation drivers <tibble>
#> 4 drought_indicator Relative wetness statistics based on NASA GR… <tibble>
#> 5 ecoregion Areal statstics of ecoregions based on TEOW <tibble>
#> 6 elevation Statistics of elevation based on NASA SRTM <tibble>
#> 7 fatalities Number of fatalities by group of conflict ba… <tibble>
#> 8 gsw_change Statistics of the surface water change layer… <tibble>
#> 9 gsw_occurrence Areal statistic of surface water based on oc… <tibble>
#> 10 gsw_recurrence Areal statistic of surface water based on re… <tibble>
#> # ℹ 15 more rows
available_indicators("treecover_area")
#> # A tibble: 1 × 3
#> name description resources
#> <chr> <chr> <list>
#> 1 treecover_area Area of forest cover by year <tibble [2 × 5]>
Say, we are interested in the treecover_area 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 treecover_area
documentation.
?treecover_area
help(treecover_area)
By inspecting the help page we learned that this indicator requires
the gfw_treecover and gfw_lossyear resources
and it requires to specify three extra arguments: the years for which to
calculate treecover, the minimum size of patches to be considered as
forest and the minimum canopy coverage of a single pixel to be
considered as forested.
With that information at hand, we can start to retrieve the required
resource. We can learn about all available resources using the
available_resources() function:
available_resources()
#> # A tibble: 23 × 5
#> name description licence source type
#> <chr> <chr> <chr> <chr> <chr>
#> 1 chirps Climate Hazards Group … CC - u… https… rast…
#> 2 esalandcover Copernicus Land Monito… CC-BY … https… rast…
#> 3 fritz_et_al Drivers of deforestati… CC-BY … https… rast…
#> 4 gfw_emissions Global Forest Watch - … CC-BY … https… rast…
#> 5 gfw_lossyear Global Forest Watch - … CC-BY … https… rast…
#> 6 gfw_treecover Global Forest Watch - … CC-BY … https… rast…
#> 7 global_surface_water_change Global Surface Water -… https:… https… rast…
#> 8 global_surface_water_occurrence Global Surface Water -… https:… https… rast…
#> 9 global_surface_water_recurrence Global Surface Water -… https:… https… rast…
#> 10 global_surface_water_seasonality Global Surface Water -… https:… https… rast…
#> # ℹ 13 more rows
available_resources("gfw_treecover")
#> # A tibble: 1 × 5
#> name description licence source type
#> <chr> <chr> <chr> <chr> <chr>
#> 1 gfw_treecover Global Forest Watch - Percentage of canopy… CC-BY … https… rast…
For the purpose of this vignette, we are going to download both, the
gfw_treecover and gfw_lossyear resources. 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 gfw_treecover
resource documentation.
?gfw_treecover
help(gfw_treecover)
?gfw_lossyear
help(gfw_lossyear)
We can now make the required resources available for our portfolio.
We will use a common interface that is used for all resources, called
get_resources(). We have to specify our portfolio object
and supply one or more resource functions with their respective
arguments. This will then download the matching resources to the output
directory specified earlier.
aoi <- get_resources(
x = aoi,
get_gfw_treecover(version = "GFC-2023-v1.11"),
get_gfw_lossyear(version = "GFC-2023-v1.11")
)
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 get_resources() function explained above. You will
have to re-run this function in every new R session, but note that data
that is already available will not be re-downloaded.
Here, we are going to calculate the treecover_area
indicator which is based on the resources from GFW. Since the resources
have 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 get_resources() beforehand.
aoi <- calc_indicators(
aoi,
calc_treecover_area(years = 2000:2023, min_size = 1, min_cover = 30)
)
#> Found a column named 'assetid'. Overwritting its values with a unique identifier.
Now let’s take a look at the results. In addition to the metadata we
are already familiar with, we see that there is an additional column
called treecover_area which contains a
tibble.
aoi
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -71.73773 ymin: 18.63179 xmax: -71.69 ymax: 18.68691
#> Geodetic CRS: WGS 84
#> # A tibble: 1 × 3
#> assetid treecover_area geom
#> <int> <list> <POLYGON [°]>
#> 1 1 <tibble [24 × 4]> ((-71.73735 18.64734, -71.71386 18.63179, -71.69 18…
The indicator is represented as a nested-list column in our
sf-object that is named alike the requested indicator. For
our single asset, this column contains a tibble with 6 rows and four
columns. Let’s have a closer look at this object
aoi$treecover_area
#> [[1]]
#> # A tibble: 24 × 4
#> datetime variable unit value
#> <dttm> <chr> <chr> <dbl>
#> 1 2000-01-01 00:00:00 treecover ha 1975.
#> 2 2001-01-01 00:00:00 treecover ha 1975.
#> 3 2002-01-01 00:00:00 treecover ha 1973.
#> 4 2003-01-01 00:00:00 treecover ha 1940.
#> 5 2004-01-01 00:00:00 treecover ha 1930.
#> 6 2005-01-01 00:00:00 treecover ha 1926.
#> 7 2006-01-01 00:00:00 treecover ha 1919.
#> 8 2007-01-01 00:00:00 treecover ha 1908.
#> 9 2008-01-01 00:00:00 treecover ha 1905.
#> 10 2009-01-01 00:00:00 treecover ha 1903.
#> # ℹ 14 more rows
The tibble follows a standard output format, which is the same for
all indicators. Each indicator is represented as a tibble with the four
columns datetime, variable, unit,
and value. In case of the treecover_area
indicator, the variable is called treecover and is
expressed in ha.
Let’s quickly visualize the results:
If you wish to change the layout of an portfolio, you can use
portfolio_long() and portfolio_wide() (see the
respective online
tutorial). 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.
geoms <- st_geometry(aoi)
portfolio_long(aoi, drop_geoms = TRUE)
#> # A tibble: 24 × 6
#> assetid indicator datetime variable unit value
#> <int> <chr> <dttm> <chr> <chr> <dbl>
#> 1 1 treecover_area 2000-01-01 00:00:00 treecover ha 1975.
#> 2 1 treecover_area 2001-01-01 00:00:00 treecover ha 1975.
#> 3 1 treecover_area 2002-01-01 00:00:00 treecover ha 1973.
#> 4 1 treecover_area 2003-01-01 00:00:00 treecover ha 1940.
#> 5 1 treecover_area 2004-01-01 00:00:00 treecover ha 1930.
#> 6 1 treecover_area 2005-01-01 00:00:00 treecover ha 1926.
#> 7 1 treecover_area 2006-01-01 00:00:00 treecover ha 1919.
#> 8 1 treecover_area 2007-01-01 00:00:00 treecover ha 1908.
#> 9 1 treecover_area 2008-01-01 00:00:00 treecover ha 1905.
#> 10 1 treecover_area 2009-01-01 00:00:00 treecover ha 1903.
#> # ℹ 14 more rows
Enabling parallel computing
{mapme.biodiversity} follows the parallel computing
paradigm of the {future}
package. That means that you as a user are in the control if and how you
would like to set up parallel processing. Since
{mapme.biodiversity} v0.7 supports parallel processing on a
nested-level. The outer level applies parallel processing to the assets
in a portfolio, the inner level to potential chunks for polygons that
are larger than the specified chunks size or the single components of
assets of type MULTIPOLYGON.
The maximum chunk size is specified in hectares via
mapme_options() and defaults to 100,000 ha. Polygons larger
than this threshold will be split into chunks of roughly the same
size.
Fine-control of parallel processing is given by using future
topologies (find more information here).
To process all assets sequentially, but allow to spawn up to 4 workers
to process chunks in parallel you might specify:
library(future)
plan(list(sequential, tweak(cluster, workers = 6)))
As another example, with the code below one would apply parallel
processing of 2 assets, with each having 4 workers available to process
chunks, thus requiring a total of 8 available cores on the host machine.
Be sure to not request more workers than available on your machine.
library(progressr)
plan(list(tweak(cluster, workers = 2), tweak(cluster, workers = 4)))
with_progress({
aoi <- calc_indicators(
aoi,
calc_treecover_area(
min_size = 1,
min_cover = 30
)
)
})
plan(sequential) # close child processes
Exporting an portfolio object
You can use the write_portfolio() function to save a
processed portfolio object to disk as a GeoPackage. This
allows sharing your data with contributors 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. You can use
read_portfolio() to read back a GeoPackage written in such
a way into R:
dsn <- tempfile(fileext = ".gpkg")
write_portfolio(x = aoi, dsn = dsn, quiet = TRUE)
from_disk <- read_portfolio(dsn, quiet = TRUE)
from_disk
#> Simple feature collection with 1 feature and 2 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -71.73773 ymin: 18.63179 xmax: -71.69 ymax: 18.68691
#> Geodetic CRS: WGS 84
#> # A tibble: 1 × 3
#> assetid treecover_area geom
#> <int> <list> <POLYGON [°]>
#> 1 1 <tibble [24 × 4]> ((-71.73735 18.64734, -71.71386 18.63179, -71.69 18…
#> [1] TRUE