---
title: "Introduction to the 'chronosphere' R package"
author: "Adam T. Kocsis, Nussaïbah B. Raja"
date: '`r Sys.Date()`'
output: pdf_document
vignette: >
  %\VignetteIndexEntry{Introduction to the 'chronosphere' R package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(chronosphere)
```

# 1. Introduction

## 1.1. Installation

To install this alpha version of the package, you must download it either from the CRAN servers or its dedicated GitHub repository (http://www.github.com/adamkocsis/chronosphere/). All minor updates will be posted on GitHub as soon as they are finished, so please check this regularly. The version on CRAN will be lagging for some time, as it takes the servers many days to process everything. All questions and bugs can be reported at the GitHub issues board (https://github.com/adamkocsis/chronosphere/issues). Instead of spending it on actual research, a tremendous amount of time was invested in making this piece of software streamlined and user-friendly. If you use a dataset of the package in a publication, please cite both its reference(s) and the `chronosphere` package itself.

After installing from CRAN, from a source, or with ``devtools::install_github()`` , you can attach the package with:

```{r package, echo= TRUE, eval=FALSE}
library(chronosphere)
```

# 1.2 General features

The purpose of the `chronosphere` project is to facilitate, streamline and speed up the finding, acquisition and importing of Earth science data in R. Although the package currently focuses on deep time datasets, the scope of the included datasets will be much larger, spanning from a single variable published as supplementary material in a journal article to GIS data or the entire output of GCM models. `chronosphere` intends to decrease the gap between research hypotheses and the finding, downloading and importing of datasets. 

# 2. RasterArray

Faster data importing and better organization represents a considerable part of this process. Spatially explicit data are excellent candidates to demonstrate how more efficient data organization can speed up research. Although R has an excellent infrastructure for handling raster data (Hijmans & van Etten, 2019), the arrangement of individual layers are rather limited, which can be a problem when a large number of layers are considered. RasterStacks and RasterBricks are very efficient for organizing RasterLayers according to a single dimension (e.g. depth for 3D variables, or time), but multidimensional structures are preferred. 

To offer a more effective solution, the `chronosphere` package includes the definition of the RasterArray S4 class. RasterArrays represent hybrids between RasterStacks and regular R arrays. In short, they are virtual arrays of RasterLayers, and can be thought of as regular arrays that include entire rasters as elements rather than single numeric, logical or character values. As regular R users are familiar with subsetting, combinations and structures of regular arrays (including formal vectors and matrices), the finding, extraction and storage of spatially explicit data is much easier in such containers.

## 2.1. Structure

RasterArrays do not directly inherit from Raster* objects of the raster package, as a considerable number of main functions differ, but they rather represent a wrapper object around regular RasterStacks - stacks of individual RasterLayers. This ensures that whenever users are unfamiliar with the methods of RasterArray class objects, they can always reduce their data to stacks or individual layers. 

The formal class RasterArray has only two slots: a stack and an index. The stack includes the Raster data in an unfolded manner, similarly to how matrices and arrays are essentially just vectors with additional attributes. The stack slot incorporates a single RasterStack object, which represents the data content of the object. The index slot, on the other hand, describes the structure of the RasterArray. It is a vector/matrix/array of integers each representing an index of the layers in the stack. The configuration (dimensions) of the index represents the entire array.

The `chronosphere` package includes two demo data sets: a set of ancient topographies (Scotese and Wright, 2018) and a time series of bio-climatic variables (annual mean temperature and precipitation) from the CHELSA project (Karger et al, 2017a, 2017b). These datasets can be attached with the data command. 
```{r dat, echo= TRUE}
data(dems)
data(clim)
```

The structure of RasterArrays can be inspected if the object's name is typed into the console:
```{r dems, echo= TRUE}
dems
```

The first part of the console output includes properties of the individual RasterLayers stored in the stack. These layers have to share essential attributes that allow them to be stored in a single stack (extent, resolution, CRS). 

The second part of the output is a visualization of the structure of the RasterArray itself. In the case of the DEMs, 10 layers are stored in the stack, each layer having its individual name (e.g. `dem_0`). It is a single dimensional array (vector), and each element has its name in the array (`0`). Differentiating between the names of layers and the names of elements allows different subsetting and replacement rules for the two, which can both be handy - depending on the needs of the user.

The structure of the RasterArray can be visualized, analyzed or processed using the the proxy object. Proxies are essentially the same as the index slots of the RasterArray, but instead of including the indices of the layers they represent, proxies include the names of the layers. These can be accessed using the `proxy()` function.

```{r dems2, echo= TRUE}
proxy(dems)
```

Proxies are displayed as the second parts of the console output when the name of the object is typed into the console (show method). 

```{r clim, echo= TRUE}
clim
```

This RasterArray has 10 rows (annual means) and two variables/columns: temperature (`bio1`) and precipitation (`bio12`). With the `proxy()` function it is easy to interact with this object, or to query or analyze it.

```{r climprox, echo= TRUE}
proxy(clim)
```

RasterArrays are fairly easy to construct: one only needs a stack of the data and a regular vector/matrix/array including integers. For instance, the dems object can be recreated from scratch without any problem.

```{r constructor, echo= TRUE}
# a stack of rasters
stackOfLayers <- dems@stack
# an index object
ind <- 1:10
names(ind) <- letters[1:10]
# a RasterArray
nra  <- RasterArray(index=ind, stack=stackOfLayers)
nra
```

The attributes of the index object are defining the structure of the RasterArray. RasterArrays can be created with the combination of individual RasterLayers (or RasterArrays) using the combine() function.

```{r comb, echo= TRUE}
# one raster
r1 <- raster()
values(r1) <- 1
# same structure, different value
r2 <-raster()
values(r2) <- 2
comb <- combine(r1, r2)
comb
```

Matrix-like RasterArrays can also be created easily with the cbind() and rbind() functions.


```{r rbind, echo= TRUE}
# bind dems to itself
cbind(dems, dems)
```


## 2.2. RasterArray attributes and function to query

Functions that query and change attributes of the RasterArray resemble general arrays more than Raster* objects. They are connected to the index slot of the RasterArray and return values accordingly. 

The number of elements represented in the RasterArray can be queried with the `length()` function:
```{r le, echo= TRUE}
length(dems)
```

This RasterArray has 10 elements. The number of column and row names can be queried in a similar way:

```{r cle, echo= TRUE}
nrow(clim)
ncol(clim)
```

These functions are summarized in the dim() function. This, however, unlike the regular `dim()` method of vectors, return the length of the RasterArray vector, rather than just `NULL`. 

```{r dime, echo= TRUE}
dim(dems)

dim(clim)
```

The organization of RasterLayers can be greatly facilitated with names. The `names()`, `colnames()`, `rownames()` and `dimnames()` functions work the same way on RasterArrays as if they were arrays of simple numeric, logical or character values. The `names()` function returns the names of individual elements of a vector-like RasterArray.
```{r demname, echo= TRUE}
names(dems)
```

The `colnames()` and `rownames()` functions are more relevant for matrix-like RasterArrays, such as `clim`.

```{r climane, echo= TRUE}
colnames(clim)
rownames(clim)
```

All name-related methods can be used for replacement as well. For instance, you can quickly rename the names of the columns of the `clim` object this way:


```{r climrename, echo= TRUE}
clim2 <- clim
colnames(clim2) <- c("temp", "prec")
clim2
```

Just as you would do with normal arrays, you can query/rewrite all names with the `dimnames()` function, which uses a list to store the names in every dimension.

```{r climrenamedim, echo= TRUE}
dimnames(clim2)[[1]] <- 1:10
clim2
```

Besides the names of the elements in the RasterArray, every layer has its own name in the stack. These can be accessed with `layers()` function:

```{r layernames, echo= TRUE}
layers(clim)
```

The total number of cells in the RasterLayer or the entire stack can accessed with the `ncell()` and `nvalues()` functions, respectively.

```{r cells, echo= TRUE}
ncell(dems)
nvalues(dems)
```

## 2.3. Subsetting and replacement
Facilitating the accession of items is the primary purpose of RasterArrays. These either focus on the layers (stack items, double bracket operator "[[") or on the elements of the RasterArray (single bracket operator "["). 

### 2.3.1 Layer selection - Double bracket [[
This form of subsetting and replacement is inherited from the RasterStack class. Individual layers can be accessed directly from the stack using either the position index, the name of the layer, or the logical value pointing to the position. Whichever is used, the RasterArray wrapper is omitted and the output will be a RasterLayer or RasterStack class object.

A single layer can be accessed using its name, regardless of its position in the RasterArray. This can be visualized either with the default `plot()` or the more general `mapplot()` function. 

```{r single, echo= TRUE, plot=TRUE, fig.height=5.5}
one <- dems[["dem_45"]]
mapplot(one, col="earth")
```

The double bracket returns a single RasterArray. Using multiple elements for subsetting will return a RasterStack:

```{r dual, echo= TRUE}
dems[[c(1,2)]]
```

These are the first two RasterLayers in the stack of the RasterArray.

Double brackets can also be used for replacements, but as this has no effect on the structure of the array, changes implemented with this method are more difficult to trace. For instance, 

```{r doublerep, echo= TRUE}
# copy
dem2 <- dems
dem2[["dem_0"]] <- dem2[["dem_5"]]
```

will rewrite the values in the first element of `dem2`, but that will not be evident in the RasterArray's structure.

```{r doublecompare, echo= TRUE}
# but these two are now the same
dem2[[1]]
dem2[[2]]
```

### 2.3.2 Single bracket

Features offered by the double bracket ("[[") operator are virtually identical with those of RasterStacks. The true utility of RasterArrays become evident with simple array-type subsetting. 

Unlike Raster* objects of the raster package, single brackets will retrieve and replace items from the RasterArray as if they were simple arrays. For example, single elements of the DEMs can be selected with the age of the DEM, passed as a character subscript.

```{r dem30, echo= TRUE}
dems["30"]
```

returning the 30Ma RasterLayer. By default, the RasterArray container is dropped, but it can be conserved if the drop argument is set to `FALSE`.
```{r demdrop, echo= TRUE}
dem30 <- dems["30", drop=FALSE]
class(dem30)
```

Accessing elements which are not present is valid for single dimensional RasterArrays (vector-like ones):

```{r beyond, echo= TRUE}
dems[4:12]
```

Missing values are legitimate parts of RasterArrays. These gaps in the data are not represented in the stacks but only in the index slots of the RasterArrays. They can be inserted or added into the layers.


```{r demna, echo= TRUE}
demna <- dems
demna[3] <- NA
```

Multidimensional subscripts work in a similar fashion. If a single layer is desired from the RasterArray, that can be accessed using the names of the margins.

```{r cellsbu, echo= TRUE, plot=TRUE, fig.height=5.5}
# character type is necessary as the rowname is "2003"
one <- clim["2003", "bio1"]
mapplot(one)
```

similarly to entire rows,

```{r entirerow, echo= TRUE}
clim["2005", ]
```

or columns

```{r entirecol, echo= TRUE}
clim[,"bio12"]
```

### 2.4 Inherited from Raster*

As the spatial information is contained entirely in the RasterStacks, a number of methods are practically inherited from RasterStack class. For instance, all RasterLayers of the RasterArray can be cropped in a single line of code.

```{r crop, echo= TRUE, plot=TRUE, fig.height=5.5}
# crop to Australia
ext <- extent(c(                
  xmin = 106.58,
  xmax = 157.82,
  ymin = -45.23,
  ymax = 1.14 
)) 

# cropping all DEMS 
au<- crop(dems, ext)

# select the first element
mapplot(au[1], col="earth")
```

If you plot all cropped DEMs, you can see how Australia drifted to its current position.

Other functions such as aggregation or resampling work in a similar way. 

```{r resam, echo= TRUE, plot=TRUE, fig.height=5.5}
template <- raster(res=5)

# resample all DEMS
coarse <- resample(dems, template)

# plot an element
mapplot(coarse["45"], col="earth")
```


# 3. Plotting
The `mapplot()` function makes it easy to produce visually pleasing plots of a Raster* object. 

## 3.1 Colour palettes
The package includes several colour palettes which can be used for plotting purposes. An additional palette option developed specifically with DEMs in mind is *earth*. This combines the *ocean* and *terra* palettes and automatically sets the breaks to differentiate between marine and terrestrial cells. An example is shown below in [RasterLayer plotting example](#RasterLayer).

```{r palette, fig.height=3, fig.width=8}
showPal()
```

## 3.2 RasterLayer {#RasterLayer}
The `mapplot()` function for RasterLayer works similar to the `plot()` function. The `mapplot` function for the Raster* objects include a default palette and omits the legend, axes and bounding box. 

```{r rlayer_plot, echo=TRUE, out.width='80%', fig.align="center"}
data(dems)
mapplot(dems[1])

#using a custom colour palette
mapplot(dems[1], col="earth", main="Using the earth palette")

```

## 3.3 RasterArray

RasterArrays can be plotted as multi-faceted plots using `mapplot()`. The `mapplot()` function keeps the structure of the RasterArray in terms of the order that the plots are generated. The plot titles are automatically generated based on the names of the layers within the RasterArray. This can be changed by passing the custom plot titles to the argument `plot.title`. **Note:** The number of plot titles provided should be equal to the number of layers in the RasterArray. 

```{r rarray_plot_def, echo=TRUE, out.width='80%', fig.align="center"}
data(dems)
mapplot(dems)

```

### 3.3.1 Adding a single consistent legend
Just like with [RasterLayers](#RasterLayer), different palettes can be used for the plots. This implements a single palette for all the plots. In addition, this can then be used to generate a single legend that ensures consistency and makes it easier to compare the figures. 
```{r rarray_plot_leg, echo=TRUE, fig.align="center"}
data(dems)
mapplot(dems, col="ocean", legend=TRUE, legend.title="DEM", 
        plot.title=proxy(dems))
```
### 3.3.2 RasterArrays with NA layers

In the event that one of the layers of the dem does not exist, i.e. has an NA value instead, the plot corresponding to this NA layer will show "Plot not available".
```{r na_layer, echo=TRUE, fig.align="center"}
data(dems)

#create NA layer
dems[5] <- NA
is.na(dems)

mapplot(dems, col="ipccTemp")

```

### 3.3.3 Number of columns and pages

The custom number of columns can also be specified through the `ncol` argument. The default number of columns is 3. 

```{r rarray_plot_ncol, echo=TRUE}
# 4 columns
data(dems)
mapplot(dems, ncol=4)
```

Sometimes, there might be too many plots to plot all on one single page. In this instance, the argument `multi` can be set to `TRUE` to allow the maps to be plotted in multiple figures. 
```{r rarray_plot_multi, echo = TRUE, eval=FALSE}
data(dems)
mapplot(dems, multi=TRUE)
```

### 3.3.4 Multiple variables in RasterArray

The above described methods can also be applied to RasterArrays containing more than one variable. However, in this instance, the number of columns automatically defaults to the number of variables in the array. Also, plot titles for each individual plot are not printed. Instead, row labels corresponding to the row names of each variable are added, and when `legend=TRUE`, the default legend title defaults to the variables names in the RasterArray. 

```{r rarray_plot_nvar, echo=TRUE, fig.width=8, fig.height=4.5}
data(clim)

mapplot(clim[1:2,], legend=TRUE)
```

The row labels and legend titles can be edited using the argument `rowlabels` and `legend.title` respectively. Different colours for each variable can also be provided. The number of provided colour palettes should be either 1 or correspond to the number of variables.  
```{r rarray_plot_rowlabs, echo=TRUE, fig.width=8, fig.height=4.5}
data(clim)

mapplot(clim[1:2,], col=c("ipccTemp", "ipccPrec"), 
        legend=TRUE, legend.title=c("Temperature", "Precipitation"))
```

# References
Hijmans, R. J., & van Etten, J. (2019). raster: Geographic Data Analysis and Modeling. Retrieved from https://cran.r-project.org/package=raster

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., … Kessler, M. (2017a). Data from: Climatologies at high resolution for the earth’s land surface areas. Dryad Digital Repository. https://doi.org/10.5061/dryad.kd1d4

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., … Kessler, M. (2017b). Climatologies at high resolution for the earth’s land surface areas. Scientific Data, 4(1), 170122. https://doi.org/10.1038/sdata.2017.122

Scotese, C. R. Wright, N. (2018). PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic. URL: https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/