In kbvernon/hiker: Model Hiking Paths Over Costly Terrains

knitr::opts_chunk$set(
  collapse = TRUE,
  strip.white = TRUE,
  comment = "#>",
  out.width = "100%",
  warning = FALSE,
  message = FALSE,
  fig.retina = 3,
  fig.asp = 0.7
)

The hiker package serves as an aid to the modeling of specifically human movement across a landscape. Much of its code was inspired by Jacob van Etten's (2017) gdistance package, though it works with modern spatial tools, specifically sf and terra.

Functions

hiker exposes two sets of functions to the user.

1. Functions to build and modify cost surfaces.
2. hf_terrain() generates a cost-surface,
3. hf_barrier() adds an impassable barrier to travel,

4. hf_channel() adds a channel of reduced cost to travel.

5. Functions to estimate costs and identify paths.

6. hf_hike() estimates the shortest or least-cost path between points,
7. hf_appraise() measures the cost of hiking the shortest path between points, and
8. hf_survey() calculates the cost of travel from a point to anywhere else.

In addition to these, there is also hf_rasterize(), which turns a cost terrain into a SpatRaster. This is mostly for visualizing with terra::plot().

Note the 'hf' prefix. Short for "hiking function," it helps avoid potential namespace collisions.

The motivation

For demonstration purposes, imagine going for a hike around the Red Butte Canyon Research Natural Area (RNA) in the Wasatch Foothills above the University of Utah. Maybe you want to know how long it would take to hike from one ridge to another or from the mouth to the crest or even just from any random point to any other location. This is the kind of problem hiker is intended to solve.

library(hiker)
library(sf)
library(terra)
library(viridis)

Building a cost-surface

To estimate travel time, or to find the least-cost path, the first step is to generate a cost surface. This is the terrain being hiked in R.

fn <- system.file("extdata/red_butte_dem.tif", package = "hiker")
red_butte_dem <- rast(fn)

terrain <- hf_terrain(red_butte_dem)

The result of hf_terrain() is a terrain object. This is just a list with information about its spatial extent (bb8 in the list) and coordinate reference system (specifically, its epsg), along with a sparse adjacency matrix representing the inverse of travel cost, what Jacob van Ettern refers to as conductance. A simple print method is available for terrains, which displays this information.

terrain

par(mfrow = c(1, 2))

plot(red_butte_dem, main = "Elevation (meters)", col = viridis(50))
plot(terrain, main = "Travel cost (seconds)", col = viridis(50))

By default, hf_terrain() implements the Campbell hiking function (CHF), but you can also request Tobler's hiking function (THF) by passing the character string 'tobler' to the hf argument,

hf_terrain(red_butte_dem, hf = "tobler")

Additional arguments can be passed down to the specific user-selected hiking function. For instance, CHF will estimate velocity relative to the deciles of the sampled hikers used in the original analysis, with lower deciles representing slower hikers, upper deciles faster hikers. In hiker, CHF defaults to the median.

hf_terrain(red_butte_dem, hf = "campbell", decile = 90)

Modifying a cost surface

The RNA has an old reservoir built by the US military in 1930 to supply Fort Douglass (now part of the University of Utah). Obviously, we wouldn't expect a hiker to defy the laws of physics by walking on water, so we need to add that information to the terrain. This is done with hf_barrier().

fn <- system.file("extdata/red_butte_reservoir.geojson", package = "hiker")
red_butte_reservoir <- read_sf(fn)

terrain2 <- hf_barrier(terrain, red_butte_reservoir)

par(mfrow = c(1, 2))

plot(terrain,  main = "No barrier", col = viridis(50))
plot(terrain2, main = "Barrier", col = viridis(50))

We can also channelize a cost surface by making it easier to pass through an area.

fn <- system.file("extdata/red_butte_road.geojson", package = "hiker")
red_butte_road <- read_sf(fn)

terrain3 <- hf_channel(terrain, 
                       channel = st_buffer(red_butte_road, dist = 100),
                       .m = 0.1)

par(mfrow = c(1, 2))

plot(terrain,  main = "No channel", col = viridis(50))
plot(terrain3, main = "Channel", col = viridis(50))

Path functions

Now, let's suppose we have these start and end points.

from <- st_sf(
  geometry = st_sfc(
    st_point(c(432000, 4514000)),
    st_point(c(434000, 4518000))
  ),
  crs = 26912
)

to <- st_sf(
  geometry = st_sfc(
    st_point(c(431000, 4515000)),
    st_point(c(434000, 4515500)),
    st_point(c(436500, 4518500))
  ),
  crs = 26912
)

Travel cost

In the simplest case, we just want to know the cost of traveling the shortest path from each start to each end point. For this, hiker provides the hf_appraise() function, as in let's appraise the cost of taking the short path.

hf_appraise(terrain, from, to)

As you see, this returns a data.frame in long format, with columns for the from-points, to-points, and cost of travel between them.

Travel path

What about those shortest paths, though? Where are they? For this, we can use the hf_hike() function, as in let's take a hike along the short paths.

short_paths <- hf_hike(terrain, from, to)

plot(red_butte_dem, col = viridis(50))

plot(st_geometry(short_paths),
     col = "white",
     lwd = 1.5,
     add = TRUE)

plot(st_geometry(from), 
     pch = 21, 
     cex = 1.3, 
     bg = "red2", 
     col = "white", 
     add = TRUE)

plot(st_geometry(to), 
     pch = 21, 
     cex = 1.3, 
     bg = "dodgerblue3", 
     col = "white", 
     add = TRUE)

Survey cost

Finally, we can ask how long it would take to get from any one point to all other locations.

accumulated_cost <- hf_survey(terrain, from)

plot(accumulated_cost, col = viridis(50))

plot(st_geometry(from), 
     pch = 21, 
     cex = 1.3, 
     bg = "red2", 
     col = "white", 
     add = TRUE)

Piping

The two functions designed to modify a cost surface, hf_barrier() and hf_channel() are kinda-sorta pipeable. The first argument to each is a data object, in this case a terrain, and they both return a modified version of the same.

library(dplyr)

red_butte_dem |> 
  hf_terrain() |> 
  hf_barrier(red_butte_reservoir) |> 
  hf_channel(red_butte_road, .m = 0.6) |> 
  hf_hike(from, to)

Working with lines and polygons

Unfortunately, it's not obvious - not obvious to me, anyway - how best to handle lines and polygons, so the path functions in hiker work exclusively with POINT and MULTIPOINT geometries. That said, there are a few, shall we say, less than desireable methods for working with other geometries.

Center

from <- st_centroid(red_butte_reservoir)

accumulated_cost <-
  red_butte_dem |>
  hf_terrain() |>
  hf_survey(from)

bb8 <- 
  red_butte_reservoir |> 
  st_buffer(dist = 100) |> 
  st_bbox() |> 
  st_as_sfc(crs = st_crs(red_butte_reservoir))

plot(crop(accumulated_cost, bb8), col = viridis(50))
plot(st_geometry(from), 
     pch = 21, 
     cex = 1.3, 
     bg = "red2", 
     col = "white", 
     add = TRUE)

Cast

You can cast a LINESTRING or POLYGON to POINT, though beware that this will get out of hand pretty quick if you have a large shape or one with densely packed vertices.

from <- st_cast(red_butte_reservoir, "POINT")

accumulated_cost <-
  red_butte_dem |>
  hf_terrain() |>
  hf_survey(from)

plot(crop(accumulated_cost, bb8), col = viridis(50))
plot(st_geometry(from), 
     pch = 21, 
     cex = 1.3, 
     bg = "red2", 
     col = "white", 
     add = TRUE)

Sample

from <-
  red_butte_reservoir |>
  st_cast("POINT") |>
  dplyr::sample_n(round(n()/10))

accumulated_cost <-
  red_butte_dem |>
  hf_terrain() |>
  hf_survey(from)

plot(crop(accumulated_cost, bb8), col = viridis(50))
plot(st_geometry(from), 
     pch = 21, 
     cex = 1.3, 
     bg = "red2", 
     col = "white", 
     add = TRUE)

kbvernon/hiker documentation built on Dec. 9, 2022, 11:16 p.m.