Least-cost paths and resistance surfaces

landscape ecology

GIS

ecology tutorial

ggplot2

Compute cost distance and least-cost paths in base R with iterative sweeps, then measure how much the answer depends on the resistance values you chose.

Author

Tidy Ecology

Published

2026-07-12

Modified

2026-10-07

Corrected 7 October 2026: the post said the sweep count confirmed that the number of sweeps follows the grid diameter rather than the resistance values, which runs sharing one ranking and one route cannot show; it now says the count follows the number of steps along the cheapest routes, which on this smooth surface is close to the grid diameter, and that a surface whose cheapest routes wind around barriers needs many more sweeps.

Corrected 30 September 2026: the path tracer stepped to the neighbour with the lowest cost distance and ignored the cost of the step, so the paths it drew were not least-cost paths, and the post concluded that the route changed with the resistance contrast. With the step cost included, the traced path costs exactly the cost distance and is the same at every contrast tried. The tracer, the section on route stability, the second figure and the thumbnail are corrected, and the post now says the contrast moved the corridor but not the route over the range tried. The post also said the route held because the best route on this field is already short; it now says the drawn route is the quality-driven one, that it would give way to a route nearer the straight line at contrasts just below the range tried, and why it stays the cheapest within that range.

A resistance surface says how costly it is to cross each cell of a landscape, and a least-cost path is the cheapest route between two points across it. The pair is the workhorse of corridor planning: give the model a source, a destination and a resistance map, and it returns a line on a map that looks authoritative.

The line is only as good as the resistance values, and those are almost always assigned by expert judgement or by a habitat model with no movement data behind it. This tutorial implements cost distance and path tracing from scratch in base R, then measures how much of the answer the resistance values decide.

Cost distance without a package

Cost distance is a shortest-path problem on a grid graph, and it does not need a priority queue if you are willing to sweep. Start with infinite distance everywhere except the source, then repeatedly relax every cell against its eight neighbours until nothing changes. Each sweep is eight shifted-matrix comparisons, which is fast in R, and the number of sweeps needed is roughly the diameter of the grid unless the cheapest routes have to wind around barriers.

library(ggplot2)

te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f",
               clay = "#b5534e", gold = "#cda23f", line = "#dad9ca",
               ink = "#16241d", paper = "#f5f4ee")

theme_te <- function() {
  theme_minimal(base_size = 12) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line(colour = "#e7e6dc"),
          plot.background = element_rect(fill = "#f5f4ee", colour = NA),
          panel.background = element_rect(fill = "#f5f4ee", colour = NA),
          plot.title = element_text(face = "bold", colour = te_pal$ink),
          axis.title = element_text(colour = "#2c3a31"))
}

smooth_field <- function(n, passes, seed) {
  set.seed(seed)
  z <- matrix(runif(n * n), n, n)
  for (p in seq_len(passes)) {
    pad <- rbind(z[1, ], z, z[n, ])
    pad <- cbind(pad[, 1], pad, pad[, n])
    z <- (pad[1:n, 2:(n + 1)] + pad[3:(n + 2), 2:(n + 1)] +
            pad[2:(n + 1), 1:n] + pad[2:(n + 1), 3:(n + 2)] +
            4 * pad[2:(n + 1), 2:(n + 1)]) / 8
  }
  (z - min(z)) / (max(z) - min(z))
}

shift_mat <- function(m, dy, dx) {
  n <- nrow(m)
  p <- ncol(m)
  out <- matrix(Inf, n, p)
  yi <- max(1, 1 - dy):min(n, n - dy)
  xi <- max(1, 1 - dx):min(p, p - dx)
  out[yi, xi] <- m[yi + dy, xi + dx]
  out
}

cost_distance <- function(res, src, max_iter = 2000) {
  d <- matrix(Inf, nrow(res), ncol(res))
  d[src] <- 0
  moves <- list(c(-1, 0), c(1, 0), c(0, -1), c(0, 1),
                c(-1, -1), c(-1, 1), c(1, -1), c(1, 1))
  for (it in seq_len(max_iter)) {
    old <- d
    for (mv in moves) {
      step <- if (sum(abs(mv)) == 2) sqrt(2) else 1
      d <- pmin(d, shift_mat(d, mv[1], mv[2]) +
                  step * (res + shift_mat(res, mv[1], mv[2])) / 2)
    }
    if (max(abs(old - d)[is.finite(old) & is.finite(d)]) < 1e-9 &&
        sum(is.finite(d)) == sum(is.finite(old))) break
  }
  attr(d, "iterations") <- it
  d
}

trace_path <- function(d_from, res, start_idx) {
  n <- nrow(d_from)
  pos <- start_idx
  path <- pos
  moves <- expand.grid(dy = -1:1, dx = -1:1)
  moves <- moves[!(moves$dy == 0 & moves$dx == 0), ]
  while (d_from[pos] > 0) {
    yy <- ((pos - 1) %% n) + 1
    xx <- ((pos - 1) %/% n) + 1
    cand <- data.frame(y = yy + moves$dy, x = xx + moves$dx,
                       len = ifelse(abs(moves$dy) + abs(moves$dx) == 2,
                                    sqrt(2), 1))
    cand <- cand[cand$y >= 1 & cand$y <= n &
                   cand$x >= 1 & cand$x <= ncol(d_from), ]
    idx <- (cand$x - 1) * n + cand$y
    # cost distance of the neighbour plus the cost of the step to it
    total <- d_from[idx] + cand$len * (res[pos] + res[idx]) / 2
    pos <- idx[which.min(total)]
    path <- c(path, pos)
  }
  path
}

path_cost <- function(path, res) {
  n <- nrow(res)
  y <- ((path - 1) %% n) + 1
  x <- ((path - 1) %/% n) + 1
  k <- seq_len(length(path) - 1)
  len <- ifelse(abs(diff(y)) + abs(diff(x)) == 2, sqrt(2), 1)
  sum(len * (res[path[k]] + res[path[k + 1]]) / 2)
}

The cost of stepping between two cells is the mean of their resistances times the step length, one for orthogonal moves and the square root of two for diagonals. That averaging convention matters: taking the destination cell’s resistance instead systematically overcounts the cost of leaving a cheap cell and undercounts the cost of leaving an expensive one, and different software packages make different choices here.

Path tracing is the reverse operation. Once cost distance from the source is known everywhere, walk back from the destination, and at each cell step to the neighbour for which its cost distance plus the cost of the step itself is smallest. Along the least-cost path that sum equals the current cell’s cost distance exactly, which gives a check worth running every time: the summed step costs of the traced path must equal the cost distance at the destination. The tempting shortcut, stepping to whichever neighbour has the lowest cost distance and ignoring what the step costs, fails that check. It happily cuts across an expensive cell to reach a neighbour that is cheap to get to from the source, and it returns a route that is not the cheapest one.

One landscape, four resistance scales

The habitat quality surface is fixed. What varies is how much more expensive the worst cell is than the best, from a factor of five to a factor of a thousand: the same ranking of the landscape, four different contrasts.

n <- 100
quality <- smooth_field(n, 8, seed = 473)
src_a <- matrix(FALSE, n, n)
src_a[8, 8] <- TRUE
src_b <- matrix(FALSE, n, n)
src_b[93, 93] <- TRUE
idx_b <- (93 - 1) * n + 93
diag_cells <- cbind(round(seq(8, 93, length.out = 200)),
                    round(seq(8, 93, length.out = 200)))

contrast_run <- function(rmax) {
  res <- 1 + (rmax - 1) * (1 - quality)
  d_a <- cost_distance(res, src_a)
  d_b <- cost_distance(res, src_b)
  corridor <- d_a + d_b
  list(rmax = rmax, res = res, cost = d_a[idx_b],
       path = trace_path(d_a, res, idx_b), corridor = corridor,
       straight_cost = mean(res[diag_cells]) * sqrt(2) * 85,
       corridor_frac = mean(corridor <= 1.05 * min(corridor)),
       iters = attr(d_a, "iterations"))
}

runs <- lapply(c(5, 20, 100, 1000), contrast_run)
tab <- do.call(rbind, lapply(runs, function(r) data.frame(
  max_resistance = r$rmax, cost_distance = round(r$cost, 1),
  straight_line_cost = round(r$straight_cost, 1),
  detour_saving = round(1 - r$cost / r$straight_cost, 4),
  path_cells = length(r$path),
  path_cost_matches = isTRUE(all.equal(path_cost(r$path, r$res), r$cost)),
  corridor_share = round(r$corridor_frac, 4), sweeps = r$iters)))
print(tab, row.names = FALSE)

 max_resistance cost_distance straight_line_cost detour_saving path_cells
              5         388.4              404.4        0.0396        103
             20        1356.7             1470.1        0.0771        103
            100        6521.0             7153.9        0.0885        103
           1000       64619.8            71096.9        0.0911        103
 path_cost_matches corridor_share sweeps
              TRUE         0.3411    104
              TRUE         0.3052    104
              TRUE         0.3040    104
              TRUE         0.3035    105

Three numbers in that table are worth reading carefully.

The detour saving is what the least-cost route buys over walking the straight line: 3.96 per cent at a contrast of five, rising to 9.11 per cent at a contrast of a thousand. On a smoothly varying surface there is simply not much to gain by going around, and the model’s own optimum is only slightly better than the naive route. That is not a defect of the method; it is information about the landscape, and it is worth computing before a corridor is drawn on the strength of the model.

The sweep count stays near 104 for a 100 by 100 grid, close to the grid diameter. That is a property of this surface rather than of the method: each sweep carries the cost distance at least one step further along every cheapest route, so the count follows the number of steps in the longest cheapest route, and on this smooth surface no cheapest route is much longer than the grid is wide. A resistance surface whose cheapest routes wind back and forth, such as a maze of barriers, would need many more sweeps.

The corridor share is the uncomfortable one. Adding cost distance from both endpoints gives, for every cell, the cost of the cheapest route through that cell, and cells within five per cent of the global minimum form the usual definition of a corridor. That set covers 34.11 per cent of the map at a contrast of five and about 30 per cent at higher contrasts. The single least-cost line is one route through a very wide band of nearly equivalent routes, and drawing it as a line hides that.

r3 <- runs[[3]]
grid_df <- data.frame(
  x = rep(seq_len(n), each = n), y = rep(seq_len(n), times = n),
  resistance = as.vector(r3$res),
  in_corridor = as.vector(r3$corridor <= 1.05 * min(r3$corridor)))
path_df <- data.frame(y = ((r3$path - 1) %% n) + 1,
                      x = ((r3$path - 1) %/% n) + 1)

ggplot(grid_df, aes(x, y)) +
  geom_raster(aes(fill = resistance)) +
  geom_raster(data = subset(grid_df, in_corridor), fill = "white",
              alpha = 0.35) +
  geom_path(data = path_df, colour = te_pal$ink, linewidth = 0.8) +
  coord_equal() +
  scale_fill_gradient(low = te_pal$paper, high = te_pal$forest,
                      name = "resistance") +
  labs(title = "The corridor is wide, the path is one line through it",
       subtitle = "shaded: routes within five per cent of the cheapest total cost",
       x = NULL, y = NULL) +
  theme_te() +
  theme(axis.text = element_blank(), panel.grid = element_blank())

A green-to-cream resistance map with a broad pale corridor running from the lower left to the upper right and a dark line threading through it. — Figure 1: Resistance surface at a contrast of one hundred, with the corridor of routes within five per cent of the cheapest total cost shaded and the least-cost path drawn.

How stable is the route?

If the resistance values are guesses, the question is whether the route survives changing them. Compare the sets of cells the paths occupy, using the share of cells shared between two paths out of the cells in either, and do the same for the corridors.

overlap <- function(a, b) length(intersect(a, b)) / length(union(a, b))
in_corr <- function(r) which(r$corridor <= 1.05 * min(r$corridor))
ov_path_5_1000 <- overlap(runs[[1]]$path, runs[[4]]$path)
ov_path_20_100 <- overlap(runs[[2]]$path, runs[[3]]$path)
ov_corr_5_1000 <- overlap(in_corr(runs[[1]]), in_corr(runs[[4]]))
res_path_100 <- mean(runs[[3]]$res[runs[[3]]$path])
res_map_100 <- mean(runs[[3]]$res)
cat("path overlap, contrast 5 against 1000:", round(ov_path_5_1000, 4),
    " contrast 20 against 100:", round(ov_path_20_100, 4), "\n")

path overlap, contrast 5 against 1000: 1  contrast 20 against 100: 1

cat("corridor overlap, contrast 5 against 1000:", round(ov_corr_5_1000, 4), "\n")

corridor overlap, contrast 5 against 1000: 0.8183

cat("mean resistance along the path at contrast 100:", round(res_path_100, 3),
    " over the whole map:", round(res_map_100, 3), "\n")

mean resistance along the path at contrast 100: 49.941  over the whole map: 58.272

Across the contrasts tried the route does not move. The path at a contrast of five shares 100 per cent of its cells with the path at a contrast of a thousand, and the table shows a path of the same length, 103 cells, at every contrast, while the cost distance grows more than a hundredfold. What does move is the corridor: the routes within five per cent of the optimum overlap by 82 per cent between the lowest and the highest contrast, because the band of near-equivalent routes narrows a little as expensive ground gets more expensive.

That stability is a property of this surface and of the range tried, not a guarantee. Resistance here is one plus a multiple of one minus quality, so the cost of a route is its length plus the contrast minus one times its summed poor quality. A low contrast gives length more say, a high one gives quality more say, and on a landscape where the shortest reasonable route and the best-quality route go different ways the two can change places. They do on this field, just below the range tried: the drawn path is noticeably longer than the straight diagonal, and at contrasts a little under five, length starts to win and the route swings towards the diagonal, while costing nearly the same as the drawn path. Within the range tried the drawn route stays the cheapest, because a route’s cost is linear in the contrast: a route that is cheapest at both five and a thousand is also cheapest at every contrast in between.

The last line points at why the saving stays modest. The mean resistance along the chosen path is 49.94 against 58.27 over the whole map: the least-cost route is cheaper than average terrain, but it cannot dodge expensive ground altogether, because on a smooth surface every route has to cross some of it.

corr_df <- do.call(rbind, lapply(c(1, 4), function(k) {
  data.frame(x = rep(seq_len(n), each = n), y = rep(seq_len(n), times = n),
             quality = as.vector(quality),
             in_corridor = as.vector(runs[[k]]$corridor <=
                                       1.05 * min(runs[[k]]$corridor)),
             contrast = paste("max resistance", runs[[k]]$rmax))
}))
corr_df$contrast <- factor(corr_df$contrast,
                           levels = c("max resistance 5", "max resistance 1000"))
paths_df <- do.call(rbind, lapply(c(1, 4), function(k) {
  p <- runs[[k]]$path
  data.frame(y = ((p - 1) %% n) + 1, x = ((p - 1) %/% n) + 1,
             contrast = paste("max resistance", runs[[k]]$rmax))
}))
paths_df$contrast <- factor(paths_df$contrast, levels = levels(corr_df$contrast))

ggplot(corr_df, aes(x, y)) +
  geom_raster(aes(fill = quality)) +
  geom_raster(data = subset(corr_df, in_corridor), fill = "white",
              alpha = 0.4) +
  geom_path(data = paths_df, colour = te_pal$ink, linewidth = 0.8) +
  facet_wrap(~ contrast) +
  coord_equal() +
  scale_fill_gradient(low = te_pal$forest, high = te_pal$paper,
                      name = "habitat quality") +
  labs(title = "The same route at both contrasts",
       subtitle = "shaded: routes within five per cent of the cheapest total cost",
       x = NULL, y = NULL) +
  theme_te() +
  theme(axis.text = element_blank(), panel.grid = element_blank())

Two panels of the same habitat quality map, each with a pale corridor from the lower left to the upper right and the same dark path through it; the corridor is slightly wider in the left panel, for a contrast of five, than in the right panel, for a contrast of a thousand. — Figure 2: The corridor of routes within five per cent of the cheapest total cost under resistance contrasts of five and one thousand, on the same habitat quality surface, with the least-cost path, which is the same at both contrasts.

What this licenses

Three statements are supported by the machinery above, and a fourth is not.

Cost distance is a legitimate, reproducible summary of a resistance surface: given the surface, the numbers are exact and the algorithm has no tuning parameters. It is a good way to compare candidate routes, to rank pairs of patches by isolation, and to produce a corridor map rather than a line.

The corridor map is more honest than the path, because it shows the width of the set of near-optimal routes. Reporting the path alone implies a precision the model does not have, and the measurement above puts a number on that: about a third of this map lies within five per cent of the optimum.

The contrast of the resistance values is a modelling decision in its own right, separate from their ranking. Here it changed the cost distances by orders of magnitude and the corridor a little, and left the route alone over the range tried, but not below it: the contrast decides how much length is traded for quality, and on this same map a contrast a little below the lowest tried would move the route. Assigning resistances on an arbitrary scale (one to ten, or one to a hundred) should therefore be varied as a sensitivity analysis rather than fixed by convention; it costs one extra run, and it tells you which of your conclusions depend on the scale.

What none of this licenses is the claim that organisms move along the least-cost path. The model optimises total accumulated cost for an omniscient walker; real animals sample locally, take risks, and sometimes cross hostile ground quickly rather than avoiding it. The next tutorial takes the opposite extreme, a walker with no plan at all, and shows that the two models disagree in a way that is itself informative.

References

Adriaensen F, Chardon JP, De Blust G, Swinnen E, Villalba S, Gulinck H, Matthysen E 2003 Landscape and Urban Planning 64(4):233-247 (10.1016/S0169-2046(02)00242-6)

Spear SF, Balkenhol N, Fortin MJ, McRae BH, Scribner K 2010 Molecular Ecology 19(17):3576-3591 (10.1111/j.1365-294X.2010.04657.x)

Sawyer SC, Epps CW, Brashares JS 2011 Journal of Applied Ecology 48(3):668-678 (10.1111/j.1365-2664.2011.01970.x)

Zeller KA, McGarigal K, Whiteley AR 2012 Landscape Ecology 27(6):777-797 (10.1007/s10980-012-9737-0)

--- title: "Least-cost paths and resistance surfaces" description: "Compute cost distance and least-cost paths in base R with iterative sweeps, then measure how much the answer depends on the resistance values you chose." date: "2026-07-12 10:00" categories: [R, landscape ecology, GIS, ecology tutorial, ggplot2] image: thumbnail.png date-modified: "2026-10-07" image-alt: "A habitat quality map with a pale corridor of near-optimal routes and one least-cost path running through it from the lower left to the upper right corner." --- *Corrected 7 October 2026: the post said the sweep count confirmed that the number of sweeps follows the grid diameter rather than the resistance values, which runs sharing one ranking and one route cannot show; it now says the count follows the number of steps along the cheapest routes, which on this smooth surface is close to the grid diameter, and that a surface whose cheapest routes wind around barriers needs many more sweeps.* *Corrected 30 September 2026: the path tracer stepped to the neighbour with the lowest cost distance and ignored the cost of the step, so the paths it drew were not least-cost paths, and the post concluded that the route changed with the resistance contrast. With the step cost included, the traced path costs exactly the cost distance and is the same at every contrast tried. The tracer, the section on route stability, the second figure and the thumbnail are corrected, and the post now says the contrast moved the corridor but not the route over the range tried. The post also said the route held because the best route on this field is already short; it now says the drawn route is the quality-driven one, that it would give way to a route nearer the straight line at contrasts just below the range tried, and why it stays the cheapest within that range.* A resistance surface says how costly it is to cross each cell of a landscape, and a least-cost path is the cheapest route between two points across it. The pair is the workhorse of corridor planning: give the model a source, a destination and a resistance map, and it returns a line on a map that looks authoritative. The line is only as good as the resistance values, and those are almost always assigned by expert judgement or by a habitat model with no movement data behind it. This tutorial implements cost distance and path tracing from scratch in base R, then measures how much of the answer the resistance values decide. ## Cost distance without a package Cost distance is a shortest-path problem on a grid graph, and it does not need a priority queue if you are willing to sweep. Start with infinite distance everywhere except the source, then repeatedly relax every cell against its eight neighbours until nothing changes. Each sweep is eight shifted-matrix comparisons, which is fast in R, and the number of sweeps needed is roughly the diameter of the grid unless the cheapest routes have to wind around barriers. ```{r setup} #| message: false #| warning: false library(ggplot2) te_pal <- list(forest = "#275139", green = "#2f8f63", sage = "#93a87f", clay = "#b5534e", gold = "#cda23f", line = "#dad9ca", ink = "#16241d", paper = "#f5f4ee") theme_te <- function() { theme_minimal(base_size = 12) + theme(panel.grid.minor = element_blank(), panel.grid.major = element_line(colour = "#e7e6dc"), plot.background = element_rect(fill = "#f5f4ee", colour = NA), panel.background = element_rect(fill = "#f5f4ee", colour = NA), plot.title = element_text(face = "bold", colour = te_pal$ink), axis.title = element_text(colour = "#2c3a31")) } ``` ```{r engine} smooth_field <- function(n, passes, seed) { set.seed(seed) z <- matrix(runif(n * n), n, n) for (p in seq_len(passes)) { pad <- rbind(z[1, ], z, z[n, ]) pad <- cbind(pad[, 1], pad, pad[, n]) z <- (pad[1:n, 2:(n + 1)] + pad[3:(n + 2), 2:(n + 1)] + pad[2:(n + 1), 1:n] + pad[2:(n + 1), 3:(n + 2)] + 4 * pad[2:(n + 1), 2:(n + 1)]) / 8 } (z - min(z)) / (max(z) - min(z)) } shift_mat <- function(m, dy, dx) { n <- nrow(m) p <- ncol(m) out <- matrix(Inf, n, p) yi <- max(1, 1 - dy):min(n, n - dy) xi <- max(1, 1 - dx):min(p, p - dx) out[yi, xi] <- m[yi + dy, xi + dx] out } cost_distance <- function(res, src, max_iter = 2000) { d <- matrix(Inf, nrow(res), ncol(res)) d[src] <- 0 moves <- list(c(-1, 0), c(1, 0), c(0, -1), c(0, 1), c(-1, -1), c(-1, 1), c(1, -1), c(1, 1)) for (it in seq_len(max_iter)) { old <- d for (mv in moves) { step <- if (sum(abs(mv)) == 2) sqrt(2) else 1 d <- pmin(d, shift_mat(d, mv[1], mv[2]) + step * (res + shift_mat(res, mv[1], mv[2])) / 2) } if (max(abs(old - d)[is.finite(old) & is.finite(d)]) < 1e-9 && sum(is.finite(d)) == sum(is.finite(old))) break } attr(d, "iterations") <- it d } trace_path <- function(d_from, res, start_idx) { n <- nrow(d_from) pos <- start_idx path <- pos moves <- expand.grid(dy = -1:1, dx = -1:1) moves <- moves[!(moves$dy == 0 & moves$dx == 0), ] while (d_from[pos] > 0) { yy <- ((pos - 1) %% n) + 1 xx <- ((pos - 1) %/% n) + 1 cand <- data.frame(y = yy + moves$dy, x = xx + moves$dx, len = ifelse(abs(moves$dy) + abs(moves$dx) == 2, sqrt(2), 1)) cand <- cand[cand$y >= 1 & cand$y <= n & cand$x >= 1 & cand$x <= ncol(d_from), ] idx <- (cand$x - 1) * n + cand$y # cost distance of the neighbour plus the cost of the step to it total <- d_from[idx] + cand$len * (res[pos] + res[idx]) / 2 pos <- idx[which.min(total)] path <- c(path, pos) } path } path_cost <- function(path, res) { n <- nrow(res) y <- ((path - 1) %% n) + 1 x <- ((path - 1) %/% n) + 1 k <- seq_len(length(path) - 1) len <- ifelse(abs(diff(y)) + abs(diff(x)) == 2, sqrt(2), 1) sum(len * (res[path[k]] + res[path[k + 1]]) / 2) } ``` The cost of stepping between two cells is the mean of their resistances times the step length, one for orthogonal moves and the square root of two for diagonals. That averaging convention matters: taking the destination cell's resistance instead systematically overcounts the cost of leaving a cheap cell and undercounts the cost of leaving an expensive one, and different software packages make different choices here. Path tracing is the reverse operation. Once cost distance from the source is known everywhere, walk back from the destination, and at each cell step to the neighbour for which its cost distance plus the cost of the step itself is smallest. Along the least-cost path that sum equals the current cell's cost distance exactly, which gives a check worth running every time: the summed step costs of the traced path must equal the cost distance at the destination. The tempting shortcut, stepping to whichever neighbour has the lowest cost distance and ignoring what the step costs, fails that check. It happily cuts across an expensive cell to reach a neighbour that is cheap to get to from the source, and it returns a route that is not the cheapest one. ## One landscape, four resistance scales The habitat quality surface is fixed. What varies is how much more expensive the worst cell is than the best, from a factor of five to a factor of a thousand: the same ranking of the landscape, four different contrasts. ```{r runs} n <- 100 quality <- smooth_field(n, 8, seed = 473) src_a <- matrix(FALSE, n, n) src_a[8, 8] <- TRUE src_b <- matrix(FALSE, n, n) src_b[93, 93] <- TRUE idx_b <- (93 - 1) * n + 93 diag_cells <- cbind(round(seq(8, 93, length.out = 200)), round(seq(8, 93, length.out = 200))) contrast_run <- function(rmax) { res <- 1 + (rmax - 1) * (1 - quality) d_a <- cost_distance(res, src_a) d_b <- cost_distance(res, src_b) corridor <- d_a + d_b list(rmax = rmax, res = res, cost = d_a[idx_b], path = trace_path(d_a, res, idx_b), corridor = corridor, straight_cost = mean(res[diag_cells]) * sqrt(2) * 85, corridor_frac = mean(corridor <= 1.05 * min(corridor)), iters = attr(d_a, "iterations")) } runs <- lapply(c(5, 20, 100, 1000), contrast_run) tab <- do.call(rbind, lapply(runs, function(r) data.frame( max_resistance = r$rmax, cost_distance = round(r$cost, 1), straight_line_cost = round(r$straight_cost, 1), detour_saving = round(1 - r$cost / r$straight_cost, 4), path_cells = length(r$path), path_cost_matches = isTRUE(all.equal(path_cost(r$path, r$res), r$cost)), corridor_share = round(r$corridor_frac, 4), sweeps = r$iters))) print(tab, row.names = FALSE) ``` Three numbers in that table are worth reading carefully. The **detour saving** is what the least-cost route buys over walking the straight line: `r sprintf("%.2f", 100 * (1 - runs[[1]]$cost / runs[[1]]$straight_cost))` per cent at a contrast of five, rising to `r sprintf("%.2f", 100 * (1 - runs[[4]]$cost / runs[[4]]$straight_cost))` per cent at a contrast of a thousand. On a smoothly varying surface there is simply not much to gain by going around, and the model's own optimum is only slightly better than the naive route. That is not a defect of the method; it is information about the landscape, and it is worth computing before a corridor is drawn on the strength of the model. The **sweep count** stays near 104 for a 100 by 100 grid, close to the grid diameter. That is a property of this surface rather than of the method: each sweep carries the cost distance at least one step further along every cheapest route, so the count follows the number of steps in the longest cheapest route, and on this smooth surface no cheapest route is much longer than the grid is wide. A resistance surface whose cheapest routes wind back and forth, such as a maze of barriers, would need many more sweeps. The **corridor share** is the uncomfortable one. Adding cost distance from both endpoints gives, for every cell, the cost of the cheapest route through that cell, and cells within five per cent of the global minimum form the usual definition of a corridor. That set covers `r sprintf("%.2f", 100 * runs[[1]]$corridor_frac)` per cent of the map at a contrast of five and about 30 per cent at higher contrasts. The single least-cost line is one route through a very wide band of nearly equivalent routes, and drawing it as a line hides that. ```{r fig-corridor} #| fig-width: 7 #| fig-height: 4.6 #| fig-cap: "Resistance surface at a contrast of one hundred, with the corridor of routes within five per cent of the cheapest total cost shaded and the least-cost path drawn." #| fig-alt: "A green-to-cream resistance map with a broad pale corridor running from the lower left to the upper right and a dark line threading through it." r3 <- runs[[3]] grid_df <- data.frame( x = rep(seq_len(n), each = n), y = rep(seq_len(n), times = n), resistance = as.vector(r3$res), in_corridor = as.vector(r3$corridor <= 1.05 * min(r3$corridor))) path_df <- data.frame(y = ((r3$path - 1) %% n) + 1, x = ((r3$path - 1) %/% n) + 1) ggplot(grid_df, aes(x, y)) + geom_raster(aes(fill = resistance)) + geom_raster(data = subset(grid_df, in_corridor), fill = "white", alpha = 0.35) + geom_path(data = path_df, colour = te_pal$ink, linewidth = 0.8) + coord_equal() + scale_fill_gradient(low = te_pal$paper, high = te_pal$forest, name = "resistance") + labs(title = "The corridor is wide, the path is one line through it", subtitle = "shaded: routes within five per cent of the cheapest total cost", x = NULL, y = NULL) + theme_te() + theme(axis.text = element_blank(), panel.grid = element_blank()) ``` ## How stable is the route? If the resistance values are guesses, the question is whether the route survives changing them. Compare the sets of cells the paths occupy, using the share of cells shared between two paths out of the cells in either, and do the same for the corridors. ```{r overlap} overlap <- function(a, b) length(intersect(a, b)) / length(union(a, b)) in_corr <- function(r) which(r$corridor <= 1.05 * min(r$corridor)) ov_path_5_1000 <- overlap(runs[[1]]$path, runs[[4]]$path) ov_path_20_100 <- overlap(runs[[2]]$path, runs[[3]]$path) ov_corr_5_1000 <- overlap(in_corr(runs[[1]]), in_corr(runs[[4]])) res_path_100 <- mean(runs[[3]]$res[runs[[3]]$path]) res_map_100 <- mean(runs[[3]]$res) cat("path overlap, contrast 5 against 1000:", round(ov_path_5_1000, 4), " contrast 20 against 100:", round(ov_path_20_100, 4), "\n") cat("corridor overlap, contrast 5 against 1000:", round(ov_corr_5_1000, 4), "\n") cat("mean resistance along the path at contrast 100:", round(res_path_100, 3), " over the whole map:", round(res_map_100, 3), "\n") ``` Across the contrasts tried the route does not move. The path at a contrast of five shares `r sprintf("%.0f", 100 * ov_path_5_1000)` per cent of its cells with the path at a contrast of a thousand, and the table shows a path of the same length, `r length(runs[[1]]$path)` cells, at every contrast, while the cost distance grows more than a hundredfold. What does move is the corridor: the routes within five per cent of the optimum overlap by `r sprintf("%.0f", 100 * ov_corr_5_1000)` per cent between the lowest and the highest contrast, because the band of near-equivalent routes narrows a little as expensive ground gets more expensive. That stability is a property of this surface and of the range tried, not a guarantee. Resistance here is one plus a multiple of one minus quality, so the cost of a route is its length plus the contrast minus one times its summed poor quality. A low contrast gives length more say, a high one gives quality more say, and on a landscape where the shortest reasonable route and the best-quality route go different ways the two can change places. They do on this field, just below the range tried: the drawn path is noticeably longer than the straight diagonal, and at contrasts a little under five, length starts to win and the route swings towards the diagonal, while costing nearly the same as the drawn path. Within the range tried the drawn route stays the cheapest, because a route's cost is linear in the contrast: a route that is cheapest at both five and a thousand is also cheapest at every contrast in between. The last line points at why the saving stays modest. The mean resistance along the chosen path is `r sprintf("%.2f", res_path_100)` against `r sprintf("%.2f", res_map_100)` over the whole map: the least-cost route is cheaper than average terrain, but it cannot dodge expensive ground altogether, because on a smooth surface every route has to cross some of it. ```{r fig-paths} #| fig-width: 7 #| fig-height: 4.2 #| fig-cap: "The corridor of routes within five per cent of the cheapest total cost under resistance contrasts of five and one thousand, on the same habitat quality surface, with the least-cost path, which is the same at both contrasts." #| fig-alt: "Two panels of the same habitat quality map, each with a pale corridor from the lower left to the upper right and the same dark path through it; the corridor is slightly wider in the left panel, for a contrast of five, than in the right panel, for a contrast of a thousand." corr_df <- do.call(rbind, lapply(c(1, 4), function(k) { data.frame(x = rep(seq_len(n), each = n), y = rep(seq_len(n), times = n), quality = as.vector(quality), in_corridor = as.vector(runs[[k]]$corridor <= 1.05 * min(runs[[k]]$corridor)), contrast = paste("max resistance", runs[[k]]$rmax)) })) corr_df$contrast <- factor(corr_df$contrast, levels = c("max resistance 5", "max resistance 1000")) paths_df <- do.call(rbind, lapply(c(1, 4), function(k) { p <- runs[[k]]$path data.frame(y = ((p - 1) %% n) + 1, x = ((p - 1) %/% n) + 1, contrast = paste("max resistance", runs[[k]]$rmax)) })) paths_df$contrast <- factor(paths_df$contrast, levels = levels(corr_df$contrast)) ggplot(corr_df, aes(x, y)) + geom_raster(aes(fill = quality)) + geom_raster(data = subset(corr_df, in_corridor), fill = "white", alpha = 0.4) + geom_path(data = paths_df, colour = te_pal$ink, linewidth = 0.8) + facet_wrap(~ contrast) + coord_equal() + scale_fill_gradient(low = te_pal$forest, high = te_pal$paper, name = "habitat quality") + labs(title = "The same route at both contrasts", subtitle = "shaded: routes within five per cent of the cheapest total cost", x = NULL, y = NULL) + theme_te() + theme(axis.text = element_blank(), panel.grid = element_blank()) ``` ## What this licenses Three statements are supported by the machinery above, and a fourth is not. Cost distance is a legitimate, reproducible summary of a resistance surface: given the surface, the numbers are exact and the algorithm has no tuning parameters. It is a good way to compare candidate routes, to rank pairs of patches by isolation, and to produce a corridor map rather than a line. The corridor map is more honest than the path, because it shows the width of the set of near-optimal routes. Reporting the path alone implies a precision the model does not have, and the measurement above puts a number on that: about a third of this map lies within five per cent of the optimum. The contrast of the resistance values is a modelling decision in its own right, separate from their ranking. Here it changed the cost distances by orders of magnitude and the corridor a little, and left the route alone over the range tried, but not below it: the contrast decides how much length is traded for quality, and on this same map a contrast a little below the lowest tried would move the route. Assigning resistances on an arbitrary scale (one to ten, or one to a hundred) should therefore be varied as a sensitivity analysis rather than fixed by convention; it costs one extra run, and it tells you which of your conclusions depend on the scale. What none of this licenses is the claim that organisms move along the least-cost path. The model optimises total accumulated cost for an omniscient walker; real animals sample locally, take risks, and sometimes cross hostile ground quickly rather than avoiding it. The next tutorial takes the opposite extreme, a walker with no plan at all, and shows that the two models disagree in a way that is itself informative. ## References Adriaensen F, Chardon JP, De Blust G, Swinnen E, Villalba S, Gulinck H, Matthysen E 2003 Landscape and Urban Planning 64(4):233-247 (10.1016/S0169-2046(02)00242-6) Spear SF, Balkenhol N, Fortin MJ, McRae BH, Scribner K 2010 Molecular Ecology 19(17):3576-3591 (10.1111/j.1365-294X.2010.04657.x) Sawyer SC, Epps CW, Brashares JS 2011 Journal of Applied Ecology 48(3):668-678 (10.1111/j.1365-2664.2011.01970.x) Zeller KA, McGarigal K, Whiteley AR 2012 Landscape Ecology 27(6):777-797 (10.1007/s10980-012-9737-0) ## Related tutorials - [Patch metrics and fragmentation](../patch-metrics-and-fragmentation/) - [Circuit theory and effective resistance](../circuit-theory-and-effective-resistance/) - [Checking a connectivity analysis](../checking-a-connectivity-analysis/) - [Isolation by distance and by resistance](../isolation-by-resistance/)