params <-
list(run_diagnostics = FALSE, run_lcp = FALSE, data_dir = "data-raw", 
    target_crs = "EPSG:32748", region_name = "Example region", 
    dem_sign_multiplier = 1L, inundation_threshold = 0L, land_threshold = 0L, 
    water_threshold = 0L, resample_method = "bilinear", road_exclude_field = "man_made", 
    road_exclude_values = "pier", grid_multiplier = 50L, dem_multiplier = 50L, 
    road_buffer_m = 2L, escape_buffer_m = 5L, final_road_buffer_m = 3L, 
    escape_roads_inset_x_m = 50L, escape_roads_inset_y_m = 50L, 
    road_aware_escape_zone = TRUE, walking_speed_mps = 1.22, 
    seed = 23401L, output_dir = "_outputs/example")

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 8,
  fig.height = 6,
  dpi = 120
)

## ----render-command, eval = FALSE---------------------------------------------
# rmarkdown::render(
#   "vignettes/diagnostic-example.Rmd",
#   params = list(
#     run_diagnostics = TRUE,
#     run_lcp = TRUE,
#     data_dir = "data-raw",
#     target_crs = "EPSG:32748",
#     dem_sign_multiplier = 1,
#     grid_multiplier = 50,
#     dem_multiplier = 50,
#     escape_roads_inset_x_m = 50,
#     escape_roads_inset_y_m = 50,
#     road_aware_escape_zone = TRUE
#   )
# )

## ----render-fast, eval = FALSE------------------------------------------------
# rmarkdown::render(
#   "vignettes/diagnostic-example.Rmd",
#   params = list(
#     run_diagnostics = TRUE,
#     run_lcp = FALSE,
#     data_dir = "data-raw"
#   )
# )

## ----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)

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_
}

## ----switches-----------------------------------------------------------------
run_diagnostics
run_lcp

## ----locate-files-------------------------------------------------------------
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

example_data_available <- all(input_files$exists)
example_data_available

## ----read-inputs, eval = run_diagnostics && example_data_available------------
# 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)
# )

## ----raw-diagnostics, eval = run_diagnostics && example_data_available--------
# 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))

## ----plot-raw-inputs, eval = run_diagnostics && example_data_available, fig.height = 7----
# 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)

## ----clean-roads, eval = run_diagnostics && example_data_available------------
# 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)

## ----prepare-tsunami-zones, eval = run_diagnostics && example_data_available----
# 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")

## ----land-water-check, eval = run_diagnostics && example_data_available, fig.height = 7----
# 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)

## ----plot-tsunami-zones, eval = run_diagnostics && example_data_available, fig.height = 7----
# 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)

## ----crop-roads-inner-extent, eval = run_diagnostics && example_data_available----
# 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-roads-inner-extent, eval = run_diagnostics && example_data_available, fig.height = 7----
# 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
# )

## ----road-aware-escape-zone, eval = run_diagnostics && example_data_available----
# 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-road-aware-zone, eval = run_diagnostics && example_data_available, fig.height = 7----
# 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)

## ----compare-escape-points, eval = run_diagnostics && example_data_available----
# 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.")

## ----plot-escape-comparison, eval = run_diagnostics && example_data_available, fig.height = 8----
# 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)

## ----assign-core-objects, eval = run_diagnostics && example_data_available----
# hazard_zone <- zones$hazard_zone
# escape_zone <- zones$escape_zone
# dem <- zones$dem
# roads <- roads_clean

## ----prepare-core-inputs, eval = run_diagnostics && example_data_available----
# 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.")

## ----make-grid, eval = run_diagnostics && example_data_available--------------
# 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-grid, eval = run_diagnostics && example_data_available--------------
# 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)

## ----make-road-mask, eval = run_diagnostics && example_data_available---------
# 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-road-mask, eval = run_diagnostics && example_data_available---------
# 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)

## ----make-road-origins, eval = run_diagnostics && example_data_available------
# 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-road-origins, eval = run_diagnostics && example_data_available------
# 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)

## ----make-conductance, eval = run_diagnostics && example_data_available-------
# 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-conductance, eval = run_diagnostics && example_data_available-------
# 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)

## ----sample-destinations, eval = run_diagnostics && example_data_available----
# 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")

## ----calc-distance, eval = run_diagnostics && example_data_available && run_lcp----
# 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-distance-points, eval = run_diagnostics && example_data_available && run_lcp----
# 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)

## ----calc-time, eval = run_diagnostics && example_data_available && run_lcp----
# 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"
#   )
# )

## ----make-polygons, eval = run_diagnostics && example_data_available && run_lcp----
# 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, eval = run_diagnostics && example_data_available && run_lcp, fig.height = 7----
# 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)

## ----run-wrapper, eval = run_diagnostics && example_data_available && run_lcp, message=FALSE,echo=FALSE, eval = FALSE----
# result <- run_evacpath(
#   hazard_zone = zones$hazard_zone,
#   escape_zone = zones$escape_zone,
#   roads = roads_clean,
#   dem = zones$dem,
#   target_crs = params$target_crs,
#   region_name = params$region_name,
#   grid_resolution = terra::res(zones$dem) * params$grid_multiplier,
#   dem_resolution = terra::res(zones$dem) * params$dem_multiplier,
#   road_buffer_m = params$road_buffer_m,
#   escape_buffer_m = params$escape_buffer_m,
#   final_road_buffer_m = params$final_road_buffer_m,
#   escape_roads_inset_x_m = params$escape_roads_inset_x_m,
#   escape_roads_inset_y_m = params$escape_roads_inset_y_m,
#   road_aware_escape_zone = params$road_aware_escape_zone,
#   escape_zone_road_buffer_m = params$road_buffer_m,
#   escape_zone_crop_buffer_m = params$final_road_buffer_m,
#   max_origins = params$max_origins,
#   max_destinations = params$max_destinations,
#   seed = params$seed,
#   walking_speed_mps = params$walking_speed_mps,
#   clip_mode = "hazard",
#   progress = TRUE,
#   progress_every = 500,
#   lcp_check_locations = FALSE
# )
# 
# result
# 
# assert_true("escape_boundary_zone" %in% names(result), "Result does not contain escape_boundary_zone.")
# assert_true("roads_for_escape" %in% names(result), "Result does not contain roads_for_escape.")
# assert_true("time_grid" %in% names(result), "Result does not contain time_grid.")
# assert_true("EvacTimeAvg" %in% names(result$time_grid), "Result time_grid lacks EvacTimeAvg.")

## ----final-map, eval = run_diagnostics && example_data_available && run_lcp, fig.width = 9, fig.height = 7, eval = FALSE----
# 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)

## ----write-outputs, eval = run_diagnostics && example_data_available && run_lcp, eval = FALSE----
# write_evac_outputs(
#   result,
#   output_dir = params$output_dir,
#   prefix = "example"
# )
# 
# list.files(params$output_dir, full.names = TRUE)