## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4)

## ----setup, message = FALSE---------------------------------------------------
library(blueterra)
library(terra)

## ----read---------------------------------------------------------------------
hitw <- read_bathy(blueterra_example("hitw"))
hoyo <- read_bathy(blueterra_example("hoyo"))
slope <- read_bathy(blueterra_example("slope"))
rectangles <- terra::vect(blueterra_example("sampling_rectangles"))

hitw_rect <- rectangles[rectangles$site_id == "hitw", ]
hoyo_rect <- rectangles[rectangles$site_id == "hoyo", ]

bathy_info(hitw)
rectangles[, c("site_id", "site_name", "feature_type")]

## ----prepare------------------------------------------------------------------
prepared <- prepare_bathy(
  hitw,
  depth_range = c(-220, -25),
  smooth = TRUE,
  smooth_window = 3
)

check_bathy_units(prepared, units = "m", positive_depth = FALSE)
range(terra::values(prepared), na.rm = TRUE)

## ----map, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Hole-in-the-Wall bathymetry with hillshade, contours, and sampling rectangle."----
plot_bathy(
  prepared,
  contours = TRUE,
  contour_interval = 25,
  vectors = hitw_rect,
  title = "Hole-in-the-Wall Bathymetry",
  subtitle = "Hillshade, contours, and sampling rectangle"
)

## ----metrics------------------------------------------------------------------
terrain <- derive_terrain(
  prepared,
  metrics = c(
    "slope", "aspect", "northness", "eastness", "tri", "rugosity",
    "bpi", "curvature", "surface_area_ratio"
  )
)

names(terrain)
assign_process_groups(terrain)
summarize_process_groups(terrain)

## ----metric-stack, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Metric stack showing slope, TRI, BPI, and curvature."----
plot_metric_stack(terrain[[c("slope_deg", "tri", "bpi_3x3", "curvature")]])

## ----zones--------------------------------------------------------------------
zone_summary <- summarize_terrain(
  terrain,
  hitw_rect,
  fun = c("mean", "sd", "min", "max")
)

zone_summary[, c("site_id", "site_name", "slope_deg_mean", "bpi_3x3_mean")]

## ----transects----------------------------------------------------------------
orientation <- estimate_surface_orientation(prepared, hitw_rect)
transects <- make_transects(hitw_rect, spacing = 75, bathy = prepared)
cross_sections <- sample_transects(transects, prepared, n = 12)

orientation
unique(as.data.frame(transects)[, c("angle_deg", "angle_source", "mean_aspect_deg")])
head(cross_sections[, c("transect_id", "distance", "bathy_m")])

## ----transect-map, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Terrain-oriented transects over hillshaded bathymetry."----
plot_transects(
  prepared,
  transects,
  color_by = "transect_id",
  show_legend = FALSE,
  contour_interval = 25,
  title = "Terrain-Oriented Transects"
)

## ----cross-section-plot, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Bathymetric cross-section profiles with bathy_m on the y-axis."----
plot_cross_sections(
  cross_sections,
  value_col = "bathy_m",
  show_legend = TRUE,
  mean_profile = TRUE,
  mean_profile_na_rm = TRUE,
  normalize_distance = FALSE,
  profile_direction = "top_to_bottom",
  title = "Bathymetric Cross-Sections",
  subtitle = "Profiles read from shallow to deep terrain"
)

## ----depth-profile, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Single bathymetric profile oriented from shallow terrain toward deeper terrain."----
one_transect <- cross_sections[
  cross_sections$transect_id == cross_sections$transect_id[1],
]

plot_depth_profile(
  one_transect,
  value_col = "bathy_m",
  profile_direction = "top_to_bottom",
  title = "Bathymetry Along One Transect",
  subtitle = "Distance is oriented from shallow to deep terrain"
)

## ----metric-profile, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Slope profile along one transect."----
metric_samples <- sample_transects(
  transects,
  c(prepared, terrain[["slope_deg"]]),
  n = 25
)
metric_one <- metric_samples[
  metric_samples$transect_id == metric_samples$transect_id[1],
]

plot_depth_profile(
  metric_one,
  depth_col = "bathy_m",
  value_col = "slope_deg",
  profile_direction = "top_to_bottom",
  title = "Slope Along Depth",
  subtitle = "Bathymetry is on the y-axis; slope is on the x-axis"
)

## ----isobaths-----------------------------------------------------------------
isobaths <- extract_isobaths(prepared, depths = c(-50, -80, -120))
corridors <- make_isobath_corridors(prepared, depths = c(-50, -80, -120), width = 5)

corridors[, c("contour_value", "depth_label", "corridor_id")]
summarize_isobath_terrain(terrain, corridors)[, c("contour_value", "slope_deg_mean", "bpi_3x3_mean")]

## ----corridor-map, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Isobath corridors over hillshaded bathymetry with source isobaths in black."----
plot_isobath_corridors(
  corridors,
  prepared,
  isobaths = isobaths,
  background_contours = FALSE,
  title = "Isobath Corridors and Source Isobaths",
  subtitle = "Corridors use a 5 m buffer around each source isobath"
)

## ----model-ready--------------------------------------------------------------
hoyo_prepared <- prepare_bathy(hoyo, depth_range = c(-220, -25), smooth = TRUE)
hoyo_metrics <- derive_terrain(hoyo_prepared, metrics = c("slope", "tri", "bpi", "curvature"))

hitw_cells <- sample_terrain_cells(
  terrain[[c("slope_deg", "tri", "bpi_3x3", "curvature")]],
  size = 45,
  method = "regular"
)
hitw_cells$site <- "Hole-in-the-Wall"

hoyo_cells <- sample_terrain_cells(
  hoyo_metrics[[c("slope_deg", "tri", "bpi_3x3", "curvature")]],
  size = 45,
  method = "regular"
)
hoyo_cells$site <- "El Hoyo"

comparison <- rbind(hitw_cells, hoyo_cells)
pca <- terrain_pca(
  comparison,
  vars = c("slope_deg", "tri", "bpi_3x3", "curvature")
)

pca$variance
pca_axis_labels(pca)

## ----pca-plot, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Terrain PCA with site-colored points and loading labels in the axis text."----
plot_process_pca(
  pca,
  color_col = "site",
  title = "Terrain PCA"
)

## ----model-matrix-------------------------------------------------------------
model_matrix <- prepare_model_matrix(
  comparison,
  response = "site",
  vars = c("slope_deg", "tri", "bpi_3x3", "curvature")
)

head(model_matrix)
terrain_correlation(comparison[, c("slope_deg", "tri", "bpi_3x3", "curvature")])

## ----custom-metric------------------------------------------------------------
slope_tri <- derive_custom_metric(
  terrain,
  name = "slope_tri_index",
  expression = quote(slope_deg * tri)
)

extended_terrain <- add_metric_layers(terrain, slope_tri)
names(extended_terrain)

## ----custom-metric-map, eval=requireNamespace("ggplot2", quietly = TRUE), fig.alt="Custom slope-TRI metric over hillshaded bathymetry."----
plot_metric(
  extended_terrain,
  metric = "slope_tri_index",
  bathy = prepared,
  contours = TRUE,
  contour_interval = 25,
  title = "Custom Slope-TRI Index"
)