---
title: "Building Geographic Spaces"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Building Geographic Spaces}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
oldpar <- par(no.readonly = TRUE)
```

```{r setup, include=F}
library(terra)
library(rTARDIS)
```

## The `geoglist` class

`rTARDIS` defined an S3 class called a `geoglist`, built from raster grids to represent a geographic space and its connections. A `geoglist` can be built from a single raster layer, or a time-ordered set of raster layers using the function `rast_to_geoglist()`. 

The structure of a `geoglist` is deliberately simple to maintain interoperability with other spatial packages like `sf` or `stars` and permit easier manipulation by the user. It consists of a list with three default entries. `$gdat` is a named numeric vector which describes the extent (`xmin`, `xmax`, `ymin`, `ymax`), the total possible cells in a single grid layer (`ncell`), the number of columns in the grid (`ncol`), and the H3 grid resolution (`hex`). In a hexagonally gridded `geoglist`, the `hex` slot in `$gdat` will record that resolution, while `ncol` will be `NA`; for a rectangular grid it is the other way around. The second slot, `$tdat`, is `NULL` be default, but can instead records the contiguous temporal boundaries of the grid layers. The final slot, `$layers`, is either a `SpatRaster` if rectangular grids are used, or a `SpatVectorCollection` if re-sampling to hexagonal grids is used.

```{r}
# convert a toy dataset to a geoglist and view its structure
rasts <- rast_to_geoglist(galapagos(), times = c(seq(2.25, 0, -0.5), 0))
rasts
```

Any links for connecting masked regions of geographic space are added as a fourth element of the `geoglist` (`$links`). In-keeping with the package architecture, this is a `SpatVector` which records the start cell ID (`$srt`), end cell ID (`$end`), parent layer (`$layer`) and great circle distance in metres (`$distance`) of each link. By default, all links generated by `rTARDIS` functions are symmetric (i.e., srt:end = 1:2 and srt:end = 2:1 will be added to a `geoglist`). As with the rest of the `geoglist`, this format has been kept deliberately simple to enable user manipulation and interoperability with other packages. 

`rTARDIS` has a custom plotting method for the `geoglist` class to display a specified layer, with options to customise plot limits, additional data (legends, axes, links), and aesthetics.

```{r, fig.width=9.5, fig.height=10, dev="svg"}
# plot the first layer of the geoglist above
plot(rasts, 1)
```

`rTARDIS` also provides a function to display global rasters on a 3D interactive sphere, based on the `rgl` plotting library. By default, all plotting functions will plot the first layer in a `geoglist`.

```{r, eval=F}
plot_globe(rasts)
```

## Masking

Masking is fundamental to distinguish between accessible and inaccessible geographic space in a raster. `rTARDIS` itself does not perform this and instead the user is directed to the `terra` function `mask()`. The function offers a few different ways of doing this. To be compatible with `rast_to_geoglist()`, however, unmasked regions must have cell values of 1 and masked regions must have their cell values set to `NA`. A couple of quick 'recipes' for creating masks for different geographic spaces are given below, using a dataset from `rTARDIS` of the topography and bathymetry of the Galapagos archipelago through time.\

### Land surface

One likely case use of `rTARDIS` is to examine movement of terrestrial species across fragmented landmasses, whether at the scale of continents or island systems. A land-sea mask is easily created as follows:

```{r, fig.width=9.5, fig.height=5, dev="svg"}
# load a dataset of the Galapagos archipelago and select the present day layer
gal <- galapagos()[[5]]

# create a land-sea mask from the archipelago raster set
gal_l <- classify(gal, matrix(c(-Inf, 0, NA, 0, Inf, 1), ncol = 3, byrow = T), right = F)

# plot original data and the land-sea mask
par(mfrow = c(1, 2))
plot(gal)
plot(gal_l)
```

### Shallow marine area

Alternatively, you might want to examine the movement of a marine species which can only occupy fairly shallow water conditions. Logically, we cannot allow it to access terrestrial areas, but neither do we want to include deeper ocean areas:

```{r, fig.width=9.5, fig.height=5, dev="svg"}
# create a shallow sea mask at -200 to 0 metres (the photic zone) from the archipelago raster set
gal_m <- classify(gal, matrix(c(-Inf, -200, NA, -200, 0, 1, 0, Inf, NA), ncol = 3, byrow = T), right = F)

# plot original data and the land-sea mask
par(mfrow = c(1, 2))
plot(gal)
plot(gal_m)
```

In the above plot, we can see that there are some parts of the coast missing. Logically, there must be shallow water here, but the resolution of the raster means than the average depth falls below our threshold. We can rectify this with a handy function for detecting cells at the boundaries of cell patches, then combine this with our bathymetric mask:

```{r, fig.width=9.5, fig.height=5, dev="svg"}
# create a land mask
land <- classify(gal, matrix(c(-Inf, 0, NA, 0, Inf, 1), ncol = 3, byrow = T), right = F)

# get all cells around the boundaries of land cell areas and find their indices
coast <- boundaries(land)
coast <- which(coast[] == 1)

# add those cell indices to our previous shallow sea mask
gal_m[coast] <- 1

# plot original data and the land-sea mask
par(mfrow = c(1, 2))
plot(gal)
plot(gal_m)
```

Much better! The mask now permits shallow marine circumnavigation of all islands. While topography and bathymetry are expected to be the primary basis for distinguishing between masked and unmasked areas, rasters with other sorts of data could be used, for example of mean annual temperature where all areas below a particular threshold for habitability are masked.

To implement a mask, we simply supply the input raster and mask to `rast_to_geoglist()`. In this instance, we will also define the temporal durations of each layer. Note that the input raster cannot contain missing values. If your dataset contains such values, then these should be infilled, potentially with an average of the surrounding cells, or a dummy value, and then excluded using the mask if desired. We will use the land mask from earlier:

```{r, fig.width=9.5, fig.height=10, dev="svg"}
rst <- rast_to_geoglist(gal, mask = gal_l)
plot(rst)
```

## Adding links

There may be many reasons for choosing to link islands. In some cases, empirical evidence may show that traversal between isolated areas of geographic space genuinely occurred. In other cases, we may want to compare the impact of different linkage structure on dispersal distances and routing solutions. Regardless of the motivation, TARDIS offers several core methods for link generation guided by geometry and a couple of simple, parsimonious assumptions.\

### Minimum spanning trees

The first method uses a Voronoi tesselation to establish 'domains' of influence that extend outwards from island areas. If two domains share a border, their islands are treated as 'in-line-of-sight' of one another and are candidates for linkage. Functions for visualising islands and their domains are available in `rTARDIS`. We will use the `geoglist` from earlier:

```{r, fig.width=9.5, fig.height=10, dev="svg"}
# detect regions (bounds.only = F means that all island cells are returned, rather than just their margins)
rg <- regions(rst, bounds.only = F)

# get the voronoi domains for those regions
dm <- domains(rst)

# plot together, cropping to the voronoi domain extent
plot(rst, xlim = ext(dm)[1:2], ylim = ext(dm)[3:4])
plot(dm, add = T)
```

The assumption is then made that movement between two islands will be easiest between their geographically closest points, designating the link position. The most conservative way of linking the islands such that any island can be reached from any other island is with a minimum spanning tree. Each island is connected to its nearest neighbour based on the set of Voronoi-derived links. Linked sets of islands are then progressively linked to their other nearest linked sets until the minimum spanning tree is complete. This is implemented in the function `link_islands()` with the argument `klink = 1`:

```{r, fig.width=9.5, fig.height=10, dev="svg"}
# minimum spanning tree linkage, then plot
rst <- link_islands(rst, klink = 1)
plot(rst, xlim = ext(dm)[1:2], ylim = ext(dm)[3:4])
```

### k-nearest neighbours

Many of the islands in the Galapagos are in line of sight of more than one island, with not much extra distance to their next-nearest neighbours. As such, we might decide that minimum spanning tree linkage is too restrictive. Instead, we can repeat the progressive linkage procedure with higher values of `klink`. Higher-degree linkage between islands is subject to the number of neighbours in the Voronoi diagram. Let's try with `klink = 2` and see what the result looks like:

```{r, fig.width=9.5, fig.height=10, dev="svg"}
# link then plot
rst <- link_islands(rst, klink = 2)
plot(rst, xlim = ext(dm)[1:2], ylim = ext(dm)[3:4])
```

Our neighbourhood of islands is now better connected, although some close islands are still not directly connected. The Voronoi diagram itself represents the highest degree of linkage possible when only single links between islands are permitted and can be returned by setting `klink = NULL`. This time we will replot with the voronoi domains to show how they relate to the full set of links:

```{r, fig.width=9.5, fig.height=10, dev="svg"}
# link then plot
rst <- link_islands(rst, klink = NULL)
plot(rst, xlim = ext(dm)[1:2], ylim = ext(dm)[3:4], lcol = "gold", lwd = 2)
plot(dm, add = T, border = "grey")
```

### Delaunay linkage

While full Voronoi linkage creates a more densely connected neighbourhood, it only produces single connections between islands, which may still be quite restrictive. Besides multiple islands lying at potentially similar distances to one another, there may be multiple, similarly spaced points between the boundaries of neighbouring islands between which linkage would be reasonable. The final core method of linkage uses a different function, `link_delaunay()` to connect individual boundary cells to their neighbours, based on Delaunay triangulation:

```{r, eval = F}
# link then plot
rst <- link_delaunay(rst)
plot(rst)
```

```{r, echo = F, fig.width=9.5, fig.height=10, dev="svg"}
# link then plot
rst$links <- readRDS(system.file("extdata", "dlinks.RData", package = "rTARDIS"))
plot(rst)
```

Once we have added links, we can re-use the `regions()` function to determine which islands are linked together by setting `use.links = T`:

```{r, fig.width=9.5, fig.height=10, dev="svg"}
# link then plot
regs <- regions(rst, use.links = T, bounds.only = F)
plot(regs)
```

This will let us know if a layer still contains isolated regions (in this case none). This could be desirable for geographic realism, although they may cause problems for downstream functions and the user should consider masking them. The output of `regions()` itself could be used for this task.\

### Interactive linkage

Beyond these core methods, `rTARDIS` offers a function to build links interactively, `click_to_link()`. This function will plot the selected layer in the RStudio window and allow the user to add `nlink` new links by clicking on pairs of cells to connect. All clicks are resolved to the nearest available boundary cells. The inverse operation is also possible. `click_to_snip()` will plot the map and links, and the user can then click either side of `nsnip` links to remove them. You can try this out with the `rst` object created above.

Ultimately, linkage in this way will not necessarily be guided by geometric principles. It is therefore down the the user to justify the creation of links in this way. It may be useful, however, for manually creating linkage scenarios that fall between between full Voronoi linkage and Delaunay linage in terms of connectivity.\

## Hexagonal gridding

Thus far, we have built our models of geographic space directly from a rectangular raster grid. Another, key feature of `rTARDIS`, however, is the potential to resample from a rectangular grid to a hexagonal grid, supported by Uber's H3 library and the interface provided by the `h3jsr` R package. Hexagonal grids can be used at any spatial scale (subject to their tessellation resolution), but are anticipated to be most useful when constructing geographic spaces at global scales in order to minimise grid distortions towards the poles. 

Hexagonal gridding is implemented in `rast_to_geoglist()` by setting `as.hex = T`. By default, the function will choose the H3 resolution level which most closely matches the resolution of the input raster (`hex = "auto"`). In practice, this grid could be a little finer than the input and result in quite a large grid. As such, you can alternatively select the resolution directly (e.g., `hex = 6`). We will use a different dataset this time, of the topography and bathymetry of the Earth in the late Cretaceous:

```{r, fig.width=9.5, fig.height=4.5, dev="svg"}
# load the dataset and select the first time slice
cret <- cretaceous()[[1]]

# create the land-sea mask
cret_l <- classify(cret, matrix(c(-Inf, 0, NA, 0, Inf, 1), ncol = 3, byrow = T), right = F)

# convert to geoglist with h = 2, and h = 3
h2 <- rast_to_geoglist(cret, cret_l, as.hex = T, hex = 2)
h3 <- rast_to_geoglist(cret, cret_l, as.hex = T, hex = 3)

# plot to compare
plot(h2)
plot(h3)
```

Islands in hexagonal grids can be linked in precisely the same way as shown above for raster grids. Let's add the Voronoi neighbourhood:

```{r}
# voronoi linkage (klink is NULL by default)
h3 <- link_islands(h3)
```

`rTARDIS` also offers a function to transform geoglists to projections other than longitude-latitude. This is primarily intended for the sake of making plots as other `rTARDIS` functions will not accept geoglists that are not longitude-latitude. As we are dealing with a global map, we will use the Equal Earth projection:

```{r, fig.width=9.5, fig.height=10, dev="svg"}
# project to Equal Earth (you can look up proj strings for other map types online)
heq <- project_geoglist(h3, crs = "+proj=eqearth")

# plot
plot(heq, bg = "aliceblue")
```

Looking pretty good! As we can see, North America and Eurasia remain connected across the antimeridian. As such, the only route into South America is via Africa, even though what will become the Gulf of Mexico is also somewhat close to the northern coast of South America. This is an instance where, if fossil evidence supported it, we might want to manually add a connection using `click_to_link()`. We might also decide that, while geometrically permissible in the Voronoi diagram, the connections between India (newly rifted from Africa) and island chains near China are unrealistically long. These could be manually deleted using `click_to_snip`. Feel free to try these operations out.\

## The `tardis` class

Once we have our geographic spaces masked and linked, we are ready to convert them to a graph ready for analysis. This is implemented using the function `build_tardis()`, using a `geoglist` that we have prepared, and returns an S3 class, `tardis`. Crucially, this function assumes that the input `geoglist` directly records physical geographic space, i.e., topography and bathymetry in metres. As such, using rasters in other values will cause funky and likely erroneous results. In the simplest case, we can build this object from a single geographic layer. We will reuse our Galapagos data from earlier:

```{r}
rtd <- build_tardis(rst)
```

As with the `geoglist`, a `tardis` object is deliberately simple to enable manipulation and interoperability. It is a list with several elements. `$gdat` records the same landscape properties as the input `geoglist`. `$tdat` and `$tlink` are both `NULL` in this case as we have only one landscape layer and so no temporal data to worry about (we will build a multi-layer `tardis` shortly). `$tgraph` is the most contrived component; it is a list of graph attributes designed for internal use elsewhere in `rTARDIS` where functions from the `cppRouting` package are called in order to eliminate the processing overhead of constructing this data from scratch. The main element of the `tardis` class, however is `$edges`. This is a matrix which records the base physical attributes of the landscape it represents:

```{r}
head(rtd$edges)
```

The matrix is an edgelist, describing the connections between cells and the properties of those edges. Columns `from` and `to` a numeric indexes that correspond to cells in the input `geoglist`. `type` denotes the link type: 0 for adjacent cells, 1 for links across masked space, and 2 for temporal links (see below). The final four columns describe the edge direction (`bearing`), horizonal distance (`hdist`), elevation change along its length (`vdist`) and the total Pythogorean distance incurred by horizonal and vertical movement (`gdist`). All edges are symmetric by default. If you examine the first two rows you will notice that `from` and `to` are inverted, with the elevation change switching the sign of `vdist` - i.e., uphill versus downhill.\

## Cell indexing

Now that we have seen the structure of a `geoglist` and a `tardis` object, we can understand how `rTARDIS` indexes cells in space and time. For a rectangular grid, the rows and columns (rows * columns) within its extent determines the number of possible cells in a given layer, `ncell`. Any coordinate can be resolved to a row-major numeric cell index with the `terra` function `cellFromXY()`, letting us know its position relative to the total cell count (e.g., cell 234 out of `ncell`). 

Things are a little more complicated for hexagonal grids. Any coordinate can be resolved to an H3 cell at a specified H3 resolution with the `h3jsr` function `point_to_cell()`, but the cell ID itself (e.g., "83dc8cfffffffff") does not tell as anything about its location relative to `ncell` as the hexagonally gridded cells do not display spatially intuitive ordering, unlike the row-major order for rectangular grid cells. Nonetheless, the order of cells for a given extent and H3 resolution is still consistent, allowing us to retrieve a character vector of cell IDs with `ncell` elements, then match an H3 cell ID against this vector to produce a consistent numeric index (e.g., 112 out of `ncell`). All cells in a hexagonally gridded `geoglist` layer follow this same consistent ordering, with masked cells represented by empty geometries and `NA` values to reduce memory footprint, but retaining a structure similar to the `data.frame` of a rectangular raster grid where all grid cells are recorded, including those with `NA` values.

This spatial consistency permits `rTARDIS` to modify cell indexes to embed temporal as well as spatial position. Just as each graph node represents a discrete portion of continuous geographic space, each graph layer represents a discrete interval of continuous time. There can only ever maximally be `1:ncell` non-NA cells in a given spatial layer. Cells `1:ncell` come from the first graph layer, but cells in the second layer are then indexed as `(1:ncell) + (ncell * 1)`. In this way, the first possible grid cell in the subsequent layer has an index of `ncell + 1`. This value falls outside of the range `1:ncell` so we know that it cannot come from the first graph layer. The pattern continues for the third layer as `(1:ncell) + (ncell * 2)` and so on. In this way, the temporal index (layer number) can be retrieved by integer division: `(cellid %/% ncell) + 1`; while the spatial index (grid cell number) can be retrieved by the modulo: `cellid %% ncell`. Some examples are given below, using our new `tardis` object

```{r}
# get a random sample of cell IDs
smp <- sample(rtd$edges[,1], 10)

# get their temporal position (i.e., layer 1 of 1)
smp %/% rtd$gdat["ncell"] + 1

# get their spatial position (raster cell number). As we expect, these values are unchanged as they come from the first (and only) layer in the geoglist
smp %% rtd$gdat["ncell"]
```

Ultimately this system means that nodes in the resulting geospatial graph, which are numbered as sequential integers (e.g., 1, 2, 3), will directly correspond to a discontinuous integer series of spatiotemporal cell indexes (e.g., 2, 315, 6554). We can see this by re-examining the `$tgraph`:

```{r}
head(rtd$tgraph$dict)
```

Here, column `ref` contains the `tardis` cell indices while column `id` is just a 0-indexed integer sequence. Ultimately, the edgelist in a `tgraph` object can be converted to graphs from a variety of other R packages for analysis of routes and distance (e.g., `igraph`, `cppRouting`, `rlemon`), cell indexes retrieved from node number, then resolved back to geographic space and temporal context as described above.\

## Temporal linkage

So far, we has discussed the idea of multiple, time-ordered geographic layers, but we have not yet implemented this in our landscape representations. `rTARDIS` makes this very simple to do with the function `build_tardis()`. All we need is a `geoglist` with multiple layers, and a vector of `nlayers + 1` positive numbers describing the time intervals which each layer covers. These intervals need not be even - the only assumption made is that all intervals are time-contiguous.

Lets revisit the Galapagos dataset in full. This is a set of 5 rasters, representing the topography and bathymetry of the archipelago at 0.5 million year (Mya) intervals from 2 million years ago to the present (i.e. slices centered on 2, 1.5, 1, 0.5 and 0 Mya). These share the same extent, and if we plot a couple together, we can clearly see the effect of the hotspot moving east to west within that frame of reference, forcing up the ocean crust and creating new land surface in an east-west manner:

```{r, echo=FALSE, fig.width=9.5, fig.height=4.5, dev="svg"}
# load data
gal <- galapagos()

# create a land-sea mask from the archipelago raster set
gal_l <- classify(gal, matrix(c(-Inf, 0, NA, 0, Inf, 1), ncol = 3, byrow = T), right = F)

# plot original data and the land-sea mask
par(mfrow = c(1, 2))
plot(gal[[1]])
plot(gal[[3]])
```

As such, geographic locations in this space remain constant, even though the geographic composition of that space is changing through time. By default, `build_tardis()` will assume that the layers in geog are ordered forwards in time and create forwards-in-time links between spatially homologous cells between time layers (e.g., between cell `234` and cell `234 + ncell`). This assumption is implemented with the default argument `tlink = 1`. If desired, this can be reversed (`tlink = 2`) or even be made bidirectional (`tlink = 3`), although there are less clear case uses for time to flow other than past to future.

We will use the multi-layer geoglist created at the start of this vignette. To build our spatiotemporal landscape, we will create our interval vector such that each interval spans 0.25 Mya around the layer mid point ages, aside from the present day which will end at 0 Mya. As time will flow forwards in this landscape and as locations remain spatially constant in the frame of reference, we do not need other arguments:

```{r}
# create our time vector
tms <- c(seq(2.25, 0, -0.5), 0)

# build a tardis graph from the masked landscape
rtd <- build_tardis(rasts)

# view the temporal information
rtd$tdat
rtd$tlink

# view the link types
unique(rtd$edges[,"type"])
```

We can now see that the `tardis` object contains our temporal data and links. If we repeat the cell indexing operations used earlier, we can more clearly see the effect of splitting a cell index into its temporal and spatial components:

```{r}
# consistent seed for the vignette output
set.seed(1)

# sample cells randomly
smp <- sample(rtd$edges[,1], 100)

# get their temporal position (multiple layers represented)
tmp <- smp %/% rtd$gdat["ncell"] + 1

# get their spatial position (raster cell number).
spt <- smp %% rtd$gdat["ncell"]

# convert spatial position to coordinates
crds <- xyFromCell(gal, spt)

# print the first 10 as examples
cbind(smp, tmp, spt, crds)[1:10,]
```

Now that we can parameterise the connectivity of a geographic space, we can think about altering this connectivity according to other spatially varying properties or calculating distances and routes within its structure using other functions from `rTARDIS`. These are covered by other vignettes.\

```{r, include = F}
par(oldpar)
```