In EhrmannS/rasterTools: Obtain and process earth observation data

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Introduction

You might find rasterTools useful when your work depends in one way or the other on earth observation data [@bush_connecting_2017]. This may be for simulating spatial patterns, building species distribution models or analysing and measuring spatial patterns/landscape features.

rasterTools does NOT intend to replace raster [@hijmans_raster_2017] but tries to present an additional workflow that is aimed at being more transparent and reproducible and easily accessible. Many of the typical steps in working with raster data are summarised into operators (functions that carry out a specific spatial operation) and are grouped according to the group of tasks. In rasterTools the four functions obtain(), generate(), modify() and measure() represent the major tasks of spatial operations. rasterTools has been conceptualised so that these core functions are modular. Each of these functions manages all the code-logic that is common for the task at hand, e.g. to obtain datasets. More specific code that manages, for instance, "loading the MODIS dataset" is oursourced to a seperate function (the oMODIS() operator). Hence, new operators to all the core functions can be devised easily.

Several operators are combined to algorithms to carry out (a sequence of) more complex spatial operations. Operators are scope specific, i.e. they are typically only called in their respective core function. An algorithm in rasterTools is given as list of lists.

algo_name <- list(list(operator1),
                  list(operator2))

Each operator represents a do.call() compliant definition (type ?do.call into the console to learn about it).

do.call(what, args, quote = FALSE, envir = parent.frame())

Which translates to

algo_name <- list(operator1 = list(what, args),
                  operator2 = list(what, args))

Its first element what = operator_name would outline the name of the function utilized in this operator (such as "oMODIS" or "oCLC"). Its second element (args) can have any number of arguments required to put together the call (such as product = "MOD17A3", layer = 2, period = c(2012, 2014) or simply period = 2006).

datasets <- list(list(operator = "oCLC", period = 2006),
                 list(operator = "oGFC", period = c(2005:2007)))

This modularity is one of the strong sides of rasterTools. It allows the user to write their own operator to be used in the respective core-function. An example would be a dataset that is not yet supported by rasterTools, to utilize a modification which is not yet included in this package or to develop new landscape metrics. All of this can easily and seamlessly be combined with the previously defined operators in the same, simple workflow.

Example workflow

In the following you will find an example for each of the core functions.

rasterTools may write large files to your harddisc and you should assign a specific directory on a storage medium that has sufficient space via setPaths(root = "/path/to/the/spatial/files").

Obtain spatial datasets

An algorithm with which one would obtain information from various spatial data sets could be:

library(rasterTools)

myDatasets <- list(list(operator = "oGFC", period = c(2006)),
                   list(operator = "oMODIS", product = "mod17a3", period = 2006,
                        layer = 2))

Most of the different spatial datasets are available in a specific file format and may have other specific properties, which require consideration. For instance, MODIS data are stored as products with several layers within each file. Consequently the product and layer need to be specified. Each of the required arguments can be found in the respective documentation, for instance with ?oMODIS.

Typically we only need a (vastly) smaller subset of an earth observation dataset and we would want to use a rectangular mask to outline this area of interest:

myMask <- rtGeoms$mask

The collection of datasets for the area of our interest can be obtained by the following code:

myData <- obtain(data = myDatasets, mask = myMask)

(output here not shown as it may download the datasets, if you do not have them yet)

Modify raster objects

A typical modification algorithm determines patches of the foreground in a raster with continuous integer values:

get_patches <- list(list(operator = "rBinarise", thresh = 30),
                    list(operator = "rPatches"))

When using this algorithm in modify() while sequential = TRUE the binarised raster will be used to derive patches. The output of rBinarise will be passed to rPatches, which determines patches from the thresholded rasters values.

However, if you want to determine patches, but also derive categories from the original raster, you have to give each operator in the algorithm a specific name. Operators which should be part of the same sub-algorithm have to have the same name:

cc_cats <- list(get_patches = list(operator = "rBinarise", thresh = 30),
                get_patches = list(operator = "rPatches"),
                get_categories = list(operator = "rCategorise", n = 5))

In this case sequential = TRUE would be ignored globally, but applied locally for those functions that share the same name (also if it would be set to FALSE in the call). You can find many more examples in the Usecases vignette.

The modifications would be carried out by:

myInput <- rtRasters$continuous
myPatches <- modify(input = myInput, by = get_patches, sequential = TRUE)
# visualise(raster::stack(myInput, myPatches))

Measure raster objects

Landscape metrics are another case of algorithms in rasterTools. We distinguish between generic and derived landscape metrics. To evaluate a generic metric, typically just the respective operator with its arguments needs to be given:

myInput <- rtRasters$categorical
myMetrics <- list(a_c = list(operator = "mArea", scale = "class"),
                  a_l = list(operator = "mArea", scale = "landscape"))
measure(input = myInput, with = myMetrics)

The so defined terms are the basis for derived metrics. Such a metric is defined simply by its mathematical equation, for example the Class proportional area is the class areas divided by the landscape area:

myMetric <- list(a_c = list(operator = "mArea", scale = "class"),
                 a_l = list(operator = "mArea", scale = "landscape"),
                 mCPA = "a_c / a_l * 100")

Several derived metrics can be defined in one algorithm, given the required terms are defined:

myMetrics <- list(a_p = list(operator = "mArea", scale = "patch"),
                  a_c = list(operator = "mArea", scale = "class"),
                  a_l = list(operator = "mArea", scale = "landscape"),
                  mCPA = "a_c / a_l * 100",
                  mLPI = "max(a_p) / a_l * 100")

The generic metrics mNumber(), mPerimeter(), mArea(), mAdjacency() and mValues() are available. Equations for the derived metrics can be found in the landscape metrics vignette (for example the largest patch index mLPI or the percentage of like adjacencies mPLA; and with increasing version number more of them). It becomes apparent that also new, possibly improved landscapes metrics can be easily prototyped and tested in this framework.

Finally, the measurement of the class proportional area and the largest patch index, in a raster where patches have been determined, would be carried out by:

measure(input = myInput, with = myMetrics)

Acknowledgements

I am grateful for financial support from the PROFOUND Cost-action, which gave me the opportunity to work in a concentrated effort a large part of the functionality. This package has been developed in support of the FunBo Project, which was made possible by the Grünewald-Zuberbier Scholarship handed out by the University of Freiburg.

Thanks are also due to Prof. Arne Pommerening who was a great source of inspiration for what rasterTools is now.