This vignette walks through a complete oystermapR workflow using a simulated survey of Example Bay — a fictional sheltered coastal inlet with realistic physical characteristics. The full pipeline is covered:
The example data files (example_bay_adcp.csv,
example_bay_soundings.xyz, and
example_bay_ctd.csv) are included in the package and
represent a simulated mid-summer survey.
Example Bay is a sheltered inlet approximately 6 km × 4 km in extent, with depths ranging from 2 m at the shoreline to around 22 m in the central channel. A six-hour ADCP transect was run on a flood tide during June. Bathymetric soundings were collected concurrently using a single-beam echosounder.
The target species is Ostrea edulis, the European Flat Oyster, whose optimal conditions include: depths of 1–20 m, current speeds of 0.05–0.4 m/s, temperatures of 5–25°C, and salinities of 25–35 PSU.
read_nortek_adcp() reads the Nortek Signature 500 merged
CSV export. It auto-detects velocity bins, spatially averages ensembles
onto a grid, and derives bed shear stress from the near-bed velocity
profile. The spatial_res argument controls the decimal
places used for lat/lon grid binning; 2 gives approximately
1 km cells, suitable for a bay-scale survey.
adcp_file <- system.file("extdata", "example_bay_adcp.csv", package = "oystermapR")
adcp <- read_nortek_adcp(
file = adcp_file,
spatial_res = 2,
verbose = TRUE
)
head(adcp[, c("lat", "lon", "current_velocity", "shear_stress")])The output contains one row per grid cell.
current_velocity is the mean near-bed flow magnitude;
shear_stress is the estimated bed shear stress in N/m²,
used by predict_oyster() for suspension-feeder scoring.
read_soundings_xyz() reads a space-delimited XYZ point
cloud, grids and averages the depths, and derives slope and rugosity via
finite differences.
xyz_file <- system.file("extdata", "example_bay_soundings.xyz", package = "oystermapR")
bathy <- read_soundings_xyz(
file = xyz_file,
spatial_res = 2,
min_soundings = 5,
verbose = TRUE
)
head(bathy[, c("lat", "lon", "depth", "slope", "roughness")])slope is the maximum downslope gradient in degrees at
each grid cell; values above ~15° typically exclude oyster settlement in
the scoring model. roughness is a dimensionless rugosity
index (1.0 = flat; higher = more complex substrate).
The ADCP and echosounder provide hydrodynamic and bathymetric
variables, but predict_oyster() can also use water-quality
data when available. The Example Bay survey included a 5 × 4 grid of CTD
casts recording temperature, salinity, chlorophyll_a,
pH, and total alkalinity, stored as a
plain CSV.
read_generic_csv() handles any tabular sensor export
with flexible column matching. Column names are recognised automatically
when they follow standard conventions
(lat/lon, temperature,
salinity, chlorophyll_a, ph,
alkalinity). For instruments that export non-standard
headers, supply an explicit col_map:
# Example with non-standard headers from a YSI EXO sonde export:
ctd <- read_generic_csv(
"ysi_export.csv",
col_map = c(lat = "GPS_Lat", lon = "GPS_Lon",
temperature = "Temp_C", salinity = "Sal_PSU",
ph = "pH_total", alkalinity = "TA_umolkg")
)For Example Bay the column names match the oystermapR conventions directly:
ctd_file <- system.file("extdata", "example_bay_ctd.csv", package = "oystermapR")
ctd <- read_generic_csv(
file = ctd_file,
spatial_res = 2,
verbose = TRUE
)
head(ctd[, c("lat", "lon", "temperature", "salinity", "chlorophyll_a",
"ph", "alkalinity")])Spatially averaging at spatial_res = 2 collapses the
four replicate casts per station to a single value per ~1 km grid cell.
The ph and alkalinity columns will be used in
Step 9 for aragonite saturation scoring.
merge_sensor_data() performs a full outer join of any
number of sensor dataframes on rounded lat/lon grid keys. Cells present
in only one source are retained with NA for the missing
variables, so no data is discarded. All three layers — ADCP, bathymetry,
and CTD — are combined in a single call:
survey <- merge_sensor_data(adcp = adcp, bathy = bathy, ctd = ctd)
cat("Merged survey:", nrow(survey), "grid cells\n")
cat("Columns:", paste(names(survey), collapse = ", "), "\n")A substrate_hardness column would typically come from a
sidescan mosaic processed via read_sonar_tif(). Here we add
a representative value for a mixed shell-gravel substrate, typical of
productive inshore oyster ground.
set.seed(101)
n <- nrow(survey)
# 0 = soft mud, 1 = hard rock; 0.3--0.7 = shell/gravel mix
survey$substrate_hardness <- round(runif(n, 0.30, 0.70), 2)qc_survey_data() applies three complementary checks to
every numeric column:
"range_fail".iqr_k × IQR from the median are flagged
"outlier" (default k = 3, Tukey’s outer fence).Setting apply_flags = TRUE replaces flagged values with
NA so they are silently skipped downstream rather than
causing erroneous scores.
survey_clean <- qc_survey_data(
df = survey,
apply_flags = TRUE,
verbose = TRUE
)
# Count any flags raised across all columns
flag_cols <- grep("^qc_flag_", names(survey_clean), value = TRUE)
n_flagged <- sum(sapply(survey_clean[flag_cols], function(x) sum(!is.na(x) & x != "pass")))
cat("Total flagged values replaced with NA:", n_flagged, "\n")The QC step is non-destructive by default
(apply_flags = FALSE) — it adds
qc_flag_<variable> columns so you can inspect which
cells were problematic before deciding whether to exclude them.
predict_oyster() applies AHP-weighted scoring rules to
each available variable, combines them into a suitability score in [0,
1], and classifies locations as High / Moderate / Low / Very Low /
Excluded.
result <- predict_oyster(
data = survey_clean,
species = "ostrea_edulis",
verbose = TRUE
)
# Summary of suitability classes
table(result$suitability_class)# Mean score and range
cat(sprintf(
"Suitability: mean = %.2f, range = %.2f -- %.2f\n",
mean(result$suitability, na.rm = TRUE),
min(result$suitability, na.rm = TRUE),
max(result$suitability, na.rm = TRUE)
))The result dataframe retains all input columns plus
suitability, suitability_class,
data_completeness (fraction of variables scored),
n_layers_scored (integer count of variables contributing at
each point), and per-variable component scores
(score_depth, score_current_velocity,
score_ph, score_omega_aragonite, etc.).
# Points with fewer scored variables may have less reliable scores
summary(result$n_layers_scored)
table(result$n_layers_scored)oystermapR includes optional risk modules that can be appended to the result. Here we add wave exposure (derived from the current speed data and fetch geometry) and a simple HAB risk score.
# Wave exposure: uses current_velocity and depth as proxies for fetch exposure
result <- score_wave_exposure(result, verbose = FALSE)
# HAB risk: without live ICES data, scores from chlorophyll_a alone
result <- score_hab_risk(result, verbose = FALSE)
cat("Wave exposure range:", round(range(result$wave_exposure, na.rm=TRUE), 3), "\n")
cat("HAB risk range: ", round(range(result$hab_risk, na.rm=TRUE), 3), "\n")export_geotiff() interpolates the suitability scores
onto a regular raster and writes a five-band GeoTIFF.
export_qml_style() writes a matching QGIS colour-ramp style
file (.qml) so the layer renders immediately with the
standard oystermapR colour scheme (red = excluded, green = high
suitability).
# Write five-band GeoTIFF and companion QGIS style file
export_geotiff(
df = result,
path = "example_bay_suitability.tif",
resolution = 0.001,
contours = TRUE
)
export_qml_style("example_bay_suitability.tif")Load example_bay_suitability.tif into QGIS via
Layer → Add Layer → Add Raster Layer, then right-click
the layer and choose Load Layer Style to apply the
.qml file.
The GeoTIFF contains five bands:
| Band | Name | Description |
|---|---|---|
| 1 | suitability |
Continuous score [0, 1] |
| 2 | excluded_mask |
1 = hard-excluded by a threshold |
| 3 | n_observations |
Survey points per raster cell |
| 4 | dist_to_nearest_m |
Distance to nearest survey point (m) |
| 5 | n_layers_scored |
Number of variables contributing to the score |
Band 5 is particularly useful as a data-coverage overlay: cells where only 2–3 variables were scored are visually distinguishable from cells with full data.
The result dataframe includes a score_<variable>
column for every variable that was scored. Comparing these helps
diagnose which environmental factor is the main limiting constraint at a
site.
score_cols <- grep("^score_", names(result), value = TRUE)
# Mean component score per variable (higher = more suitable)
col_means <- sort(colMeans(result[score_cols], na.rm = TRUE))
print(round(col_means, 3))Variables scoring consistently below 0.5 are the main limiting factors for Ostrea edulis at this site. Scores near 1.0 indicate that variable is not constraining growth.
pH and aragonite saturation state (Ω_arag) are scored for all 17
species. If ph and alkalinity are present in
the merged survey data, predict_oyster() computes
omega_aragonite automatically via
calculate_aragonite() before the scoring step — no manual
pre-processing is needed.
# Verify aragonite was auto-calculated and scored
"omega_aragonite" %in% names(result) # column present
"score_ph" %in% names(result) # pH scored
"score_omega_aragonite" %in% names(result) # aragonite scored
# Distribution of omega_aragonite across Example Bay
summary(result$omega_aragonite)If your sensor does not record alkalinity, you can estimate Ω_arag from typical open-ocean values for your region and supply it as a pre-computed column, or omit it — the scoring model will redistribute its weight across remaining variables.
# Manual calculation: sensors gave pH only, alkalinity approximated
df$alkalinity <- 2300 # µmol/kg — representative NE Atlantic value
df$omega_aragonite <- calculate_aragonite(
pH = df$ph,
alkalinity = df$alkalinity,
temperature = df$temperature,
salinity = df$salinity
)variable_impact() summarises the contribution of each
environmental variable to the suitability score across the dataset. It
is the primary QA tool for understanding model behaviour and planning
future surveys.
The output columns are:
norm_weight_pct — the variable’s share
of the total AHP weight at a typical fully-observed point.mean_score — average per-variable
score across non-excluded points. Values below ~0.5 identify
environmental bottlenecks at this site.mean_contribution —
norm_weight_pct × mean_score / 100. The net weighted
contribution to the suitability score. Highest-impact variables have
both high weight and high score.pct_coverage — percentage of
non-excluded points where this variable was observed (non-NA). Low
values indicate data sparsity.# Variables scoring below 0.5 are potential site limiters
impact[impact$mean_score < 0.5, c("variable", "norm_weight_pct", "mean_score")]
# Variables with sparse data coverage
impact[impact$pct_coverage < 80, c("variable", "pct_coverage")]Use sort_by = "pct_coverage" to prioritise sensor
deployment for the next survey, or sort_by = "mean_score"
to focus on the worst-performing variables.
area_summary() converts the point-based result into
habitat area estimates at sub-hectare resolution. This is especially
useful for restoration planning, where the difference between 200 m² and
800 m² of suitable habitat is operationally significant.
# Auto-estimate cell size from median nearest-neighbour spacing
s <- area_summary(result, verbose = TRUE)The function prints a concise summary to the console. The returned list contains three elements:
# Per-class breakdown
s$class_summary[, c("class", "area_m2", "area_ha", "pct_total_area",
"mean_suitability")]
# Totals
s$total[c("surveyed_area_m2", "suitable_area_m2", "pct_suitable", "cell_size_m")]
# Contiguous patches of High + Moderate suitability
head(s$patches)
# Patches meeting the OSPAR 100 m² viable area threshold
viable <- s$patches[s$patches$viable, ]
cat(nrow(viable), "viable patches; largest:", round(max(viable$area_m2)), "m²\n")For ROV or AUV surveys with a known fixed grid resolution, supply it explicitly to avoid auto-estimation error:
# 5 m ROV grid survey
s5m <- area_summary(result, cell_size_m = 5, viable_area_m2 = 100)
# 25 m ADCP trackline survey, larger minimum viable unit for production scale
s25m <- area_summary(result, cell_size_m = 25, viable_area_m2 = 500)plot_tolerance() draws the mathematical scoring function
for any scored variable directly from the species tolerance parameters —
no dataset required. This is useful for QA (confirming threshold values
match published literature), stakeholder reporting (“here is exactly
what the model rewards”), and creating publication-quality figures.
# All four seasons overlaid on a single plot
plot_tolerance("ostrea_edulis", "temperature", season = "all")# Ocean acidification variables
plot_tolerance("ostrea_edulis", "ph")
plot_tolerance("ostrea_edulis", "omega_aragonite")# Dissolved oxygen: compare tolerance between species
plot_tolerance("ostrea_edulis", "dissolved_oxygen")
plot_tolerance("crassostrea_iredalei", "dissolved_oxygen")Curves are colour-coded: green background = optimal zone, orange =
acceptable/poor, red = hard-excluded. Dashed red vertical lines show the
exclusion thresholds. The function returns a ggplot2 object
invisibly, so it can be saved or modified:
p <- plot_tolerance("ostrea_edulis", "salinity")
ggplot2::ggsave("salinity_tolerance_O_edulis.png", p, width = 8, height = 5)For data-driven response curves showing how the model actually
responded to observed covariate values in this survey, use
sensitivity_analysis() instead.
validate_against_records() to calculate AUC, TSS, F1, and
Brier score against model predictions.update_species_tolerances() can refine the O.
edulis scoring parameters from field observations at this
site."magallana_gigas" and call
compare_species() to see which species is better suited to
each grid cell. Useful for restoration target selection.sensitivity_analysis() produces partial dependence curves
showing how suitability varies with each variable at the observed
covariate values in this survey (data-driven, contrast to
plot_tolerance() which is parameter-driven).fetch_live_environmental_data() can replace the simulated
temperature and salinity with real-time CMEMS model output or ICES
observational data; see ?oystermapR_live_config for
credential setup.