In andresgmejiar/lbmech: A Package to Study the Mechanics of Landscape and Behavior
options(encoding = 'UTF-8')
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
echo = TRUE,
encoding = 'UTF-8',
dev = "ragg_png"
)
options('digits' = 3)
In this chapter, I introduce the lbmech package for R's cost distance tools, originally the core of the lbmech package's functionality. I derive cost functions for the movement of an arbitrary animal from first principles of energetics and mechanical physics---not from regressed functions of energetic expenditure with limited applicability across taxa and behavioral modes. This allows us to ensure that the cost model is generalizable, and allows us to separate out three independent types of energy expenditure involved in locomotive transport: (1) kinematic/mechanical work for moving things; (2) kinematic/mechanical work against gravity for lifting things; and (3) basal metabolic expenditure required to keep the animal alive. In the second section. I provide a brief example of the analytical framework using GPS, topographic, and ocean current data from the Illes Balears [@Clementi2021;@Farr2007;@Lera2017].
1.1 Computational Considerations
lbmech is a geospatial package originally designed for least-cost path analysis in R employing differential surfaces representing the cost to move between any adjacent locations on the landscape. It contains additional tools to calculate time- and energy-based costs-of-travel for humans and animals moving across the landscape, including built-in statistical tools to derive the cost functions themselves from GPS data. Moreover, the package contains a number of additional features in development related to estimating potential energetic net primary productivity on the landscape from tabulated data; the demographic implications of adaptive processes; and geostatistical tools to measure inequality. However those functions are independent of the cost-distance functions of lbmech that were developed in response to questions arising from the original purpose and will be discussed in other vignettes.
The general philosophy behind the least-cost aspects of the package is that energetic and least-cost path analyses should always be simple, the most time-taking parts should be done at most once, and ideally that costs should be rooted in empirical reality. In general terms, both the cost-distance functions of lbmech and the deterministic functions of the R library gdistance [@vanEtten2017] provide similar capabilities but with notable differences in computation and ease-of-use. Both employ differential surfaces representing the cost to move between any adjacent locations on the landscape, and both store these costs in ways that enumerate the possible transitions. Both then perform transformations on the cost values before generating a network graph object and using igraphs highly-efficient distances() function [@Csardi2006] to perform distance calculations before outputting a matrix or raster. gdistance stores movement costs as their reciprocal (\textit{conductance} , i.e. $1/\textrm{resistance}$ in a sparse matrix representing every possible transition between two cells on a raster. The use of conductance allows for most cell-to-cell transitions to be zero, since impossible transitions due to distance have a conductance of $1/\infty = 0$. This, in turn, allows the sparse matrix to be relatively small on the disk or memory given that zero values are not stored. However, the use of conductance makes it cumbersome to employ functions where $f(0) \neq 0$ requiring the use of index masking. This in turn encounters integer overflow errors with "large" datasets that are nonetheless of a necessary size for most reasonable purposes. Moreover, the use of sparse matrices requires that the object be recreated every time it is modified, greatly limiting I/O speed and severely affecting complex calculations.
lbmech stores data and performs all linear algebra directly on resistance values using the data.table package [@Dowle2023]. Not only is this more intuitive, but it greatly simplifies the syntax necessary for many types of algebraically simple operations. lbmech stores each possible movement as its own unique row, with entries for a from node, a to node, and either the difference in or raw final and initial values of the raster encountered during the transition. Nodes are named after the coordinates of the raster cell to which they correspond, and are stored as character strings in the form 'x,y'. This allows for (1) in-place modification of objects, greatly increasing processing speed; (2) bidirectional raster analysis describing accumulated costs to and from a node, and (3) additive, nonlinear, and multivariate transformations of large rasters and independent considerations without running into integer overflow limits. Since lbmech largely functions as a wrapper for applying data.table, igraph, and fst functions to terra SpatRaster raster objects, it's easy to generate any arbitrary cost function using data.table syntax.
1.2 The Energetics of Locomotion
In order to move across the terrestrial landscape, two principal types of work must be performed. The first is kinetic work, required to move an object at a particular speed, $K = \frac{1}{2}mv^2$ for an arbitrary particle:
```{=tex}
\begin{equation}
K = \frac{1}{4} m v^2 \frac{\ell_s}{L}
\end{equation}
for human locomotion as derived by @Kuo2007: 654---where $m$ is an object's mass, $v$ it's speed, $\ell_s$ is the average step length and $L$ the average length of a leg; and for quadrupeds $K = m(0.685v + 0.072)$ as derived experimentally by @Heglund1982. LBMech also provides alternate energetic models including one proposed by @Pontzer2023 for any arbitrary animal $K = 8m^{-0.34} + 50\left(1 + \sin(2\theta - 74^\circ)\right)m^{-0.12}$, with $\theta$ in degrees), and one based on a simplified harmonic oscillator based on variation about an average velocity ($K = 2mv^2\gamma$, where $\gamma$ is the fractional difference between the maximum/minimum and average velocities).
For the purposes of simplicity, the rest of the derivation will be using the appropriate parameters for humans according to Kuo but the steps are equivalent for the derivation of quadrupeds. The second type of work is work against gravity, required to lift an object a height $h$ ($U = mgh$), where $g$ is the magnitude of the acceleration due to gravity, and $h$ the height an object is raised. Importantly, this latter form of work cost is incurred only when moving uphill. Summing the two components results in the total amount of work performed to move across the landscape:
```{=tex}
\begin{equation}
\tag{1}
W = \frac{1}{4} m v^2 \frac{\ell_s}{L} + \begin{cases}
mgh & dz \geq 0 \\
0 & dz < 0 \\\end{cases}
\end{equation}
To obtain the expected velocity when moving across a landscape of variable topography, we can employ Tobler's Hiking function, which predicts that velocity decreases exponentially as the slope departs from an 'optimum' slope at which an traveler's velocity can be maximized:
```{=tex}
\begin{equation}
\tag{2}
\frac{d\ell}{dt} = v_{\textrm{max}} e^{-k | \frac{dz}{d\ell} - s |}
\end{equation}
where $\ell$ is horizontal movement, $z$ vertical movement, $v_{\textrm{max}}$ the maximum walking speed, $s$ the slope of maximum speed and $k$ how sensitive speed variation is to slope. Canonical applications of the Tobler function for humans employ $v_{\textrm{max}} = 1.5$ m/s, $k = 3.5$, and $s = - 0.05$ = $\tan(-2.83^\circ)$. However, the `getVelocity()` function in the `lbmech` package allows researchers to obtain the appropriate variables for their subject species of interest from high-resolution GPS tracking data.
Humans and animals, however, are far from ideal engines. Generally, for every calorie of work expended humans require five calories of dietary energy resulting in an efficiency factor of $\varepsilon = 0.2$. Dividing Eq. 1 by $\varepsilon$, and plugging in Tobler's Hiking Function (Eq. 2) for $v$ results in the total amount of energy used by a human to travel between two points on the landscape. This is our function for kinematic energy expenditure.
People also consume energy just to stay alive as defined by our base metabolic rate (BMR, represented by $\frac{dH}{dt}$). Assuming a relatively constant amount of daily physical activity and by extension a relatively constant basal metabolic rate, we can model the total amount of energy used by a human to travel between two points on the landscape and stay alive, by adding $H$ to our equation:
```{=tex}
\begin{equation}
\tag{3}
E = \frac{1}{4} m v^2 \frac{\ell_s}{L} + \frac{dH}{dt} dt + \begin{cases} mgh & dz \geq 0 \\ 0 & dz < 0 \\\end{cases}
\end{equation}
Both of the energetic equations derived above display anisotropic behavior due to dependence on $\theta = \arctan \left(\frac{dz}{dx}\right)$. However, only the potential energy term is linearly anisotropic whereas the kinetic energy term is dependent on slope through an exponential term to the second power. As such, the ArcGIS Distance toolkit is insufficient for true energetic analyses. This incapacity is largely because it---like other least-cost tools perform the calculations directly upon the localized value at a particular cell. Rather, since this is a problem of performing a line integral over a differential surface, an algorithm that performs its calculations based on the difference between individual cells is required.
The figure below plots each of the above cost functions (for time, kinetic energy, gravitational potential energy, mechanical/kinematic energy, and net metabolic energy) versus slope, assuming the canonical parameters for Tobler's Hiking Function and a human male. $k = 3.5$, $s = -0.05$ $m = 68$ kg, $v_\textrm{max} = 1.5$ m/s, $L = 0.75$ m, $\ell_s = 1.6$ m, $g = 9.81$ m/s$^2$, $\frac{dH}{dt} = 93$ J/s, and $\varepsilon = 0.2$; (a) Tobler's Hiking Function inverted for Pace $t$; (b) Kinetic Energy, $K$; (c) Work against Gravitational Potential Energy, $U$; (e) Kinematic/Mechanical Energy Performed, $W = K + U$; (f) Net Metabolic Energy Consumed, $E = \frac{W}{\varepsilon} + \frac{dH}{dt}t$.
library(data.table)
library(ggplot2)
funs <- data.table(Slope = seq(-0.25,0.25,length.out = 1000)
)[, `(a) Pace (s/m)` := 1.5 * exp(3.5 * abs(Slope + 0.05))
][, `(b) Kinteic Energy (J/m)` := 1 / 2 * 68 * `(a) Pace (s/m)`^-2
][, `(c) Potential Energy (J/m)` := fifelse(Slope > 0, 68 * 9.81 * Slope, 0),
][, `(d) Mechanical Work (J/m)` := `(b) Kinteic Energy (J/m)` + `(c) Potential Energy (J/m)`
][, `(e) Net Metabolic Work (J/m)` := `(d) Mechanical Work (J/m)`/0.2 + 72 * `(a) Pace (s/m)`]
f <- melt(funs,id.vars = 'Slope', variable.name = 'Cost')
ggplot(f, aes(x = Slope, y = value, color = Cost)) +
geom_line(show.legend=FALSE, linewidth = 1) + facet_wrap(Cost ~ ., scales = 'free', ncol = 2) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5, size = 13, face = 'bold'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.text.x = element_text(size = 12),
strip.background =element_rect(fill=NA,linewidth=NA)) +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
ggtitle("Cost versus Slope using Canonical Parameters\nFor Various Costs") + ylab("Cost")
A least-cost path algorithm is going to preferentially select paths whose segments tend towards the lowest values on these plots. When minimizing the amount of time (Fig. 1a, and considering the canonical parameters), slopes closer towards the slope of fastest travel $s = - 0.05$ are preferred, with costs increasing exponentially as one moves away from that value. When minimizing the energy of motion (kinetic energy, Fig. 1b), the most extreme slopes will be preferred since movement is slowest and kinetic energy is minimized. When minimizing the work performed against gravitational potential energy (Fig. 1c), any negative slope is equally preferred, with positive slopes increasing in cost linearly. When minimizing the kinetic energy, steeper slopes are preferred downhill mildly over flat-ish slopes, with a notable avoidance of slopes tending around $s = -0.05$, with slopes linearly---but rapidly---increasing in positive slopes. In other words, downhill movement is primarily kinetic energy-limited, whereas uphill movement is mainly gravitational potential energy-limited. The situation is similar for net metabolic energy, but here additional penalties are incurred for overly-slow routes that \textit{may} be more mechanically efficient but take long enough that metabolic processes incur costs. Thus, the most efficient slopes are those that occur at the global minimum at $\frac{dz}{d\ell} = -0.0818 = \arctan(-4.7^\circ)$.
One final consideration that to date no other least-cost toolkit rigorously assesses is potential movement on and across water. Movement considerations across bodies of water are conceptually simpler than movement across land. The water's surface can be represented by a two-dimensional vector field with values at each point representing water's surface velocity $\vec{v}\textrm{water}(x,y)$, with still water having a velocity of $\vec{v}\textrm{water} = 0$. The velocity of someone traveling on water then is simply:
```{=tex}
\begin{equation}
\tag{4}
\vec{v}(x,y) = \vec{v}\textrm{water} + \vec{v}\textrm{travel},
\end{equation}
where $\vec{v}_\textrm{travel}(x,y) = v_\textrm{travel} \hat{\ell}$ is the travel speed of the method of travel (such as swimming or rowing a boat) over still water, and $\hat{\ell}$ the intended direction of travel. If $\vec{v} < 0$, then travel in that direction is impossible. For the purposes of simplicity, we assume that the caloric cost is proportional to $v_\textrm{travel}^2$.
## 1.3 Estimating Velocity
The GPS dataset consists of 15,370 usable GPX tracks obtained from <https://wikiloc.com>, and was first described by @Lera2017 who used it to identify seasonal variability in hiking trail activity and network organization. Given a directory with `.gpx` files, the tracks can be imported as a data.table using `lbmech::importGPX()`:
```r
dir = tempdir()
getData('baleares-gps', dir = dir)
# Detect files with 'gpx' extension in working directory
gpx <- list.files(dir, recursive = TRUE, full.names = TRUE, pattern = ".gpx$")
name <- stringr::str_extract(gpx,'(?<=/)[A-Za-z0-9]+(?=.gpx)')
# Import with lbmech::importGPX(), this will take a while
gpx <- lapply(gpx,importGPX)
names(gpx) <- name
# Join into a single object
gpx <- rbindlist(gpx,idcol = TRUE)
Since GPS tracks can be sampled at a per-length rate smaller than the size of even the smallest raster pixels we might employ, GPS tracks are downsampled to a comparable rate---about a pixel length per GPS sample. For an expected pixel size of 50 m and an approximate maximum speed of 1.5 m/s, t_step = 50/1.5. We can use the lbmech::downsampleXYZ() function for this, applying it on all GPS observations with an associated time stamp (filtering for !is.na(t)):
gpx <- downsampleXYZ(gpx[!is.na(t)], t_step = 50/1.5,
t = 't', x = 'long', y = 'lat', z = 'z',
ID = '.id')
library(lbmech)
dir <- "D:/lbmech/Altar"
rd <- "D:/lbmech/Altar/Baleares"
gpx <- fst::read_fst(paste0(rd,"/gpx_filt.fst"),as.data.table=TRUE)
The xyz data is then passed to lbmech::getVelocity(), which for each sequence of GPS points (1) calculates the changes in elevation $\frac{dz}{d\ell}$ and planimetric speed $\frac{d\ell}{dt}$, and (2) performs a nonlinear quantile regression
to get a function of the form Tobler, which returns a list with the nonlinear model, the model parameters, and the transformed data.
velocity_gps <- getVelocity(data = gpx[y < 90], # Filter for y,
degs = TRUE, # Lat/Long; geodesic correction
tau_vmax = 0.95, # Quantile for v_max
tau_nlrq = 0.50, # Quantile for nlrq
v_lim = 3) # Filter for dl/dt
print(velocity_gps[1:(length(velocity_gps)-1)])
Finally, the lbmech::plotVelocity() function can be used to plot the log-transformed probability of a given speed at a given slope versus the regressed function:
plotVelocity(velocity_gps)
1.4 Energetic Cost-Distance Analysis
1.4.1 Procuring Data
Given the large amounts of data downloaded, processed, and re-used at different parts of the workflow, it is highly recommended to define a consistent directory for the workflow particularly during model-building stages.
# Define a working directory
rd <- "Baleares"
if (!dir.exists(rd)){
dir.create(rd)
}
lbmech can automatically download terrain data, but we also want to consider ocean movement. Here, we'll import the Mediterranean Sea Physics Analysis and Forecast (CMEMS MED-Currents, EAS6 system) [@Clementi2021] for ocean currents around two of the principal Balearic Islands---Mallorca and Menorca---on 13 June 2022 and define our region of maximum movement (our world) around that raster's extent:
# Projection will be UTM 31N WGS1984
proj <- "EPSG:32631"
# Import an ocean water surface velocity dataset
currents <- project(getData('baleares-currents'),
proj)
# Define region of maximum possible movement
region <- as.polygons(ext(currents),proj)
library(lbmech)
library(data.table)
library(ggplot2)
library(tidyterra)
library(ggspatial)
library(pandoc)
library(png)
library(grid)
library(magick)
library(cowplot)
library(ggquiver)
setwd("D:/lbmech/Altar/vignettes")
rd <- "D:/lbmech/Altar/Baleares"
# Projection will be UTM 31N WGS1984
proj <- "EPSG:32631"
# Import Ocean current data
currents <- rast(paste0(rd,"/Ocean_Currents.nc"))
currents <- project(currents,proj)
currentU_vals <- values(currents[[1]])
currentV_vals <- values(currents[[2]])
i <- !is.na(currentU_vals)
coords <- as.data.table(xyFromCell(currents, 1:ncell(currents)))
coords <- coords[as.logical(i)]
currentU_vals <- currentU_vals[as.logical(i)]
currentV_vals <- currentV_vals[as.logical(i)]
tpsU <- fields::Tps(coords, currentU_vals)
tpsV <- fields::Tps(coords, currentV_vals)
currents[[1]] <- terra::interpolate(rast(currents[[1]]),tpsU)
currents[[2]] <- terra::interpolate(rast(currents[[2]]),tpsV)
# Define region of maximum possible movement
region <- as.polygons(ext(currents),proj)
region$ID <- 1
# Define the raster properties
z_fix <- fix_z(proj = proj, res = 50)
# Make a grid
grid <- makeGrid(dem = region, nx = 15, ny = 15,
sources = TRUE, zoom = 11,
overlap = 0.05)
# Import list of locations on Mallorca, Menorca and Cabrera
pobles <- getData('baleares-places')
pobles <- project(pobles,proj)
muns <- paste0(rd,"/Borders.shp")
if (!file.exists(muns)){
download.file("https://stacks.stanford.edu/file/druid:sp807yd2954/data.zip?download=true",
paste0(tempdir(),"/SpainBorders.zip"),mode='wb')
unzip(paste0(tempdir()),overwrite=TRUE,junkpaths = TRUE)
muns <- list.files(paste0(tempdir(),"/SpainBorders"),pattern=".shp$",full.names = TRUE)
muns <- vect(muns)
muns <- project(muns,proj)
muns <- crop(muns, currents)
writeVector(muns,outline)
}
muns <- vect(muns)
plotreg <- as.polygons(ext(buffer(muns,12500)),crs = crs(muns))
insetExt <- sf::st_read("./Outputs/Baleares_Region.kml")
insetExt <- project(vect(insetExt),proj)
islands <- aggregate(muns)
countries <- list.files("D:/Census/Global_Divisions/",pattern="_0.shp",full.names=TRUE)
countries <- lapply(countries,vect)
countries <- do.call(rbind,countries)
countries <- project(countries,proj)
countries <- crop(countries,ext(insetExt))
states <- list.files("D:/Census/Global_Divisions/",pattern="_1.shp",full.names=TRUE)
states <- lapply(states,vect)
states <- do.call(rbind,states)
states <- project(states,proj)
states <- crop(states,ext(insetExt))
dem <- rast(elevatr::get_elev_raster(raster::raster(project(currents,"+proj=longlat")), z = 9))
dem[dem < 0] <- NA
dem <- project(dem,proj)
dem <- crop(dem,plotreg)
slope <- terrain(dem,v='slope',unit='radians')
aspect <- terrain(dem,v='aspect',unit='radians')
# Calculate a hillshade raster
hill <- shade(slope,aspect)
names(hill) <- "shades"
# Get a pallete of greys for hillshading
pal_greys <- hcl.colors(1000, "Grays")
# Round the hillshade values to the nearest factor
hill <- (hill - minmax(hill)['min',]) / (minmax(hill)['max',] - minmax(hill)['min',])
hill <- ceiling(hill * (length(pal_greys)-1)) + 1
# Get colors for the elevation DEM subplot
grad_hypso <- hypso.colors2(10, "dem_poster")
names(dem) <- "Elevation"
# Transform currents into data.table to turn into vector arrows
field <- project(currents, fix_z(proj,res=8000),align=TRUE)
field <- crop(field,plotreg)
field <- as.data.table(rastToTable(field))
names(field)[1:2] <- c("u","v")
field <- na.omit(field)
dem[dem == 0] <- NA
hill[is.na(dem)] <- NA
topo <- ggplot() +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, linewidth = 0.75,
col = 'lightblue', inherit.aes=FALSE) +
geom_spatraster(data= dem, na.value = NULL) +
scale_fill_gradientn(" \nElevació\n(m)",
colours = hypso.colors(15,'dem_poster'),
limits = c(0,unlist(global(dem,'max',na.rm=TRUE))),
na.value = "transparent") +
theme_bw() +
theme(legend.position = "top",
legend.text = element_text(angle = 45, size = 15),
legend.title = element_text(size=20),
strip.text = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = '#005EB8')
) +
guides(fill = guide_colorbar(barwidth = 60, barheight = 1.5, alpha = 0.7,
label.position = "top",
override.aes = list(alpha = 0.7))) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4,
maxcell = Inf, na.value = NULL) +
geom_spatvector(data = muns, fill=NA,
size=1,col = '#EEEEEE') +
geom_spatvector(data = islands,
linewidth = 2.5, col = '#454545',fill=NA) +
geom_spatvector(data = pobles, color = 'red', size = 5) +
ggrepel::geom_label_repel(data = pobles, mapping = aes(x = geom(pobles)[,'x'],
y = geom(pobles)[,'y'],
label = pobles$Location),
size = 7) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0)) +
annotation_scale(location = 'br', height = unit(1,'cm'),text_cex = 3) +
annotation_north_arrow(location = 'tl',
style = north_arrow_nautical,
height = unit(4.5,'cm'),
width = unit(4.5,'cm'))
political <- ggplot(states) +
geom_spatvector(fill = '#396b27',
linewidth = 0.5, col = '#D0D0D0') +
geom_spatvector(data = countries,
linewidth = 1, col = 'black', fill= NA) +
geom_spatvector(data = plotreg,
linewidth = 2, col = 'red', fill = NA) +
theme_bw()+
theme(axis.text.y=element_blank(), #remove y axis labels
axis.ticks.y=element_blank(), #remove y axis ticks
axis.text.x=element_blank(), #remove x axis labels
axis.ticks.x=element_blank(),
panel.background = element_rect(fill='#bce6f9'),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
scale_x_continuous(expand=c(0,0)) +
scale_y_continuous(expand = c(0,0))
header <- ggdraw() + draw_image(paste0(rd,"/RegionMap_Header.png"))
side <- plot_grid(header,political,nrow=2, rel_heights = c(0.66,0.33))
gg <- plot_grid(topo,side,ncol=2,rel_widths = c(0.66,0.33))
gg
# pdf(file = "D:/lbmech/Figs/Map.pdf",
# height = 9.5, width =20)
# gg
# dev.off()
1.4.2 Defining a 'world'
As with any raster operation, a consistent projection and grid needs to be defined to snap/fix all sources. However, unlike most raster packages lbmech does not require foreknowledge of the spatial extents; only the resolution and offsets. A raster can be provided, or one can be made with the lbmech::fix_z() function:
# Define the raster properties
z_fix <- fix_z(proj = proj,
res = 50)
Although a world can be defined directly based on a digital elevation model, generally it is easier to define a polygon coincident with the coverage of a digital elevation model and an attribute pointing to a URL/filepath with the source. The lbmech::makeGrid() can make such a polygon for a raster or filepath input, while using a polygon will induce future functions to download SRTM data as-needed. Similar polygons are frequently distributed by state GIS agencies (e.g. Pennsylvania; Massachusetts):
# Make a grid
grid <- makeGrid(dem = region, # Input; here polygon for SRTM
nx = 15, ny = 15, # Cols/Rows to divide the polygon
sources = TRUE, # Just crop/divide, or point to source?
zoom = 11, # Zoom level for SRTM data
overlap = 0.05) # Fraction overlap between adjacent tiles
A world is formally defined as a directory containing all possible constraints on movement, fixed to a particular spatial projection. Given the very large areas and relatively-small spatial resolutions often employed, the data can rapidly become too large to compute all at once in the memory. The defineWorld() function divides the potential world of movement into individual, slightly-overlapping sectors that are calculated and read-in only as needed. It itself generally is not used to perform calculations.
# Resample currents so they're at a more realistic resolution
currents_high <- project(currents,
disagg(currents, ceiling(res(currents)/50)),
method = 'cubicspline')
# Define world of motion within a workspace
defineWorld(source = grid, # Elevation source, like makeGrid() output
grid = grid, # How to divide up the world
directions = 8, # How adjacency between cells are defined
neighbor_distance = 10, # Overlap between tiles in addition to ^
cut_slope = 0.5, # Maximum traversible slope
water = currents_high, # Water velocity source
priority = 'land', # If data for land or water, who wins
z_min = 0, # Minimum elevation, below NA in land
z_fix = z_fix, # Grid with defined projection, resolution
dist = 'karney', # Geodesic correction method
dir = rd, # Working directory
overwrite=FALSE)
1.4.3 Estmiating Costs
Once a world has been defined, a cost function can be applied using the lbmech::calculateCosts() function. The lbmech::energyCosts() function is provided to perform least-time, least-work, and least-energy analyses as in the above sections of this manuscript. As with lbmech::defineWorld(), it itself is not generally used to perform calculations---just to define the values.
calculateCosts(costFUN = energyCosts,
dir = rd, # Working directory with world
method = "kuo", # Method to calculate work
water = TRUE, # Consider water, or only land?
v_max = velocity_gps$vmax, # Max walking speed (m/s)
k = velocity_gps$k, # Slope sensitivity (dimensionless)
s = velocity_gps$s, # Slope of fastest motion (dimlss.)
row_speed = 1.8, # Speed over water (m/s)
row_work = 581, # Work over water (J/s)
m = 68, # Mass (kg)
BMR = 72, # Basal metabolic rate (J/s)
l_s = 1.8, # Stride length (m)
L = 0.8) # Leg length
For a given (set of) cost function(s) and origin(s)/destination(s), getCosts() will identify the needed sectors, make sure the data has been downloaded (and do so if not), perform the calculations, and save the cell-wise transition cost tables to the working directory. That way, if they are needed for another calculation in the future they won't need to be fetched and prepossessed again.
# Import locations on Mallorca, Menorca, and Cabrera
pobles <- project(getData("baleares-places"), region)
reg <- as.polygons(ext(buffer(pobles, 10000)),crs = crs(proj))
costs <- getCosts(region = reg, # Area of maximum movement
from = pobles, # Origins
to = NULL, # Destinations
costname = 'energyCosts', # Name of costFUN above
id = 'Location', # Column with origin names
dir = rd, # Directory with world
output = c('object','file'), # Save files
destination = 'all') # Distance to all points in region
1.4.4 Paths and Corridors
Likewise, least-cost paths can be computed given a set of nodes:
# Find all least-cost paths between Palma de Mallorca and Maó
paths <- getPaths(region = reg,
nodes = pobles,
costname = 'energyCosts',
order = c("Palma de Mallorca","Maó"),
id = 'Location',
dir = rd)
reg <- as.polygons(ext(buffer(pobles, 10000)),crs = crs(proj))
# Import list of locations on Mallorca, Menorca and Cabrera
pobles <- fread("D:/lbmech/Altar/vignettes/localitats.csv")
# Convert to SpatVector
pobles <- vect(pobles, geom=c('long','lat'), crs = '+proj=longlat')
pobles <- project(pobles,proj)
# Uncomment if you've never run this section before before, since the figure
# differs from the example given in the text for clarity
## Get combinations of all paths to and from Palma de Mallorca
# path <- pobles$Location[pobles$Location != 'Palma de Mallorca']
# path <- c(rbind(rep("Palma de Mallorca",length(path)+1),
# path))
# path <- path[1:(length(path)-1)]
## Calculate the paths
# palma_paths <- getPaths(region = reg,
# nodes = pobles,
# costname = 'energyCosts',
# order = path,
# id = 'Location',
# dir = rd)
# Import the least-cost paths
palma_paths <- vect("D:/lbmech/altar/baleares/world/costrasters/PalmaPaths.shp")
# Convert currents to a (smaller) data.table
field <- project(currents, fix_z(proj,res=8000),align=TRUE)
field <- crop(field,reg)
field <- as.data.table(rastToTable(field))
names(field)[1:2] <- c("u","v")
field <- na.omit(field)
# Get even coarser DEM for hillshading
dem <- rast(elevatr::get_elev_raster(raster::raster(project(currents,"+proj=longlat")), z = 8))
dem[dem < 0] <- 0
dem <- project(dem,proj)
dem <- crop(dem,reg)
# Calculate aspect, slope
slope <- terrain(dem,v='slope',unit='radians')
aspect <- terrain(dem,v='aspect',unit='radians')
# Calculate a hillshade raster
hill <- shade(slope,aspect)
names(hill) <- "shades"
# Get a pallette of greys for hillshading
pal_greys <- hcl.colors(1000, "Grays")
# Round the hillshade values to the nearest factor
hill <- (hill - minmax(hill)['min',]) / (minmax(hill)['max',] - minmax(hill)['min',])
hill <- ceiling(hill * (length(pal_greys)-1)) + 1
# Import the cost rsters
costs <- rast(paste0(rd,"/World/CostRasters/Pobles.tif"))
# Keep only the ones involving Palma
palma <- costs[[stringr::str_detect(names(costs),"Palma de Mallorca")]]
# The time ones
palma_dt <- palma[[stringr::str_detect(names(palma),"dt_")]]
# SpatVector Collection with to, from rasters for time
palma_dt_paths <- svc(palma_paths[stringr::str_detect(palma_paths$segment,"Mallorca_to_") &
palma_paths$cost=='dt',],
palma_paths[stringr::str_detect(palma_paths$segment,"_to_Palma") &
palma_paths$cost=='dt',])
# ggplot call for the from raster with time, involving the islands, hillshade
# cost raster, cost paths, and ocean currents
ggtime_from <- ggplot(islands) +
geom_spatraster(data = palma_dt[[1]]/3600) +
scale_fill_viridis_c(expression(paste(atop("Hores",
italic("Hours")," (hr)"))),
option='magma',direction = -1,
limits = rowMeans(minmax(palma_dt))/3600) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = palma_dt[[1]],binwidth = 3600) +
geom_spatvector(data = palma_dt_paths[[1]],
linewidth = 1.4, col = 'yellow') +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "right",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
axis.title.y = element_blank(),
strip.text.x = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()) +
xlab(expression(paste("Temps des de Palma — ",italic("Time from Palma")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0)) +
ylab("")
# Now likewise for time to
ggtime_to <- ggplot(islands) +
geom_spatraster(data = palma_dt[[2]]/3600) +
scale_fill_viridis_c(expression(paste(atop("Hores",
italic("Hours")," (hr)"))),
option='magma',direction = -1,
limits = rowMeans(minmax(palma_dt))/3600) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = palma_dt[[2]],binwidth = 3600) +
geom_spatvector(data = palma_dt_paths[[2]],
linewidth = 1.4, col = 'yellow') +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "none",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
) +
xlab(expression(paste("Temps cap a Palma — ",italic("Time to Palma")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0)) +
ylab("")
# And likewise for work
palma_dW <- palma[[stringr::str_detect(names(palma),"dW_")]]
palma_dW_paths <- svc(palma_paths[stringr::str_detect(palma_paths$segment,"Mallorca_to_") &
palma_paths$cost=='dW_l',],
palma_paths[stringr::str_detect(palma_paths$segment,"_to_Palma") &
palma_paths$cost=='dW_l',])
ggwork_from <- ggplot(islands) +
geom_spatraster(data = palma_dW[[1]]/1000/4.184) +
scale_fill_viridis_c(expression("(kcal)"),
option='magma',direction = -1,
limits = rowMeans(minmax(palma_dW))/1000/4.184) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = palma_dW[[1]],binwidth =500000*4.184) +
geom_spatvector(data = palma_dW_paths[[1]],
linewidth = 1.4, col = 'yellow') +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "right",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
axis.title.y = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) + ylab("") +
xlab(expression(paste("Treball des de Palma — ",italic("Work from Palma")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0))
ggwork_to <- ggplot(islands) +
geom_spatraster(data = palma_dW[[2]]/1000/4.184) +
scale_fill_viridis_c(expression("(kcal)"),
option='magma',direction = -1,
limits = rowMeans(minmax(palma_dW))/1000/4.184) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = palma_dW[[2]],binwidth =500000*4.184) +
geom_spatvector(data = palma_dW_paths[[2]],
linewidth = 1.4, col = 'yellow') +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "none",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
axis.title.y = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) + ylab("") +
xlab(expression(paste("Treball cap a Palma — ",italic("Work to Palma")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0))
# And finally energy
palma_dE <- palma[[stringr::str_detect(names(palma),"dE_")]]
palma_dE_paths <- svc(palma_paths[stringr::str_detect(palma_paths$segment,"Mallorca_to_") &
palma_paths$cost=='dE_l',],
palma_paths[stringr::str_detect(palma_paths$segment,"_to_Palma") &
palma_paths$cost=='dE_l',])
ggenergy_from <- ggplot(islands) +
geom_spatraster(data = palma_dE[[1]]/1000/4.184) +
scale_fill_viridis_c(expression("(kcal)"),
option='magma',direction = -1,
limits = rowMeans(minmax(palma_dE))/1000/4.184) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = palma_dE[[1]],binwidth =1250000*4.184) +
geom_spatvector(data = palma_dE_paths[[1]],
linewidth = 1.4, col = 'yellow') +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "right",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
axis.title.y = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) + ylab("") +
xlab(expression(paste("Energia des de Palma — ",italic("Energy from Palma")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0))
ggenergy_to <- ggplot(islands) +
geom_spatraster(data = palma_dE[[2]]/1000/4.184) +
scale_fill_viridis_c(expression("(kcal)"),
option='magma',direction = -1,
limits = rowMeans(minmax(palma_dE))/1000/4.184) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = palma_dE[[2]],binwidth =1250000*4.184) +
geom_spatvector(data = palma_dE_paths[[2]],
linewidth = 1.4, col = 'yellow') +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "none",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
axis.title.y = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) + ylab("") +
xlab(expression(paste("Energia cap a Palma — ",italic("Energy to Palma")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0))
# Final figure is the above 6 with ggarrange. Can't do it by facets since the scales are all different.
egg::ggarrange(ggtime_to,ggtime_from,
ggwork_to,ggwork_from,
ggenergy_to,ggenergy_from,
nrow=3)
Corridors describe the minimum cost required to take a detour to a given location on a path between two or more given points relative to the least cost path between those points. These can be calculated from the output cost rasters using the lbmech::makeCorridor() function.
corr <- makeCorridor(rasters = costs,
order = rev(c("Castell de Cabrera","Sa Calobra","Maó")))
ggtime_corr <- ggplot(islands) +
geom_spatraster(data = corr[[1]]/3600) +
scale_fill_viridis_c(expression("(hrs)"),
option='magma',direction = -1,
limits = rowMeans(minmax(corr[[1]]))/3600) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = corr[[1]]/3600,binwidth = 0.5) +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "none",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
axis.title.y = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) + ylab("") +
xlab(expression(paste("Temps de desviació — ",italic("Time of detour")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0))
ggwork_corr <- ggplot(islands) +
geom_spatraster(data = corr[[2]]/4.184/1000) +
scale_fill_viridis_c(expression("(kcal)"),
option='magma',direction = -1,
limits = rowMeans(minmax(corr[[2]]/4.184/1000))) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = corr[[2]]/4.184/1000,binwidth = 500) +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "none",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
axis.title.y = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) + ylab("") +
xlab(expression(paste("Treball de desviació — ",italic("Work of detour")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0))
ggenergy_corr <- ggplot(islands) +
geom_spatraster(data = corr[[3]]/4.184/1000) +
scale_fill_viridis_c(expression("(kcal)"),
option='magma',direction = -1,
limits = rowMeans(minmax(corr[[3]]/4.184/1000))) +
geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)],
alpha=0.4, maxcell = Inf) +
geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE,
col = 'white', inherit.aes=FALSE) +
geom_spatvector(
linewidth = 1.25, col = '#000000',fill=NA) +
geom_spatraster_contour(data = corr[[3]]/4.184/1000,binwidth = 1000) +
geom_spatvector(data = pobles,
size = 3, color = 'red') +
theme_bw() +
theme(legend.position = "none",
legend.text = element_text(angle = 45, size = 10),
legend.title = element_text(size=12, hjust=0),
strip.text.x = element_blank(),
axis.title.y = element_blank(),
strip.background =element_rect(fill=NA,linewidth=NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
) + ylab("") +
xlab(expression(paste("Energia de desviació — ",italic("Energy of detour")))) +
scale_x_continuous(position = 'top', expand=c(0,0)) +
scale_y_continuous(position = 'right', expand=c(0,0))
cowplot::plot_grid(ggtime_corr, ggwork_corr, ggenergy_corr,
ncol = 1)
1.5 Influence
Different measures of distance in different topologies are often taken to be measures of influence. In other words, nearer things are more likely to influence an observation than farther things. LBMech provides tools to represent these relationships as various types of neighbor weights (for influence between observations) and Voronoi polygons (for the relative influence of observations everywhere else in space).
1.5.1 Neighbor Weights
Distance matrices can be obtained from the getCosts() function by setting the destination parameter to 'pairs':
dists <- getCosts(region = reg, # Area of maximum movement
from = pobles, # Origins
costname = 'energyCosts', # Name of costFUN above
costFUN = energyCosts,
id = 'Location', # Column with origin names
dir = rd, # Directory with world
destination = 'pairs') # Distance to all points in region
We can try and visually compare the different topologies using the normMDS() function, which performs multidimensional scaling on each distance matrix and rescales/orients them so they are in units of---and directionally relative to---a common path:
uv <- normMDS(dists, # List of distance matrices
xy = pobles, # Cartesian coordinates (optional)
id = 'Location', # Unique ID column
ax = c('Palma de Mallorca','Inca'), # Reference path
chir = 'Sa Calobra') # Third location to preserve chirality
colors <- c(ref = 'black', dt = 'blue', dW_l = 'darkgreen', dE_l = 'darkred')
uvplot <- rbindlist(uv, idcol = TRUE)
ax <- c('Palma de Mallorca','Inca')
axplot <- uvplot[Location %in% ax & .id == 'ref']
uvplot[Location %in% ax & .id != 'ref', Location := '']
ggplot(data = uvplot, mapping = aes(x = u, y = v, label = Location, color = as.factor(.id))) +
scale_color_manual('Distance', values = colors) +
geom_point() +
ggrepel::geom_text_repel(size = 2.5) +
geom_point(data = axplot, mapping = aes(x = u, y = v),
color = 'black', fill = 'yellow', shape = 23, size = 2.5) +
ggtitle('Different implied 2D topologies based on different costs',
subtitle = 'MDS output aligned and scaled such that\nthe paths between Palma and Inca are equivalent') +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5)) +
coord_equal()
lbmech provides makeWeights which can convert any distance matrix (geographic, or cost distance-based) into various types of spatial weights matrices. At present, it supports:
- Cutoff distances and distance-decay formulas for determining group membership
- Any arbitrary distance-decay formula
- Rank-distance in lieu of raw distance
- Tools to deal with singularities arising from distance-decay
- Row standardization (where all rows sum to one)
- Fuzzy row standardization (where entries are a proportion of the largest value of a row)
- Binary weights matrices depending on group membership if provided only a vector of group IDs
weights <- list()
# The three nearest neighbors time-wise are definitely neighbors to an observation
# Everyone else definitely isn't.
weights[[1]] <- makeWeights(dists$dt, bw = 3)
# All observations within 20 MJ of work are neighbors, with row-standardized weights
# based on inverse distance, assuming self-weights are equivalent to 20 MJ
weights[[2]] <- makeWeights(dists$dW_l, bw = 2e7, offset = 2e7,
mode = 'fixed',
weighting = 'distance',
row.stand = TRUE)
# All observations might be neighbors.
# The weights represent the probability that they're neighbors.
# The probability is based on the inverse distance square of the rank of the energy distance
weights[[3]] <- makeWeights(dists$dE_l,
weighting = 'rank',
FUN = function(x) 1/x^2,
row.stand = 'fuzzy')
# Neighbors are those locations on the same island
weights[[4]] <- makeWeights(pobles$Island, ID = pobles$Location)
library(ggraph)
set.seed(43884)
w <- as.matrix(weights[[1]])
w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE)
w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w)))
ggw1 <- ggraph(w, algorithm = 'graphopt', 'igraph') +
geom_point(aes(x = x, y = y)) +
geom_edge_link0() +
ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3)
w <- as.matrix(weights[[2]])
w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE)
w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w)))
ggw2 <- ggraph(w, algorithm = 'kk', 'igraph') +
geom_point(aes(x = x, y = y)) +
geom_edge_link0() +
ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3)
w <- as.matrix(weights[[3]])
w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE)
w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w)))
ggw3 <- ggraph(w, algorithm = 'gem', 'igraph') +
geom_point(aes(x = x, y = y)) +
geom_edge_link0(edge_alpha = igraph::as_data_frame(w)$weight) +
ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3)
w <- as.matrix(weights[[4]])
w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE)
w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w)))
ggw4 <- ggraph(w, algorithm = 'gem', 'igraph') +
geom_point(aes(x = x, y = y)) +
geom_edge_link0() +
ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3)
cowplot::plot_grid(ggw1,ggw2,ggw3,ggw4,ncol=2)
Note that in the top-left subplot, only edges that exist in both directions are plotted.
1.5.2 Catchments and Voronoi Polygons
Voronoi polygons allocate all land to the nearest observation. These are most commonly employed using planar geographic distance as the metric and generally assume equal influence between all observations, such as with the voronoi() function from terra:
v <- voronoi(pobles)
Sometimes the assumption of equal influence does not hold. There exist a number of ways to differentially weight the influence any one point has versus another. Multiplicatively-weighted diagrams are particularly useful when a point's ability to influence space versus a different point is constant. For example, if we wanted to determine the catchment of two towns based on who can get to each point the fastest---with Town 1 possessing a car capable of traveling twice as fast as Town 2---we might consider a Multiplicatively Weighed Voronoi diagram giving Town 1 a weight of $w_1 = 2$ and Town 2 a weight of $w_2 = 1$. In short, Multiplicatively Weighted Voronoi diagrams transform distance to each point such that:
{=tex}
D = w_i * \ell_i
where $D$ is the modified distance, $w_i$ the weight of point $i$, and $\ell$ the raw distance to point $i$. lbmech provides a fully vectorized function to calculate such catchments for planar, spherical, and Earth-geoid topologies through the lbmech::mwVoronoi function:
While lbmech::mwVoronoi is marginally less efficient thant terra::voronoi for equal weighting schemes, it is capable of employing multiplicative weights and performing geodesic corrections.
Likewise, the Voronoi concept can be extended to any arbitrary topology---such as the time-, work- and energy-based distances obtained above. These are often known as 'catchments' in the literature. As long as we have previously calculated the cost rasters for the desired locations and costs, we can obtain the catchments with the lbmech::makeCatchment function:
References
andresgmejiar/lbmech documentation built on Feb. 2, 2025, 12:37 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
options(encoding = 'UTF-8') knitr::opts_chunk$set( collapse = TRUE, comment = "#>", echo = TRUE, encoding = 'UTF-8', dev = "ragg_png" ) options('digits' = 3)
In this chapter, I introduce the lbmech package for R's cost distance tools, originally the core of the lbmech package's functionality. I derive cost functions for the movement of an arbitrary animal from first principles of energetics and mechanical physics---not from regressed functions of energetic expenditure with limited applicability across taxa and behavioral modes. This allows us to ensure that the cost model is generalizable, and allows us to separate out three independent types of energy expenditure involved in locomotive transport: (1) kinematic/mechanical work for moving things; (2) kinematic/mechanical work against gravity for lifting things; and (3) basal metabolic expenditure required to keep the animal alive. In the second section. I provide a brief example of the analytical framework using GPS, topographic, and ocean current data from the Illes Balears [@Clementi2021;@Farr2007;@Lera2017].
1.1 Computational Considerations
lbmech is a geospatial package originally designed for least-cost path analysis in R employing differential surfaces representing the cost to move between any adjacent locations on the landscape. It contains additional tools to calculate time- and energy-based costs-of-travel for humans and animals moving across the landscape, including built-in statistical tools to derive the cost functions themselves from GPS data. Moreover, the package contains a number of additional features in development related to estimating potential energetic net primary productivity on the landscape from tabulated data; the demographic implications of adaptive processes; and geostatistical tools to measure inequality. However those functions are independent of the cost-distance functions of lbmech that were developed in response to questions arising from the original purpose and will be discussed in other vignettes.
The general philosophy behind the least-cost aspects of the package is that energetic and least-cost path analyses should always be simple, the most time-taking parts should be done at most once, and ideally that costs should be rooted in empirical reality. In general terms, both the cost-distance functions of lbmech and the deterministic functions of the R library gdistance [@vanEtten2017] provide similar capabilities but with notable differences in computation and ease-of-use. Both employ differential surfaces representing the cost to move between any adjacent locations on the landscape, and both store these costs in ways that enumerate the possible transitions. Both then perform transformations on the cost values before generating a network graph object and using igraphs highly-efficient distances() function [@Csardi2006] to perform distance calculations before outputting a matrix or raster. gdistance stores movement costs as their reciprocal (\textit{conductance} , i.e. $1/\textrm{resistance}$ in a sparse matrix representing every possible transition between two cells on a raster. The use of conductance allows for most cell-to-cell transitions to be zero, since impossible transitions due to distance have a conductance of $1/\infty = 0$. This, in turn, allows the sparse matrix to be relatively small on the disk or memory given that zero values are not stored. However, the use of conductance makes it cumbersome to employ functions where $f(0) \neq 0$ requiring the use of index masking. This in turn encounters integer overflow errors with "large" datasets that are nonetheless of a necessary size for most reasonable purposes. Moreover, the use of sparse matrices requires that the object be recreated every time it is modified, greatly limiting I/O speed and severely affecting complex calculations.
lbmech stores data and performs all linear algebra directly on resistance values using the data.table package [@Dowle2023]. Not only is this more intuitive, but it greatly simplifies the syntax necessary for many types of algebraically simple operations. lbmech stores each possible movement as its own unique row, with entries for a from node, a to node, and either the difference in or raw final and initial values of the raster encountered during the transition. Nodes are named after the coordinates of the raster cell to which they correspond, and are stored as character strings in the form 'x,y'. This allows for (1) in-place modification of objects, greatly increasing processing speed; (2) bidirectional raster analysis describing accumulated costs to and from a node, and (3) additive, nonlinear, and multivariate transformations of large rasters and independent considerations without running into integer overflow limits. Since lbmech largely functions as a wrapper for applying data.table, igraph, and fst functions to terra SpatRaster raster objects, it's easy to generate any arbitrary cost function using data.table syntax.
1.2 The Energetics of Locomotion
In order to move across the terrestrial landscape, two principal types of work must be performed. The first is kinetic work, required to move an object at a particular speed, $K = \frac{1}{2}mv^2$ for an arbitrary particle:
```{=tex} \begin{equation} K = \frac{1}{4} m v^2 \frac{\ell_s}{L} \end{equation}
for human locomotion as derived by @Kuo2007: 654---where $m$ is an object's mass, $v$ it's speed, $\ell_s$ is the average step length and $L$ the average length of a leg; and for quadrupeds $K = m(0.685v + 0.072)$ as derived experimentally by @Heglund1982. LBMech also provides alternate energetic models including one proposed by @Pontzer2023 for any arbitrary animal $K = 8m^{-0.34} + 50\left(1 + \sin(2\theta - 74^\circ)\right)m^{-0.12}$, with $\theta$ in degrees), and one based on a simplified harmonic oscillator based on variation about an average velocity ($K = 2mv^2\gamma$, where $\gamma$ is the fractional difference between the maximum/minimum and average velocities). For the purposes of simplicity, the rest of the derivation will be using the appropriate parameters for humans according to Kuo but the steps are equivalent for the derivation of quadrupeds. The second type of work is work against gravity, required to lift an object a height $h$ ($U = mgh$), where $g$ is the magnitude of the acceleration due to gravity, and $h$ the height an object is raised. Importantly, this latter form of work cost is incurred only when moving uphill. Summing the two components results in the total amount of work performed to move across the landscape: ```{=tex} \begin{equation} \tag{1} W = \frac{1}{4} m v^2 \frac{\ell_s}{L} + \begin{cases} mgh & dz \geq 0 \\ 0 & dz < 0 \\\end{cases} \end{equation}
To obtain the expected velocity when moving across a landscape of variable topography, we can employ Tobler's Hiking function, which predicts that velocity decreases exponentially as the slope departs from an 'optimum' slope at which an traveler's velocity can be maximized:
```{=tex} \begin{equation} \tag{2} \frac{d\ell}{dt} = v_{\textrm{max}} e^{-k | \frac{dz}{d\ell} - s |} \end{equation}
where $\ell$ is horizontal movement, $z$ vertical movement, $v_{\textrm{max}}$ the maximum walking speed, $s$ the slope of maximum speed and $k$ how sensitive speed variation is to slope. Canonical applications of the Tobler function for humans employ $v_{\textrm{max}} = 1.5$ m/s, $k = 3.5$, and $s = - 0.05$ = $\tan(-2.83^\circ)$. However, the `getVelocity()` function in the `lbmech` package allows researchers to obtain the appropriate variables for their subject species of interest from high-resolution GPS tracking data.
Humans and animals, however, are far from ideal engines. Generally, for every calorie of work expended humans require five calories of dietary energy resulting in an efficiency factor of $\varepsilon = 0.2$. Dividing Eq. 1 by $\varepsilon$, and plugging in Tobler's Hiking Function (Eq. 2) for $v$ results in the total amount of energy used by a human to travel between two points on the landscape. This is our function for kinematic energy expenditure.
People also consume energy just to stay alive as defined by our base metabolic rate (BMR, represented by $\frac{dH}{dt}$). Assuming a relatively constant amount of daily physical activity and by extension a relatively constant basal metabolic rate, we can model the total amount of energy used by a human to travel between two points on the landscape and stay alive, by adding $H$ to our equation:
```{=tex}
\begin{equation}
\tag{3}
E = \frac{1}{4} m v^2 \frac{\ell_s}{L} + \frac{dH}{dt} dt + \begin{cases} mgh & dz \geq 0 \\ 0 & dz < 0 \\\end{cases}
\end{equation}
Both of the energetic equations derived above display anisotropic behavior due to dependence on $\theta = \arctan \left(\frac{dz}{dx}\right)$. However, only the potential energy term is linearly anisotropic whereas the kinetic energy term is dependent on slope through an exponential term to the second power. As such, the ArcGIS Distance toolkit is insufficient for true energetic analyses. This incapacity is largely because it---like other least-cost tools perform the calculations directly upon the localized value at a particular cell. Rather, since this is a problem of performing a line integral over a differential surface, an algorithm that performs its calculations based on the difference between individual cells is required.
The figure below plots each of the above cost functions (for time, kinetic energy, gravitational potential energy, mechanical/kinematic energy, and net metabolic energy) versus slope, assuming the canonical parameters for Tobler's Hiking Function and a human male. $k = 3.5$, $s = -0.05$ $m = 68$ kg, $v_\textrm{max} = 1.5$ m/s, $L = 0.75$ m, $\ell_s = 1.6$ m, $g = 9.81$ m/s$^2$, $\frac{dH}{dt} = 93$ J/s, and $\varepsilon = 0.2$; (a) Tobler's Hiking Function inverted for Pace $t$; (b) Kinetic Energy, $K$; (c) Work against Gravitational Potential Energy, $U$; (e) Kinematic/Mechanical Energy Performed, $W = K + U$; (f) Net Metabolic Energy Consumed, $E = \frac{W}{\varepsilon} + \frac{dH}{dt}t$.
library(data.table) library(ggplot2) funs <- data.table(Slope = seq(-0.25,0.25,length.out = 1000) )[, `(a) Pace (s/m)` := 1.5 * exp(3.5 * abs(Slope + 0.05)) ][, `(b) Kinteic Energy (J/m)` := 1 / 2 * 68 * `(a) Pace (s/m)`^-2 ][, `(c) Potential Energy (J/m)` := fifelse(Slope > 0, 68 * 9.81 * Slope, 0), ][, `(d) Mechanical Work (J/m)` := `(b) Kinteic Energy (J/m)` + `(c) Potential Energy (J/m)` ][, `(e) Net Metabolic Work (J/m)` := `(d) Mechanical Work (J/m)`/0.2 + 72 * `(a) Pace (s/m)`] f <- melt(funs,id.vars = 'Slope', variable.name = 'Cost') ggplot(f, aes(x = Slope, y = value, color = Cost)) + geom_line(show.legend=FALSE, linewidth = 1) + facet_wrap(Cost ~ ., scales = 'free', ncol = 2) + theme_bw() + theme(plot.title = element_text(hjust = 0.5, size = 13, face = 'bold'), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), strip.text.x = element_text(size = 12), strip.background =element_rect(fill=NA,linewidth=NA)) + scale_x_continuous(expand = c(0,0)) + scale_y_continuous(expand = c(0,0)) + ggtitle("Cost versus Slope using Canonical Parameters\nFor Various Costs") + ylab("Cost")
A least-cost path algorithm is going to preferentially select paths whose segments tend towards the lowest values on these plots. When minimizing the amount of time (Fig. 1a, and considering the canonical parameters), slopes closer towards the slope of fastest travel $s = - 0.05$ are preferred, with costs increasing exponentially as one moves away from that value. When minimizing the energy of motion (kinetic energy, Fig. 1b), the most extreme slopes will be preferred since movement is slowest and kinetic energy is minimized. When minimizing the work performed against gravitational potential energy (Fig. 1c), any negative slope is equally preferred, with positive slopes increasing in cost linearly. When minimizing the kinetic energy, steeper slopes are preferred downhill mildly over flat-ish slopes, with a notable avoidance of slopes tending around $s = -0.05$, with slopes linearly---but rapidly---increasing in positive slopes. In other words, downhill movement is primarily kinetic energy-limited, whereas uphill movement is mainly gravitational potential energy-limited. The situation is similar for net metabolic energy, but here additional penalties are incurred for overly-slow routes that \textit{may} be more mechanically efficient but take long enough that metabolic processes incur costs. Thus, the most efficient slopes are those that occur at the global minimum at $\frac{dz}{d\ell} = -0.0818 = \arctan(-4.7^\circ)$.
One final consideration that to date no other least-cost toolkit rigorously assesses is potential movement on and across water. Movement considerations across bodies of water are conceptually simpler than movement across land. The water's surface can be represented by a two-dimensional vector field with values at each point representing water's surface velocity $\vec{v}\textrm{water}(x,y)$, with still water having a velocity of $\vec{v}\textrm{water} = 0$. The velocity of someone traveling on water then is simply:
```{=tex} \begin{equation} \tag{4} \vec{v}(x,y) = \vec{v}\textrm{water} + \vec{v}\textrm{travel}, \end{equation}
where $\vec{v}_\textrm{travel}(x,y) = v_\textrm{travel} \hat{\ell}$ is the travel speed of the method of travel (such as swimming or rowing a boat) over still water, and $\hat{\ell}$ the intended direction of travel. If $\vec{v} < 0$, then travel in that direction is impossible. For the purposes of simplicity, we assume that the caloric cost is proportional to $v_\textrm{travel}^2$. ## 1.3 Estimating Velocity The GPS dataset consists of 15,370 usable GPX tracks obtained from <https://wikiloc.com>, and was first described by @Lera2017 who used it to identify seasonal variability in hiking trail activity and network organization. Given a directory with `.gpx` files, the tracks can be imported as a data.table using `lbmech::importGPX()`: ```r dir = tempdir() getData('baleares-gps', dir = dir) # Detect files with 'gpx' extension in working directory gpx <- list.files(dir, recursive = TRUE, full.names = TRUE, pattern = ".gpx$") name <- stringr::str_extract(gpx,'(?<=/)[A-Za-z0-9]+(?=.gpx)') # Import with lbmech::importGPX(), this will take a while gpx <- lapply(gpx,importGPX) names(gpx) <- name # Join into a single object gpx <- rbindlist(gpx,idcol = TRUE)
Since GPS tracks can be sampled at a per-length rate smaller than the size of even the smallest raster pixels we might employ, GPS tracks are downsampled to a comparable rate---about a pixel length per GPS sample. For an expected pixel size of 50 m and an approximate maximum speed of 1.5 m/s, t_step = 50/1.5. We can use the lbmech::downsampleXYZ() function for this, applying it on all GPS observations with an associated time stamp (filtering for !is.na(t)):
gpx <- downsampleXYZ(gpx[!is.na(t)], t_step = 50/1.5, t = 't', x = 'long', y = 'lat', z = 'z', ID = '.id')
library(lbmech) dir <- "D:/lbmech/Altar" rd <- "D:/lbmech/Altar/Baleares" gpx <- fst::read_fst(paste0(rd,"/gpx_filt.fst"),as.data.table=TRUE)
The xyz data is then passed to lbmech::getVelocity(), which for each sequence of GPS points (1) calculates the changes in elevation $\frac{dz}{d\ell}$ and planimetric speed $\frac{d\ell}{dt}$, and (2) performs a nonlinear quantile regression
to get a function of the form Tobler, which returns a list with the nonlinear model, the model parameters, and the transformed data.
velocity_gps <- getVelocity(data = gpx[y < 90], # Filter for y, degs = TRUE, # Lat/Long; geodesic correction tau_vmax = 0.95, # Quantile for v_max tau_nlrq = 0.50, # Quantile for nlrq v_lim = 3) # Filter for dl/dt print(velocity_gps[1:(length(velocity_gps)-1)])
Finally, the lbmech::plotVelocity() function can be used to plot the log-transformed probability of a given speed at a given slope versus the regressed function:
plotVelocity(velocity_gps)
1.4 Energetic Cost-Distance Analysis
1.4.1 Procuring Data
Given the large amounts of data downloaded, processed, and re-used at different parts of the workflow, it is highly recommended to define a consistent directory for the workflow particularly during model-building stages.
# Define a working directory rd <- "Baleares" if (!dir.exists(rd)){ dir.create(rd) }
lbmech can automatically download terrain data, but we also want to consider ocean movement. Here, we'll import the Mediterranean Sea Physics Analysis and Forecast (CMEMS MED-Currents, EAS6 system) [@Clementi2021] for ocean currents around two of the principal Balearic Islands---Mallorca and Menorca---on 13 June 2022 and define our region of maximum movement (our world) around that raster's extent:
# Projection will be UTM 31N WGS1984 proj <- "EPSG:32631" # Import an ocean water surface velocity dataset currents <- project(getData('baleares-currents'), proj) # Define region of maximum possible movement region <- as.polygons(ext(currents),proj)
library(lbmech) library(data.table) library(ggplot2) library(tidyterra) library(ggspatial) library(pandoc) library(png) library(grid) library(magick) library(cowplot) library(ggquiver) setwd("D:/lbmech/Altar/vignettes") rd <- "D:/lbmech/Altar/Baleares" # Projection will be UTM 31N WGS1984 proj <- "EPSG:32631" # Import Ocean current data currents <- rast(paste0(rd,"/Ocean_Currents.nc")) currents <- project(currents,proj) currentU_vals <- values(currents[[1]]) currentV_vals <- values(currents[[2]]) i <- !is.na(currentU_vals) coords <- as.data.table(xyFromCell(currents, 1:ncell(currents))) coords <- coords[as.logical(i)] currentU_vals <- currentU_vals[as.logical(i)] currentV_vals <- currentV_vals[as.logical(i)] tpsU <- fields::Tps(coords, currentU_vals) tpsV <- fields::Tps(coords, currentV_vals) currents[[1]] <- terra::interpolate(rast(currents[[1]]),tpsU) currents[[2]] <- terra::interpolate(rast(currents[[2]]),tpsV) # Define region of maximum possible movement region <- as.polygons(ext(currents),proj) region$ID <- 1 # Define the raster properties z_fix <- fix_z(proj = proj, res = 50) # Make a grid grid <- makeGrid(dem = region, nx = 15, ny = 15, sources = TRUE, zoom = 11, overlap = 0.05) # Import list of locations on Mallorca, Menorca and Cabrera pobles <- getData('baleares-places') pobles <- project(pobles,proj) muns <- paste0(rd,"/Borders.shp") if (!file.exists(muns)){ download.file("https://stacks.stanford.edu/file/druid:sp807yd2954/data.zip?download=true", paste0(tempdir(),"/SpainBorders.zip"),mode='wb') unzip(paste0(tempdir()),overwrite=TRUE,junkpaths = TRUE) muns <- list.files(paste0(tempdir(),"/SpainBorders"),pattern=".shp$",full.names = TRUE) muns <- vect(muns) muns <- project(muns,proj) muns <- crop(muns, currents) writeVector(muns,outline) } muns <- vect(muns) plotreg <- as.polygons(ext(buffer(muns,12500)),crs = crs(muns)) insetExt <- sf::st_read("./Outputs/Baleares_Region.kml") insetExt <- project(vect(insetExt),proj) islands <- aggregate(muns) countries <- list.files("D:/Census/Global_Divisions/",pattern="_0.shp",full.names=TRUE) countries <- lapply(countries,vect) countries <- do.call(rbind,countries) countries <- project(countries,proj) countries <- crop(countries,ext(insetExt)) states <- list.files("D:/Census/Global_Divisions/",pattern="_1.shp",full.names=TRUE) states <- lapply(states,vect) states <- do.call(rbind,states) states <- project(states,proj) states <- crop(states,ext(insetExt)) dem <- rast(elevatr::get_elev_raster(raster::raster(project(currents,"+proj=longlat")), z = 9)) dem[dem < 0] <- NA dem <- project(dem,proj) dem <- crop(dem,plotreg) slope <- terrain(dem,v='slope',unit='radians') aspect <- terrain(dem,v='aspect',unit='radians') # Calculate a hillshade raster hill <- shade(slope,aspect) names(hill) <- "shades" # Get a pallete of greys for hillshading pal_greys <- hcl.colors(1000, "Grays") # Round the hillshade values to the nearest factor hill <- (hill - minmax(hill)['min',]) / (minmax(hill)['max',] - minmax(hill)['min',]) hill <- ceiling(hill * (length(pal_greys)-1)) + 1 # Get colors for the elevation DEM subplot grad_hypso <- hypso.colors2(10, "dem_poster") names(dem) <- "Elevation" # Transform currents into data.table to turn into vector arrows field <- project(currents, fix_z(proj,res=8000),align=TRUE) field <- crop(field,plotreg) field <- as.data.table(rastToTable(field)) names(field)[1:2] <- c("u","v") field <- na.omit(field) dem[dem == 0] <- NA hill[is.na(dem)] <- NA topo <- ggplot() + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, linewidth = 0.75, col = 'lightblue', inherit.aes=FALSE) + geom_spatraster(data= dem, na.value = NULL) + scale_fill_gradientn(" \nElevació\n(m)", colours = hypso.colors(15,'dem_poster'), limits = c(0,unlist(global(dem,'max',na.rm=TRUE))), na.value = "transparent") + theme_bw() + theme(legend.position = "top", legend.text = element_text(angle = 45, size = 15), legend.title = element_text(size=20), strip.text = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_rect(fill = '#005EB8') ) + guides(fill = guide_colorbar(barwidth = 60, barheight = 1.5, alpha = 0.7, label.position = "top", override.aes = list(alpha = 0.7))) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf, na.value = NULL) + geom_spatvector(data = muns, fill=NA, size=1,col = '#EEEEEE') + geom_spatvector(data = islands, linewidth = 2.5, col = '#454545',fill=NA) + geom_spatvector(data = pobles, color = 'red', size = 5) + ggrepel::geom_label_repel(data = pobles, mapping = aes(x = geom(pobles)[,'x'], y = geom(pobles)[,'y'], label = pobles$Location), size = 7) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) + annotation_scale(location = 'br', height = unit(1,'cm'),text_cex = 3) + annotation_north_arrow(location = 'tl', style = north_arrow_nautical, height = unit(4.5,'cm'), width = unit(4.5,'cm')) political <- ggplot(states) + geom_spatvector(fill = '#396b27', linewidth = 0.5, col = '#D0D0D0') + geom_spatvector(data = countries, linewidth = 1, col = 'black', fill= NA) + geom_spatvector(data = plotreg, linewidth = 2, col = 'red', fill = NA) + theme_bw()+ theme(axis.text.y=element_blank(), #remove y axis labels axis.ticks.y=element_blank(), #remove y axis ticks axis.text.x=element_blank(), #remove x axis labels axis.ticks.x=element_blank(), panel.background = element_rect(fill='#bce6f9'), panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + scale_x_continuous(expand=c(0,0)) + scale_y_continuous(expand = c(0,0)) header <- ggdraw() + draw_image(paste0(rd,"/RegionMap_Header.png")) side <- plot_grid(header,political,nrow=2, rel_heights = c(0.66,0.33)) gg <- plot_grid(topo,side,ncol=2,rel_widths = c(0.66,0.33)) gg # pdf(file = "D:/lbmech/Figs/Map.pdf", # height = 9.5, width =20) # gg # dev.off()
1.4.2 Defining a 'world'
As with any raster operation, a consistent projection and grid needs to be defined to snap/fix all sources. However, unlike most raster packages lbmech does not require foreknowledge of the spatial extents; only the resolution and offsets. A raster can be provided, or one can be made with the lbmech::fix_z() function:
# Define the raster properties z_fix <- fix_z(proj = proj, res = 50)
Although a world can be defined directly based on a digital elevation model, generally it is easier to define a polygon coincident with the coverage of a digital elevation model and an attribute pointing to a URL/filepath with the source. The lbmech::makeGrid() can make such a polygon for a raster or filepath input, while using a polygon will induce future functions to download SRTM data as-needed. Similar polygons are frequently distributed by state GIS agencies (e.g. Pennsylvania; Massachusetts):
# Make a grid grid <- makeGrid(dem = region, # Input; here polygon for SRTM nx = 15, ny = 15, # Cols/Rows to divide the polygon sources = TRUE, # Just crop/divide, or point to source? zoom = 11, # Zoom level for SRTM data overlap = 0.05) # Fraction overlap between adjacent tiles
A world is formally defined as a directory containing all possible constraints on movement, fixed to a particular spatial projection. Given the very large areas and relatively-small spatial resolutions often employed, the data can rapidly become too large to compute all at once in the memory. The defineWorld() function divides the potential world of movement into individual, slightly-overlapping sectors that are calculated and read-in only as needed. It itself generally is not used to perform calculations.
# Resample currents so they're at a more realistic resolution currents_high <- project(currents, disagg(currents, ceiling(res(currents)/50)), method = 'cubicspline') # Define world of motion within a workspace defineWorld(source = grid, # Elevation source, like makeGrid() output grid = grid, # How to divide up the world directions = 8, # How adjacency between cells are defined neighbor_distance = 10, # Overlap between tiles in addition to ^ cut_slope = 0.5, # Maximum traversible slope water = currents_high, # Water velocity source priority = 'land', # If data for land or water, who wins z_min = 0, # Minimum elevation, below NA in land z_fix = z_fix, # Grid with defined projection, resolution dist = 'karney', # Geodesic correction method dir = rd, # Working directory overwrite=FALSE)
1.4.3 Estmiating Costs
Once a world has been defined, a cost function can be applied using the lbmech::calculateCosts() function. The lbmech::energyCosts() function is provided to perform least-time, least-work, and least-energy analyses as in the above sections of this manuscript. As with lbmech::defineWorld(), it itself is not generally used to perform calculations---just to define the values.
calculateCosts(costFUN = energyCosts, dir = rd, # Working directory with world method = "kuo", # Method to calculate work water = TRUE, # Consider water, or only land? v_max = velocity_gps$vmax, # Max walking speed (m/s) k = velocity_gps$k, # Slope sensitivity (dimensionless) s = velocity_gps$s, # Slope of fastest motion (dimlss.) row_speed = 1.8, # Speed over water (m/s) row_work = 581, # Work over water (J/s) m = 68, # Mass (kg) BMR = 72, # Basal metabolic rate (J/s) l_s = 1.8, # Stride length (m) L = 0.8) # Leg length
For a given (set of) cost function(s) and origin(s)/destination(s), getCosts() will identify the needed sectors, make sure the data has been downloaded (and do so if not), perform the calculations, and save the cell-wise transition cost tables to the working directory. That way, if they are needed for another calculation in the future they won't need to be fetched and prepossessed again.
# Import locations on Mallorca, Menorca, and Cabrera pobles <- project(getData("baleares-places"), region) reg <- as.polygons(ext(buffer(pobles, 10000)),crs = crs(proj)) costs <- getCosts(region = reg, # Area of maximum movement from = pobles, # Origins to = NULL, # Destinations costname = 'energyCosts', # Name of costFUN above id = 'Location', # Column with origin names dir = rd, # Directory with world output = c('object','file'), # Save files destination = 'all') # Distance to all points in region
1.4.4 Paths and Corridors
Likewise, least-cost paths can be computed given a set of nodes:
# Find all least-cost paths between Palma de Mallorca and Maó paths <- getPaths(region = reg, nodes = pobles, costname = 'energyCosts', order = c("Palma de Mallorca","Maó"), id = 'Location', dir = rd)
reg <- as.polygons(ext(buffer(pobles, 10000)),crs = crs(proj)) # Import list of locations on Mallorca, Menorca and Cabrera pobles <- fread("D:/lbmech/Altar/vignettes/localitats.csv") # Convert to SpatVector pobles <- vect(pobles, geom=c('long','lat'), crs = '+proj=longlat') pobles <- project(pobles,proj) # Uncomment if you've never run this section before before, since the figure # differs from the example given in the text for clarity ## Get combinations of all paths to and from Palma de Mallorca # path <- pobles$Location[pobles$Location != 'Palma de Mallorca'] # path <- c(rbind(rep("Palma de Mallorca",length(path)+1), # path)) # path <- path[1:(length(path)-1)] ## Calculate the paths # palma_paths <- getPaths(region = reg, # nodes = pobles, # costname = 'energyCosts', # order = path, # id = 'Location', # dir = rd) # Import the least-cost paths palma_paths <- vect("D:/lbmech/altar/baleares/world/costrasters/PalmaPaths.shp") # Convert currents to a (smaller) data.table field <- project(currents, fix_z(proj,res=8000),align=TRUE) field <- crop(field,reg) field <- as.data.table(rastToTable(field)) names(field)[1:2] <- c("u","v") field <- na.omit(field) # Get even coarser DEM for hillshading dem <- rast(elevatr::get_elev_raster(raster::raster(project(currents,"+proj=longlat")), z = 8)) dem[dem < 0] <- 0 dem <- project(dem,proj) dem <- crop(dem,reg) # Calculate aspect, slope slope <- terrain(dem,v='slope',unit='radians') aspect <- terrain(dem,v='aspect',unit='radians') # Calculate a hillshade raster hill <- shade(slope,aspect) names(hill) <- "shades" # Get a pallette of greys for hillshading pal_greys <- hcl.colors(1000, "Grays") # Round the hillshade values to the nearest factor hill <- (hill - minmax(hill)['min',]) / (minmax(hill)['max',] - minmax(hill)['min',]) hill <- ceiling(hill * (length(pal_greys)-1)) + 1 # Import the cost rsters costs <- rast(paste0(rd,"/World/CostRasters/Pobles.tif")) # Keep only the ones involving Palma palma <- costs[[stringr::str_detect(names(costs),"Palma de Mallorca")]] # The time ones palma_dt <- palma[[stringr::str_detect(names(palma),"dt_")]] # SpatVector Collection with to, from rasters for time palma_dt_paths <- svc(palma_paths[stringr::str_detect(palma_paths$segment,"Mallorca_to_") & palma_paths$cost=='dt',], palma_paths[stringr::str_detect(palma_paths$segment,"_to_Palma") & palma_paths$cost=='dt',]) # ggplot call for the from raster with time, involving the islands, hillshade # cost raster, cost paths, and ocean currents ggtime_from <- ggplot(islands) + geom_spatraster(data = palma_dt[[1]]/3600) + scale_fill_viridis_c(expression(paste(atop("Hores", italic("Hours")," (hr)"))), option='magma',direction = -1, limits = rowMeans(minmax(palma_dt))/3600) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = palma_dt[[1]],binwidth = 3600) + geom_spatvector(data = palma_dt_paths[[1]], linewidth = 1.4, col = 'yellow') + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "right", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), axis.title.y = element_blank(), strip.text.x = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + xlab(expression(paste("Temps des de Palma — ",italic("Time from Palma")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) + ylab("") # Now likewise for time to ggtime_to <- ggplot(islands) + geom_spatraster(data = palma_dt[[2]]/3600) + scale_fill_viridis_c(expression(paste(atop("Hores", italic("Hours")," (hr)"))), option='magma',direction = -1, limits = rowMeans(minmax(palma_dt))/3600) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = palma_dt[[2]],binwidth = 3600) + geom_spatvector(data = palma_dt_paths[[2]], linewidth = 1.4, col = 'yellow') + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "none", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), ) + xlab(expression(paste("Temps cap a Palma — ",italic("Time to Palma")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) + ylab("") # And likewise for work palma_dW <- palma[[stringr::str_detect(names(palma),"dW_")]] palma_dW_paths <- svc(palma_paths[stringr::str_detect(palma_paths$segment,"Mallorca_to_") & palma_paths$cost=='dW_l',], palma_paths[stringr::str_detect(palma_paths$segment,"_to_Palma") & palma_paths$cost=='dW_l',]) ggwork_from <- ggplot(islands) + geom_spatraster(data = palma_dW[[1]]/1000/4.184) + scale_fill_viridis_c(expression("(kcal)"), option='magma',direction = -1, limits = rowMeans(minmax(palma_dW))/1000/4.184) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = palma_dW[[1]],binwidth =500000*4.184) + geom_spatvector(data = palma_dW_paths[[1]], linewidth = 1.4, col = 'yellow') + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "right", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), axis.title.y = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) + ylab("") + xlab(expression(paste("Treball des de Palma — ",italic("Work from Palma")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) ggwork_to <- ggplot(islands) + geom_spatraster(data = palma_dW[[2]]/1000/4.184) + scale_fill_viridis_c(expression("(kcal)"), option='magma',direction = -1, limits = rowMeans(minmax(palma_dW))/1000/4.184) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = palma_dW[[2]],binwidth =500000*4.184) + geom_spatvector(data = palma_dW_paths[[2]], linewidth = 1.4, col = 'yellow') + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "none", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), axis.title.y = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) + ylab("") + xlab(expression(paste("Treball cap a Palma — ",italic("Work to Palma")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) # And finally energy palma_dE <- palma[[stringr::str_detect(names(palma),"dE_")]] palma_dE_paths <- svc(palma_paths[stringr::str_detect(palma_paths$segment,"Mallorca_to_") & palma_paths$cost=='dE_l',], palma_paths[stringr::str_detect(palma_paths$segment,"_to_Palma") & palma_paths$cost=='dE_l',]) ggenergy_from <- ggplot(islands) + geom_spatraster(data = palma_dE[[1]]/1000/4.184) + scale_fill_viridis_c(expression("(kcal)"), option='magma',direction = -1, limits = rowMeans(minmax(palma_dE))/1000/4.184) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = palma_dE[[1]],binwidth =1250000*4.184) + geom_spatvector(data = palma_dE_paths[[1]], linewidth = 1.4, col = 'yellow') + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "right", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), axis.title.y = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) + ylab("") + xlab(expression(paste("Energia des de Palma — ",italic("Energy from Palma")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) ggenergy_to <- ggplot(islands) + geom_spatraster(data = palma_dE[[2]]/1000/4.184) + scale_fill_viridis_c(expression("(kcal)"), option='magma',direction = -1, limits = rowMeans(minmax(palma_dE))/1000/4.184) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = palma_dE[[2]],binwidth =1250000*4.184) + geom_spatvector(data = palma_dE_paths[[2]], linewidth = 1.4, col = 'yellow') + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "none", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), axis.title.y = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) + ylab("") + xlab(expression(paste("Energia cap a Palma — ",italic("Energy to Palma")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0))
# Final figure is the above 6 with ggarrange. Can't do it by facets since the scales are all different. egg::ggarrange(ggtime_to,ggtime_from, ggwork_to,ggwork_from, ggenergy_to,ggenergy_from, nrow=3)
Corridors describe the minimum cost required to take a detour to a given location on a path between two or more given points relative to the least cost path between those points. These can be calculated from the output cost rasters using the lbmech::makeCorridor() function.
corr <- makeCorridor(rasters = costs, order = rev(c("Castell de Cabrera","Sa Calobra","Maó")))
ggtime_corr <- ggplot(islands) + geom_spatraster(data = corr[[1]]/3600) + scale_fill_viridis_c(expression("(hrs)"), option='magma',direction = -1, limits = rowMeans(minmax(corr[[1]]))/3600) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = corr[[1]]/3600,binwidth = 0.5) + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "none", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), axis.title.y = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) + ylab("") + xlab(expression(paste("Temps de desviació — ",italic("Time of detour")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) ggwork_corr <- ggplot(islands) + geom_spatraster(data = corr[[2]]/4.184/1000) + scale_fill_viridis_c(expression("(kcal)"), option='magma',direction = -1, limits = rowMeans(minmax(corr[[2]]/4.184/1000))) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = corr[[2]]/4.184/1000,binwidth = 500) + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "none", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), axis.title.y = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) + ylab("") + xlab(expression(paste("Treball de desviació — ",italic("Work of detour")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) ggenergy_corr <- ggplot(islands) + geom_spatraster(data = corr[[3]]/4.184/1000) + scale_fill_viridis_c(expression("(kcal)"), option='magma',direction = -1, limits = rowMeans(minmax(corr[[3]]/4.184/1000))) + geom_spatraster(data = hill, fill = pal_greys[as.matrix(hill)], alpha=0.4, maxcell = Inf) + geom_quiver(field, mapping = aes(x=x,y=y,u=u,v=v),center=TRUE, col = 'white', inherit.aes=FALSE) + geom_spatvector( linewidth = 1.25, col = '#000000',fill=NA) + geom_spatraster_contour(data = corr[[3]]/4.184/1000,binwidth = 1000) + geom_spatvector(data = pobles, size = 3, color = 'red') + theme_bw() + theme(legend.position = "none", legend.text = element_text(angle = 45, size = 10), legend.title = element_text(size=12, hjust=0), strip.text.x = element_blank(), axis.title.y = element_blank(), strip.background =element_rect(fill=NA,linewidth=NA), panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) + ylab("") + xlab(expression(paste("Energia de desviació — ",italic("Energy of detour")))) + scale_x_continuous(position = 'top', expand=c(0,0)) + scale_y_continuous(position = 'right', expand=c(0,0)) cowplot::plot_grid(ggtime_corr, ggwork_corr, ggenergy_corr, ncol = 1)
1.5 Influence
Different measures of distance in different topologies are often taken to be measures of influence. In other words, nearer things are more likely to influence an observation than farther things. LBMech provides tools to represent these relationships as various types of neighbor weights (for influence between observations) and Voronoi polygons (for the relative influence of observations everywhere else in space).
1.5.1 Neighbor Weights
Distance matrices can be obtained from the getCosts() function by setting the destination parameter to 'pairs':
dists <- getCosts(region = reg, # Area of maximum movement from = pobles, # Origins costname = 'energyCosts', # Name of costFUN above costFUN = energyCosts, id = 'Location', # Column with origin names dir = rd, # Directory with world destination = 'pairs') # Distance to all points in region
We can try and visually compare the different topologies using the normMDS() function, which performs multidimensional scaling on each distance matrix and rescales/orients them so they are in units of---and directionally relative to---a common path:
uv <- normMDS(dists, # List of distance matrices xy = pobles, # Cartesian coordinates (optional) id = 'Location', # Unique ID column ax = c('Palma de Mallorca','Inca'), # Reference path chir = 'Sa Calobra') # Third location to preserve chirality
colors <- c(ref = 'black', dt = 'blue', dW_l = 'darkgreen', dE_l = 'darkred') uvplot <- rbindlist(uv, idcol = TRUE) ax <- c('Palma de Mallorca','Inca') axplot <- uvplot[Location %in% ax & .id == 'ref'] uvplot[Location %in% ax & .id != 'ref', Location := ''] ggplot(data = uvplot, mapping = aes(x = u, y = v, label = Location, color = as.factor(.id))) + scale_color_manual('Distance', values = colors) + geom_point() + ggrepel::geom_text_repel(size = 2.5) + geom_point(data = axplot, mapping = aes(x = u, y = v), color = 'black', fill = 'yellow', shape = 23, size = 2.5) + ggtitle('Different implied 2D topologies based on different costs', subtitle = 'MDS output aligned and scaled such that\nthe paths between Palma and Inca are equivalent') + theme(plot.title = element_text(hjust = 0.5), plot.subtitle = element_text(hjust = 0.5)) + coord_equal()
lbmech provides makeWeights which can convert any distance matrix (geographic, or cost distance-based) into various types of spatial weights matrices. At present, it supports:
- Cutoff distances and distance-decay formulas for determining group membership
- Any arbitrary distance-decay formula
- Rank-distance in lieu of raw distance
- Tools to deal with singularities arising from distance-decay
- Row standardization (where all rows sum to one)
- Fuzzy row standardization (where entries are a proportion of the largest value of a row)
- Binary weights matrices depending on group membership if provided only a vector of group IDs
weights <- list() # The three nearest neighbors time-wise are definitely neighbors to an observation # Everyone else definitely isn't. weights[[1]] <- makeWeights(dists$dt, bw = 3) # All observations within 20 MJ of work are neighbors, with row-standardized weights # based on inverse distance, assuming self-weights are equivalent to 20 MJ weights[[2]] <- makeWeights(dists$dW_l, bw = 2e7, offset = 2e7, mode = 'fixed', weighting = 'distance', row.stand = TRUE) # All observations might be neighbors. # The weights represent the probability that they're neighbors. # The probability is based on the inverse distance square of the rank of the energy distance weights[[3]] <- makeWeights(dists$dE_l, weighting = 'rank', FUN = function(x) 1/x^2, row.stand = 'fuzzy') # Neighbors are those locations on the same island weights[[4]] <- makeWeights(pobles$Island, ID = pobles$Location)
library(ggraph) set.seed(43884) w <- as.matrix(weights[[1]]) w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE) w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w))) ggw1 <- ggraph(w, algorithm = 'graphopt', 'igraph') + geom_point(aes(x = x, y = y)) + geom_edge_link0() + ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3) w <- as.matrix(weights[[2]]) w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE) w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w))) ggw2 <- ggraph(w, algorithm = 'kk', 'igraph') + geom_point(aes(x = x, y = y)) + geom_edge_link0() + ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3) w <- as.matrix(weights[[3]]) w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE) w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w))) ggw3 <- ggraph(w, algorithm = 'gem', 'igraph') + geom_point(aes(x = x, y = y)) + geom_edge_link0(edge_alpha = igraph::as_data_frame(w)$weight) + ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3) w <- as.matrix(weights[[4]]) w <- igraph::graph_from_adjacency_matrix(w, mode = 'directed', weighted = TRUE) w <- igraph::set_vertex_attr(w, 'name',value = names(igraph::V(w))) ggw4 <- ggraph(w, algorithm = 'gem', 'igraph') + geom_point(aes(x = x, y = y)) + geom_edge_link0() + ggrepel::geom_label_repel(aes(x = x, y = y, label = name), size = 3) cowplot::plot_grid(ggw1,ggw2,ggw3,ggw4,ncol=2)
Note that in the top-left subplot, only edges that exist in both directions are plotted.
1.5.2 Catchments and Voronoi Polygons
Voronoi polygons allocate all land to the nearest observation. These are most commonly employed using planar geographic distance as the metric and generally assume equal influence between all observations, such as with the voronoi() function from terra:
v <- voronoi(pobles)
Sometimes the assumption of equal influence does not hold. There exist a number of ways to differentially weight the influence any one point has versus another. Multiplicatively-weighted diagrams are particularly useful when a point's ability to influence space versus a different point is constant. For example, if we wanted to determine the catchment of two towns based on who can get to each point the fastest---with Town 1 possessing a car capable of traveling twice as fast as Town 2---we might consider a Multiplicatively Weighed Voronoi diagram giving Town 1 a weight of $w_1 = 2$ and Town 2 a weight of $w_2 = 1$. In short, Multiplicatively Weighted Voronoi diagrams transform distance to each point such that:
{=tex}
D = w_i * \ell_i
where $D$ is the modified distance, $w_i$ the weight of point $i$, and $\ell$ the raw distance to point $i$. lbmech provides a fully vectorized function to calculate such catchments for planar, spherical, and Earth-geoid topologies through the lbmech::mwVoronoi function:
While lbmech::mwVoronoi is marginally less efficient thant terra::voronoi for equal weighting schemes, it is capable of employing multiplicative weights and performing geodesic corrections.
Likewise, the Voronoi concept can be extended to any arbitrary topology---such as the time-, work- and energy-based distances obtained above. These are often known as 'catchments' in the literature. As long as we have previously calculated the cost rasters for the desired locations and costs, we can obtain the catchments with the lbmech::makeCatchment function:
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.