Building Geographic Spaces
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.
# 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
#> $gdat
#> xmin xmax ymin ymax ncell
#> -9.200481e+01 -8.800481e+01 -2.001090e+00 9.989105e-01 1.200000e+05
#> ncol hex
#> 4.000000e+02 NA
#>
#> $tdat
#> [1] 2.25 1.75 1.25 0.75 0.25 0.00
#>
#> $layers
#> class : SpatRaster
#> size : 300, 400, 5 (nrow, ncol, nlyr)
#> resolution : 0.01, 0.01 (x, y)
#> extent : -92.00481, -88.00481, -2.00109, 0.9989105 (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 (EPSG:4326)
#> source : galapagos.tif
#> names : layer, layer, layer, layer, layer
#> min values : -3949.31958, -3949.176758, -3948.405762, -3947.85498, -3740.720947
#> max values : 1424.979004, 1090.23645, 971.93042, 1169.806519, 1634.660767
#>
#> attr(,"class")
#> [1] "geoglist"
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.
# 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.
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:
# 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:
# 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:
# 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:
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:
# 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:
# minimum spanning tree linkage, then plot
rst <- link_islands(rst, klink = 1)
#> Resolving layers [1/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:
# link then plot
rst <- link_islands(rst, klink = 2)
#> Resolving layers [1/1]
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:
# link then plot
rst <- link_islands(rst, klink = NULL)
#> Resolving layers [1/1]
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:
# link then plot
rst <- link_delaunay(rst)
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:
# 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:
# 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)
#> Resampling layer [1/1] | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
#>
h3 <- rast_to_geoglist(cret, cret_l, as.hex = T, hex = 3)
#> Resampling layer [1/1] | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%
#>
# plot to compare
plot(h2)
Islands in hexagonal grids can be linked in precisely the same way as
shown above for raster grids. Let’s add the Voronoi neighbourhood:
# voronoi linkage (klink is NULL by default)
h3 <- link_islands(h3)
#> Resolving layers [1/1]
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:
# 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:
rtd <- build_tardis(rst)
#> Linking geog [1/1]
#> Building graph
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:
head(rtd$edges)
#> from to type bearing hdist vdist gdist
#> [1,] 14523 14524 0 89.99994 1113.127 41.86443 1113.914
#> [2,] 14524 14523 0 -89.99994 1113.127 -41.86443 1113.914
#> [3,] 14523 14923 0 180.00000 1105.744 32.61698 1106.225
#> [4,] 14923 14523 0 0.00000 1105.744 -32.61698 1106.225
#> [5,] 14523 14924 0 134.80927 1568.988 80.35590 1571.045
#> [6,] 14924 14523 0 -45.19062 1568.988 -80.35590 1571.045
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
# 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
#> [1] 1 1 1 1 1 1 1 1 1 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"]
#> [1] 26556 44066 45278 76453 65298 77047 60567 38054 75447 63369
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:
head(rtd$tgraph$dict)
#> ref id
#> 1 14523 0
#> 2 14524 1
#> 3 14923 2
#> 4 14924 3
#> 5 14525 4
#> 6 14925 5
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:
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:
# 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)
#> Linking geog [1/5]Linking geog [2/5]Linking geog [3/5]Linking geog [4/5]Linking geog [5/5]
#> Building graph
# view the temporal information
rtd$tdat
#> [1] 2.25 1.75 1.25 0.75 0.25 0.00
rtd$tlink
#> [1] 1
# view the link types
unique(rtd$edges[,"type"])
#> [1] 0 2
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:
# 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,]
#> smp tmp spt x y
#> [1,] 289881 3 49881 -89.19981 -0.2460895
#> [2,] 248763 3 8763 -88.37981 0.7839105
#> [3,] 107603 1 107603 -91.97981 -1.6960895
#> [4,] 433697 4 73697 -91.03981 -0.8460895
#> [5,] 153506 2 33506 -88.94981 0.1639105
#> [6,] 398951 4 38951 -90.49981 0.0239105
#> [7,] 546612 5 66612 -89.88981 -0.6660895
#> [8,] 533765 5 53765 -90.35981 -0.3460895
#> [9,] 268807 3 28807 -91.93981 0.2739105
#> [10,] 586754 5 106754 -88.46981 -1.6660895
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.