---
title: "Get started with blueterra"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Get started with blueterra}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

```{r setup, message = FALSE}
library(blueterra)
library(terra)
```

`blueterra` works from bathymetric and elevation rasters already available to
the analyst. The package keeps the workflow close to `terra`: rasters remain
`SpatRaster` objects, sampling zones and transects remain `SpatVector` objects,
and extracted values become ordinary tables for summaries and models.

<img src="../man/figures/study-area-pr-southwest-shelf-margin.png" width="100%" alt="Study-area context for the example data along the southwest Puerto Rico shelf margin." />

The examples use reduced bathymetry and sampling rectangles from
Hole-in-the-Wall, El Hoyo, and a broader slope clip along the southwest Puerto
Rico shelf margin near La Parguera, Puerto Rico. The files are compact, but
they retain real relief, depth gradients, slope breaks, and sampling geometry.

## Read Example Data

```{r 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")]
```

`bathy_info()` is a first check on raster dimensions, value range, resolution,
extent, and CRS. These fields affect slope, focal-window metrics, buffers, and
distance along transects.

## Prepare the Bathymetry

The example rasters store bathymetry as negative elevation. `prepare_bathy()`
preserves that convention unless a conversion is requested explicitly.

```{r 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)
```

```{r 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"
)
```

Hillshade is used here as visual relief. It helps the reader see escarpments
and local relief, but it is not a predictor unless the analyst explicitly
derives and includes it.

## Derive Metrics and Process Groups

Terrain metrics describe different aspects of the same bathymetric surface.
Slope and aspect describe local orientation; northness and eastness convert
aspect to linear components; TRI and rugosity describe local relief
variability; BPI describes relative position; curvature summarizes local
surface bending.

```{r 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)
```

```{r 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")]])
```

Process groups keep related derivatives together. They are interpretation
categories for terrain form, not direct measurements of currents, sediment
transport, habitat condition, or ecological response.

## Summarize Sampling Rectangles

Sampling rectangles provide a compact example of zone-based extraction. Each
polygon is treated as a spatial sampling frame, and summary statistics are
calculated from raster cells inside that frame.

```{r 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")]
```

Set `exact = TRUE` for boundary-sensitive polygon summaries when the optional
`exactextractr` dependency is available. Exact raster--polygon intersections
then supply coverage fractions used to weight means, population standard
deviations, medians, sums, and effective cell counts; minima and maxima use
positively intersected cells.

## Transects and Cross-Sections

Transects convert the raster surface into profiles. When `bathy` is supplied,
`make_transects()` can estimate a terrain-oriented line angle from local aspect
instead of requiring the analyst to choose a fixed direction by hand.

```{r 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",
  "orientation_resultant_length"
)])
head(cross_sections[, c("transect_id", "distance", "bathy_m")])
```

Automatic orientation metadata include `orientation_resultant_length`. Values
near one indicate a concentrated aspect direction; values near zero signal
cancelling aspect vectors and an unreliable mean angle. In that case, use a
manual or bounding-box orientation rather than overinterpreting the automatic
direction.

```{r 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"
)
```

```{r 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"
)
```

The y-axis is explicitly set to `bathy_m`. This prevents transect metadata such
as width, angle, or offset from being mistaken for the raster value. Distance
is oriented from the top or shallow endpoint toward the bottom or deeper
endpoint, and the plotted distance is reset to zero after trimming empty
profile ends. The mean cross-section is averaged from interpolated transects
with missing interpolated values removed by default, so longer profiles still
contribute beyond the shortest transect.

```{r 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 profiles can use the same transect geometry. When both a bathymetry
column and a metric column are present, `plot_depth_profile()` places the
metric on the x-axis and bathymetry or depth on the y-axis. This is useful for
showing how terrain attributes change along the depth gradient.

```{r 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"
)
```

## Isobath Corridors

Isobath corridors summarize terrain along depth horizons. The source isobaths
are shown in black so the reader can see which contour each corridor buffers.
Here `width = 5` is a one-sided buffer distance, producing a nominal 10 m
full-width corridor.

```{r 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", "buffer_distance",
  "nominal_corridor_width", "overlap_policy"
)]
summarize_isobath_terrain(terrain, corridors)[, c("contour_value", "slope_deg_mean", "bpi_3x3_mean")]
```

```{r 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 = "5 m is the one-sided buffer distance (10 m nominal full width)"
)
```

Corridors are independent buffers and may overlap. A cell in an overlap can be
included in more than one corridor summary, so corridor summaries are not
mutually exclusive or additive.

## PCA and Model-Ready Tables

Cell samples and extracted summaries can be sent directly into exploratory
models. PCA is useful for checking whether terrain metrics form separable
gradients across sites or sampling frames.

```{r 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)
```

```{r 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"
)
```

The same sampled table can be converted to a model matrix after selecting the
terrain variables that belong in the analysis design.

```{r 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 Metrics

Project-specific metrics can be added when they share the same grid geometry
as the terrain stack. The expression below combines local gradient and relief
variability into a compact demonstration index.

```{r 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)
```

```{r 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"
)
```