This vignette is a diagnostic walkthrough for the example tsunami-evacuation workflow. It is intentionally more detailed than a normal user vignette because its purpose is to test every major evacpath function independently before trusting the high-level run_evacpath() wrapper.

The workflow builds on open-source least-cost path methods for evacuation planning (Cordero et al. 2025) and uses leastcostpath for least-cost path and movement-potential modeling (Lewis 2023; Lewis 2021).

1. Setup

# This loader works both when the package is installed and during local package development.
if (requireNamespace("evacpath", quietly = TRUE)) {
  library(evacpath)
} else if (requireNamespace("devtools", quietly = TRUE) && file.exists("DESCRIPTION")) {
  devtools::load_all(".")
} else if (requireNamespace("devtools", quietly = TRUE) && file.exists(file.path("evacpath", "DESCRIPTION"))) {
  devtools::load_all("evacpath")
} else {
  stop("Install evacpath or run this from the package root/parent directory.", call. = FALSE)
}

library(terra)
#> terra 1.9.1

run_diagnostics <- isTRUE(params$run_diagnostics)
run_lcp <- isTRUE(params$run_lcp)

assert_true <- function(x, message) {
  if (!isTRUE(x)) {
    stop(message, call. = FALSE)
  }
  invisible(TRUE)
}

is_spatvector <- function(x) inherits(x, "SpatVector")
is_spatraster <- function(x) inherits(x, "SpatRaster")

count_features <- function(x) {
  if (inherits(x, "SpatRaster")) return(terra::ncell(x))
  if (inherits(x, "SpatVector")) return(nrow(x))
  NA_integer_
}

The diagnostic switches are:

run_diagnostics
#> [1] FALSE
run_lcp
#> [1] FALSE

2. Locate example input files

The diagnostic expects these files:

data-raw/dem.tif
data-raw/rds.gpkg
data-raw/tsunami_inundation_depth.tif

The function below searches for them from the package root, the parent folder, or installed package extdata.

find_input_file <- function(filename, data_dir = params$data_dir) {
  candidate_dirs <- unique(c(
    data_dir,
    file.path("evacpath", data_dir),
    file.path(getwd(), data_dir),
    file.path(getwd(), "evacpath", data_dir),
    system.file("extdata", package = "evacpath")
  ))

  candidate_dirs <- candidate_dirs[nzchar(candidate_dirs)]
  candidates <- file.path(candidate_dirs, filename)
  hits <- candidates[file.exists(candidates)]

  if (length(hits) == 0L) {
    return(NA_character_)
  }

  normalizePath(hits[1])
}

dem_file <- find_input_file("dem.tif")
roads_file <- find_input_file("rds.gpkg")
inundation_file <- find_input_file("tsunami_inundation_depth.tif")

input_files <- data.frame(
  layer = c("DEM", "roads", "inundation depth"),
  file = c(dem_file, roads_file, inundation_file),
  exists = file.exists(c(dem_file, roads_file, inundation_file))
)

input_files
#>              layer
#> 1              DEM
#> 2            roads
#> 3 inundation depth
#>                                                                                                                                  file
#> 1                      /private/var/folders/7j/dr505g_j3zd9z6m9qdykzc4w0000gn/T/RtmprRFJQI/Rinst33ee1efff34d/evacpath/extdata/dem.tif
#> 2                     /private/var/folders/7j/dr505g_j3zd9z6m9qdykzc4w0000gn/T/RtmprRFJQI/Rinst33ee1efff34d/evacpath/extdata/rds.gpkg
#> 3 /private/var/folders/7j/dr505g_j3zd9z6m9qdykzc4w0000gn/T/RtmprRFJQI/Rinst33ee1efff34d/evacpath/extdata/tsunami_inundation_depth.tif
#>   exists
#> 1   TRUE
#> 2   TRUE
#> 3   TRUE

example_data_available <- all(input_files$exists)
example_data_available
#> [1] TRUE

If example_data_available is FALSE, check params or use the packaged files in inst/extdata/.

3. Read the spatial inputs with `read_spatial()`

read_spatial() accepts a file path or an existing terra object. Raster formats are read with terra::rast() and vector formats are read with terra::vect().

dem_raw <- read_spatial(dem_file)
roads_raw <- read_spatial(roads_file)
inundation_raw <- read_spatial(inundation_file)

assert_true(is_spatraster(dem_raw), "DEM was not read as a SpatRaster.")
assert_true(is_spatvector(roads_raw), "Roads were not read as a SpatVector.")
assert_true(is_spatraster(inundation_raw), "Inundation input was not read as a SpatRaster.")

list(
  dem_class = class(dem_raw),
  roads_class = class(roads_raw),
  inundation_class = class(inundation_raw)
)

3.1 Raw input diagnostics

Before modeling, inspect CRS, extent, resolution, raster value ranges, and road attributes.

cat("DEM CRS:\n", terra::crs(dem_raw), "\n\n")
cat("Roads CRS:\n", terra::crs(roads_raw), "\n\n")
cat("Inundation CRS:\n", terra::crs(inundation_raw), "\n\n")

cat("DEM resolution:\n")
print(terra::res(dem_raw))

cat("Inundation resolution:\n")
print(terra::res(inundation_raw))

cat("DEM cell count:\n")
print(terra::ncell(dem_raw))

cat("Road feature count:\n")
print(nrow(roads_raw))

cat("Road attribute names:\n")
print(names(roads_raw))

cat("DEM range:\n")
print(terra::global(dem_raw, fun = range, na.rm = TRUE))

cat("Inundation range:\n")
print(terra::global(inundation_raw, fun = range, na.rm = TRUE))

oldpar <- par(mfrow = c(1, 2), mar = c(3, 3, 3, 5), las = 1)
plot(dem_raw, main = "Raw DEM")
plot(inundation_raw, main = "Raw inundation depth")
par(oldpar)

4. Clean roads with `clean_roads()`

This is where region-specific road filtering belongs. In the packaged example, piers are excluded because they can create misleading road-water intersections.

road_exclude <- NULL
if (!is.null(params$road_exclude_field) && nzchar(params$road_exclude_field)) {
  road_exclude <- list(
    field = params$road_exclude_field,
    values = params$road_exclude_values
  )
}

roads_clean <- clean_roads(
  roads = roads_raw,
  exclude = road_exclude,
  target_crs = params$target_crs
)

assert_true(is_spatvector(roads_clean), "clean_roads() did not return a SpatVector.")
assert_true(nrow(roads_clean) > 0, "clean_roads() returned no roads.")

nrow(roads_raw)
nrow(roads_clean)

5. Prepare tsunami zones with `prepare_tsunami_zones()`

This is the critical tsunami-specific step. The old script created a land mask, water mask, land-only TEZ, and then a TEZ-plus-water layer. The package makes that separation explicit.

zones <- prepare_tsunami_zones(
  inundation = inundation_raw,
  dem = dem_raw,
  target_crs = params$target_crs,
  inundation_threshold = params$inundation_threshold,
  land_threshold = params$land_threshold,
  water_threshold = params$water_threshold,
  dem_sign_multiplier = params$dem_sign_multiplier,
  resample_method = params$resample_method,
  as_polygon = TRUE,
  dissolve = TRUE
)

zones

assert_true(is_spatvector(zones$hazard_zone), "zones$hazard_zone is not a SpatVector.")
assert_true(is_spatvector(zones$escape_zone), "zones$escape_zone is not a SpatVector.")
assert_true(is_spatraster(zones$hazard_raster), "zones$hazard_raster is not a SpatRaster.")
assert_true(is_spatraster(zones$escape_raster), "zones$escape_raster is not a SpatRaster.")

cat("Land-only hazard-zone polygons:", nrow(zones$hazard_zone), "\n")
cat("Water-combined escape-zone polygons:", nrow(zones$escape_zone), "\n")

5.1 Check the land/water classification

If land and water are reversed, change dem_sign_multiplier from 1 to -1 or adjust the thresholds.

oldpar <- par(mfrow = c(1, 2), mar = c(3, 3, 3, 5), las = 1)
plot(zones$land_mask, col = "grey85", legend = FALSE, main = "DEM-derived land mask")
plot(zones$water_mask, col = "lightskyblue1", legend = FALSE, main = "DEM-derived water mask")
par(oldpar)

5.2 Check the two tsunami zones

The land-only zone is for origins and final mapped evacuation time. The water-combined zone is the base zone for escape-point detection.

oldpar <- par(mfrow = c(1, 2), mar = c(3, 3, 3, 5), las = 1)
plot(zones$hazard_zone, col = "tomato", border = NA, main = "hazard_zone: land-only TEZ")
plot(zones$escape_zone, col = "grey90", border = NA, main = "escape_zone: TEZ + water")
plot(zones$hazard_zone, add = TRUE, col = "tomato", border = NA)
par(oldpar)

6. Crop roads to an inner extent before escape-point detection

If the road dataset extends beyond the reliable hazard-zone coverage, roads can intersect the artificial rectangular study-area boundary. Those intersections are not real escape points. The fix is to crop the roads used for escape-point detection to a slightly inset extent. The x and y insets are controlled separately.

roads_for_escape <- crop_roads_to_inner_extent(
  roads = roads_clean,
  zone = zones$escape_zone,
  inset_x_m = params$escape_roads_inset_x_m,
  inset_y_m = params$escape_roads_inset_y_m
)

assert_true(is_spatvector(roads_for_escape), "roads_for_escape is not a SpatVector.")
assert_true(nrow(roads_for_escape) > 0, "roads_for_escape has no features.")

cat("Clean roads:", nrow(roads_clean), "\n")
cat("Roads used for escape-point detection:", nrow(roads_for_escape), "\n")

plot(zones$escape_zone, col = "grey92", border = "grey60", main = "Roads cropped to inset extent for escape detection")
plot(roads_clean, add = TRUE, col = adjustcolor("red", alpha.f = 0.25), lwd = 0.35)
plot(roads_for_escape, add = TRUE, col = "black", lwd = 0.45)
legend(
  "topright",
  legend = c("Original cleaned roads", "Roads used for escape detection"),
  col = c(adjustcolor("red", alpha.f = 0.5), "black"),
  lty = 1,
  lwd = 1,
  bty = "n",
  cex = 0.8
)

7. Combine buffered roads with the escape zone

This is the second major correction. The buffered roads should be combined with the tsunami escape zone before escape points are generated. Otherwise, bridge, causeway, or walkway segments over water can be lost because the land and water masks do not always represent the accessible transportation surface.

The package helper is make_road_aware_escape_zone().

escape_boundary_zone <- make_road_aware_escape_zone(
  escape_zone = zones$escape_zone,
  roads = roads_for_escape,
  road_buffer_m = params$road_buffer_m,
  crop_buffer_m = params$final_road_buffer_m,
  include_base_zone = TRUE
)

assert_true(is_spatvector(escape_boundary_zone), "escape_boundary_zone is not a SpatVector.")
assert_true(nrow(escape_boundary_zone) > 0, "escape_boundary_zone has no features.")

nrow(escape_boundary_zone)

plot(escape_boundary_zone, col = "grey90", border = NA, main = "Road-aware escape-boundary zone")
plot(zones$hazard_zone, add = TRUE, col = adjustcolor("tomato", alpha.f = 0.45), border = NA)
plot(roads_for_escape, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.45), lwd = 0.35)

8. Compare escape-point methods

This section intentionally compares three methods:

Wrong: land-only hazard_zone.
Better: TEZ + water escape_zone.
Best for this workflow: road-aware escape boundary plus inset-cropped roads.

escape_wrong <- try(
  find_escape_points(
    hazard_zone = zones$hazard_zone,
    roads = roads_clean
  ),
  silent = TRUE
)

escape_water <- find_escape_points(
  hazard_zone = zones$escape_zone,
  roads = roads_for_escape
)

escape_points <- find_escape_points(
  hazard_zone = escape_boundary_zone,
  roads = roads_for_escape
)

comparison <- data.frame(
  method = c(
    "wrong: land-only hazard zone",
    "better: TEZ + water, inset roads",
    "preferred: road-aware escape zone, inset roads"
  ),
  escape_point_count = c(
    if (inherits(escape_wrong, "try-error")) NA_integer_ else nrow(escape_wrong),
    nrow(escape_water),
    nrow(escape_points)
  )
)
comparison

assert_true(nrow(escape_points) > 0, "Preferred escape-point calculation returned no points.")

oldpar <- par(mfrow = c(1, 3))

plot(zones$hazard_zone, col = "grey90", border = "grey40", main = "Wrong\nland-only boundary")
plot(roads_clean, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.35), lwd = 0.35)
if (!inherits(escape_wrong, "try-error")) {
  plot(escape_wrong, add = TRUE, pch = 22, bg = "red", col = "black", cex = 0.45)
}

plot(zones$escape_zone, col = "grey90", border = "grey40", main = "Better\nTEZ + water")
plot(roads_for_escape, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.35), lwd = 0.35)
plot(escape_water, add = TRUE, pch = 22, bg = "red", col = "black", cex = 0.45)

plot(escape_boundary_zone, col = "grey90", border = "grey40", main = "Preferred\nroad-aware boundary")
plot(roads_for_escape, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.35), lwd = 0.35)
plot(escape_points, add = TRUE, pch = 22, bg = "red", col = "black", cex = 0.45)

par(oldpar)

From this point forward, the diagnostic uses:

hazard_zone          = zones$hazard_zone
escape_zone          = zones$escape_zone
roads_for_escape     = inset-cropped roads
escape_boundary_zone = road-aware TEZ + water boundary
escape_points        = preferred escape points

hazard_zone <- zones$hazard_zone
escape_zone <- zones$escape_zone
dem <- zones$dem
roads <- roads_clean

9. Prepare core inputs with `prepare_evac_inputs()`

prepare_evac_inputs() reads, cleans, projects, and converts the land-only hazard_zone into a polygon if needed. It does not create the road-aware escape boundary; that is handled separately above.

inputs <- prepare_evac_inputs(
  hazard_zone = hazard_zone,
  roads = roads,
  dem = dem,
  target_crs = params$target_crs,
  hazard_as_polygon = TRUE,
  dissolve_hazard = TRUE,
  road_exclude = NULL
)

hazard_zone <- inputs$hazard_zone
roads <- inputs$roads
dem <- inputs$dem

assert_true(is_spatvector(hazard_zone), "Prepared hazard_zone is not a SpatVector.")
assert_true(is_spatvector(roads), "Prepared roads are not a SpatVector.")
assert_true(is_spatraster(dem), "Prepared DEM is not a SpatRaster.")

10. Create the evacuation grid with `make_evac_grid()`

The evacuation grid should be created from the land-only hazard_zone, not the water-combined or road-aware escape zone.

grid_resolution <- terra::res(dem) * params$grid_multiplier

evac_grid <- make_evac_grid(
  hazard_zone = hazard_zone,
  resolution = grid_resolution
)

assert_true(is_spatvector(evac_grid), "make_evac_grid() did not return a SpatVector.")
assert_true(nrow(evac_grid) > 0, "make_evac_grid() returned no grid cells.")

cat("Grid resolution:", paste(grid_resolution, collapse = " x "), "\n")
cat("Grid cells:", nrow(evac_grid), "\n")

plot(hazard_zone, col = "grey95", border = "grey50", main = "Evacuation grid over land-only TEZ")
plot(evac_grid, add = TRUE, border = adjustcolor("black", alpha.f = 0.25), col = NA)
plot(roads, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.45), lwd = 0.35)

11. Create the road mask with `make_road_mask()`

make_road_mask() buffers the roads and preferred escape points, then combines them. The mask constrains least-cost movement to the road/pathway network while still allowing access around escape points.

road_mask_parts <- make_road_mask(
  roads = roads,
  escape_points = escape_points,
  road_buffer_m = params$road_buffer_m,
  escape_buffer_m = params$escape_buffer_m,
  return_components = TRUE
)

road_mask <- road_mask_parts$mask
roads_buffer <- road_mask_parts$roads_buffer
escape_buffer <- road_mask_parts$escape_buffer

assert_true(is_spatvector(road_mask), "Road mask is not a SpatVector.")
assert_true(is_spatvector(roads_buffer), "roads_buffer is not a SpatVector.")
assert_true(is_spatvector(escape_buffer), "escape_buffer is not a SpatVector.")

plot(hazard_zone, col = "grey95", border = "grey60", main = "Road-constrained movement mask")
plot(road_mask, add = TRUE, col = adjustcolor("grey30", alpha.f = 0.35), border = NA)
plot(escape_points, add = TRUE, pch = 22, bg = "red", col = "black", cex = 0.45)

12. Create road origin points with `make_road_origins()`

Origins should be road-based points inside the land-only inundation zone.

road_points <- make_road_origins(
  evac_grid = evac_grid,
  roads_buffer = roads_buffer,
  hazard_zone = hazard_zone,
  max_origins = params$max_origins,
  seed = params$seed
)

assert_true(is_spatvector(road_points), "make_road_origins() did not return a SpatVector.")
assert_true(nrow(road_points) > 0, "make_road_origins() returned no origins.")

cat("Road origin points retained:", nrow(road_points), "\n")

plot(hazard_zone, col = "grey95", border = "grey60", main = "Road origin points inside land-only TEZ")
plot(roads, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.35), lwd = 0.35)
plot(road_points, add = TRUE, pch = 21, bg = "dodgerblue", col = "black", cex = 0.45)
plot(escape_points, add = TRUE, pch = 22, bg = "red", col = "black", cex = 0.45)

13. Create the conductance surface with `make_conductance_surface()`

The conductance surface is built from the DEM but masked to the road-constrained movement area. This is the expensive part. dem_resolution can be coarsened for testing.

dem_resolution <- terra::res(dem) * params$dem_multiplier

conductance <- make_conductance_surface(
  dem = dem,
  road_mask = road_mask,
  resolution = dem_resolution
)

cs_raster <- leastcostpath::rasterise(conductance)

assert_true(is_spatraster(cs_raster), "Conductance rasterization failed.")

cat("Conductance raster resolution:", paste(terra::res(cs_raster), collapse = " x "), "\n")
cat("Conductance raster cells:", terra::ncell(cs_raster), "\n")

plot(cs_raster, main = "Road-constrained conductance surface")
plot(hazard_zone, add = TRUE, border = "black", col = NA)
plot(roads, add = TRUE, col = adjustcolor("grey20", alpha.f = 1), lwd = 0.35)
plot(escape_points, add = TRUE, pch = 22, bg = "red", col = "black", cex = 0.45)

14. Calculate least-cost distance with `calc_min_distance_to_safety()`

This section is optional because it can be slow for large study areas. The packaged example area is small, so the default parameters use all available origins and destinations.

escape_points_ltd <- escape_points
if (!is.null(params$max_destinations) && params$max_destinations < nrow(escape_points_ltd)) {
  set.seed(params$seed)
  escape_points_ltd <- escape_points_ltd[sort(sample.int(nrow(escape_points_ltd), params$max_destinations)), ]
}

cat("Escape points available:", nrow(escape_points), "\n")
cat("Escape points used in diagnostic LCP:", nrow(escape_points_ltd), "\n")

distance_points <- calc_min_distance_to_safety(
  cs = conductance,
  origins = road_points,
  destinations = escape_points_ltd,
  include_destinations = TRUE,
  progress = TRUE,
  progress_every = 500,
  check_locations = FALSE
)

assert_true(is_spatvector(distance_points), "calc_min_distance_to_safety() did not return a SpatVector.")
assert_true("distance" %in% names(distance_points), "distance_points does not contain a distance field.")

summary(distance_points$distance)
table(distance_points$type)

plot(hazard_zone, col = "grey95", border = "grey60", main = "Minimum least-cost distance points")
plot(roads, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.30), lwd = 0.35)
plot(distance_points, "distance", add = TRUE)

15. Convert distance to time with `calc_evac_time()`

calc_evac_time() converts distance to time using a walking speed. The default paper-style value is 1.22 m/s, but local planning scenarios should modify it.

example_distances <- c(0, 100, 500, 1000, 2000)
data.frame(
  distance_m = example_distances,
  time_min = calc_evac_time(
    distance_m = example_distances,
    walking_speed_mps = params$walking_speed_mps,
    units = "minutes"
  )
)

16. Create evacuation polygons with `make_evac_polygons()`

make_evac_polygons() converts the distance points to Voronoi polygons, calculates evacuation time, and clips the result to the selected output area. For the broad time grid, the clip area should be the land-only hazard_zone.

evac_polygons <- make_evac_polygons(
  distance_points = distance_points,
  clip_area = hazard_zone,
  walking_speed_mps = params$walking_speed_mps,
  region_name = params$region_name
)

assert_true(is_spatvector(evac_polygons), "make_evac_polygons() did not return a SpatVector.")
assert_true("EvacTimeAvg" %in% names(evac_polygons), "EvacTimeAvg field is missing.")

summary(evac_polygons$EvacTimeAvg)

plot(evac_polygons, "EvacTimeAvg", border = NA, main = "Evacuation-time polygons")
plot(roads, add = TRUE, col = adjustcolor("grey20", alpha.f = 0.30), lwd = 0.35)
plot(escape_points_ltd, add = TRUE, pch = 22, bg = "red", col = "black", cex = 0.45)

17. Run the high-level wrapper with `run_evacpath()`

After the individual steps work, the same logic should run through run_evacpath(). The key tsunami-specific arguments are:

hazard_zone = zones$hazard_zone
escape_zone = zones$escape_zone
escape_roads_inset_x_m = 50
escape_roads_inset_y_m = 50
road_aware_escape_zone = TRUE

18. Final map check

This map uses land/water as background, the land-only evacuation-time polygons as the main layer, and the preferred escape points as squares.

evac_map <- result$time_grid
roads_plot <- crop(result$roads, evac_map)
escape_plot <- crop(result$escape_points, evac_map)

dem_bg <- result$dem
if (!same.crs(dem_bg, evac_map)) {
  dem_bg <- project(dem_bg, crs(evac_map))
}

e <- ext(evac_map)
xpad <- (xmax(e) - xmin(e)) * 0.06
ypad <- (ymax(e) - ymin(e)) * 0.06
e2 <- ext(xmin(e) - xpad, xmax(e) + xpad, ymin(e) - ypad, ymax(e) + ypad)

dem_bg <- crop(dem_bg, e2)
land_water <- ifel(dem_bg > params$land_threshold, 2, 1)
names(land_water) <- "land_water"

bg_cols <- c(
  adjustcolor("lightskyblue1", alpha.f = 0.65),
  adjustcolor("grey92", alpha.f = 1.00)
)
time_cols <- hcl.colors(100, "YlOrRd", rev = TRUE)

oldpar <- par(mar = c(4, 4, 3, 5), mgp = c(2.2, 0.7, 0), las = 1, bg = "white", xpd = FALSE)

plot(
  land_water,
  ext = e2,
  col = bg_cols,
  legend = FALSE,
  axes = TRUE,
  box = FALSE,
  main = "Modeled Pedestrian Evacuation Time",
  cex.main = 1.1
)

plot(
  evac_map,
  "EvacTimeAvg",
  add = TRUE,
  col = time_cols,
  border = NA,
  plg = list(title = "Minutes", cex = 0.75, title.cex = 0.85, x = "right")
)

plot(roads_plot, add = TRUE, col = adjustcolor("grey25", alpha.f = 0.30), lwd = 0.30)
plot(escape_plot, add = TRUE, pch = 22, bg = adjustcolor("red", alpha.f = 0.75), col = "black", cex = 0.38, lwd = 0.35)

legend(
  "topright",
  legend = c("Escape / safety points", "Roads", "Land", "Water"),
  pch = c(22, NA, 15, 15),
  pt.bg = c("red", NA, "grey92", "lightskyblue1"),
  col = c("black", adjustcolor("grey25", alpha.f = 0.5), "grey92", "lightskyblue1"),
  lty = c(NA, 1, NA, NA),
  lwd = c(NA, 1, NA, NA),
  pt.cex = c(0.8, NA, 1.1, 1.1),
  bty = "n",
  cex = 0.7,
  inset = 0.02
)
par(oldpar)

19. Write outputs with `write_evac_outputs()`

write_evac_outputs() writes the major spatial outputs from a completed result object. This should only run after the diagnostic workflow succeeds.

write_evac_outputs(
  result,
  output_dir = params$output_dir,
  prefix = "example"
)

list.files(params$output_dir, full.names = TRUE)

20. Troubleshooting checklist

If escape points look wrong, check these first:

Did you create both hazard_zone and escape_zone with prepare_tsunami_zones()?
Did you crop the roads used for escape detection with crop_roads_to_inner_extent()?
Did you combine buffered roads with the escape zone using make_road_aware_escape_zone()?
Did you use the road-aware boundary for find_escape_points()?
Does the DEM-derived land/water mask look correct?
If land/water are reversed, did you set dem_sign_multiplier = -1?
Are piers, ferry routes, or non-road features still present in the road layer?
Is the target CRS projected in meters?
Are you using dem_resolution to coarsen the conductance raster for testing, not only grid_resolution?

The expected tsunami pattern is:

zones <- prepare_tsunami_zones(inundation, dem, target_crs = target_crs)

roads_for_escape <- crop_roads_to_inner_extent(
  roads = roads,
  zone = zones$escape_zone,
  inset_x_m = 50,
  inset_y_m = 50
)

escape_boundary_zone <- make_road_aware_escape_zone(
  escape_zone = zones$escape_zone,
  roads = roads_for_escape,
  road_buffer_m = 2,
  crop_buffer_m = 3
)

escape_points <- find_escape_points(
  hazard_zone = escape_boundary_zone,
  roads = roads_for_escape
)

result <- run_evacpath(
  hazard_zone = zones$hazard_zone,
  escape_zone = zones$escape_zone,
  roads = roads,
  dem = zones$dem,
  target_crs = target_crs,
  escape_roads_inset_x_m = 50,
  escape_roads_inset_y_m = 50,
  road_aware_escape_zone = TRUE
)

References

Cordero, E., Ruiz Vélez, R., Huérfano Moreno, V., Sherman, C., 2025. Enhancing tsunami evacuation strategies in Puerto Rico using open-source least-cost path analysis. J. Disaster Sci. Manag. 1, 18. https://doi.org/10.1007/s44367-025-00018-y

Lewis, J., 2023. leastcostpath: Modelling Pathways and Movement Potential Within a Landscape.

Lewis, J., 2021. Probabilistic Modelling for Incorporating Uncertainty in Least Cost Path Results: a Postdictive Roman Road Case Study. Journal of Archaeological Method and Theory 28, 911-924. https://doi.org/10.1007/s10816-021-09522-w

Evacuation workflow: function-by-function diagnostics

1. Setup

2. Locate example input files

3. Read the spatial inputs with read_spatial()

3.1 Raw input diagnostics

4. Clean roads with clean_roads()

5. Prepare tsunami zones with prepare_tsunami_zones()

5.1 Check the land/water classification

5.2 Check the two tsunami zones

6. Crop roads to an inner extent before escape-point detection

7. Combine buffered roads with the escape zone

8. Compare escape-point methods

9. Prepare core inputs with prepare_evac_inputs()

10. Create the evacuation grid with make_evac_grid()

11. Create the road mask with make_road_mask()

12. Create road origin points with make_road_origins()

13. Create the conductance surface with make_conductance_surface()

14. Calculate least-cost distance with calc_min_distance_to_safety()

15. Convert distance to time with calc_evac_time()

16. Create evacuation polygons with make_evac_polygons()

17. Run the high-level wrapper with run_evacpath()