Predict and Map Oyster Growth Suitability from Environmental Data

Predicts spatial suitability for oyster growth from environmental survey data using rank-reciprocal weighted scoring. Users supply sensor data from Acoustic Doppler Current Profilers (ADCP), Conductivity-Temperature-Depth (CTD) sensors, bathymetric sonar, and sidescan sonar, specify a target species, and receive per-location suitability scores, a five-band 'GeoTIFF' heatmap for 'QGIS', contour lines, and a formatted PDF or HTML report. Supports seventeen species across global aquaculture regions, including Ostrea edulis, Magallana gigas, Crassostrea virginica, Crassostrea hongkongensis, and thirteen further species; see list_species(). Includes ocean acidification scoring via in-house aragonite saturation state (Omega_arag) calculation using Lueker et al. (2000) and Mucci (1983) equilibrium constants (no external dependencies), variable impact diagnostics (variable_impact()), fine-scale habitat area analysis for restoration reporting in m2 with contiguous patch identification (area_summary()), tolerance curve visualisation (plot_tolerance()), season-aware scoring, tidal height correction, Bayesian tolerance parameter updating from field observations, spatial block cross-validation (Roberts et al., 2017, ), permutation variable importance, wave exposure and sediment stability modules, Harmful Algal Bloom (HAB) risk and anthropogenic disturbance scoring with optional live International Council for the Exploration of the Sea (ICES) data integration, hybrid larval dispersal connectivity scoring (union-find Gaussian kernel plus optional 'OpenDrift' or Finite Volume Community Ocean Model ('FVCOM') connectivity matrix), batch multi-species comparison, and evidence-weighted scoring in which each variable's weight is its published rank multiplied by the measured reliability and realised spread of that layer, with the fraction of the species' evidence base the survey observed reported alongside the score rather than normalised away, and lineage-specific tolerance profiles for related populations, in which thinly sampled lineages borrow strength from relatives across a genetic distance matrix only where tolerance is shown to track relatedness, with inverse mapping of where each lineage is more likely to occur and which source population suits a site.

oystermapR

Predict and map oyster growth suitability from environmental survey data

oystermapR takes tabular sensor data — from ADCPs, CTDs, bathymetric sonar, or standard CSV files — applies species-specific rank-reciprocal weighted scoring rules (each variable's weight is the inverse of its importance rank), and returns per-location suitability scores alongside a five-band GeoTIFF heatmap ready to load in QGIS.

Supported species

Installation

# From CRAN (once accepted)
install.packages("oystermapR")

# Development version from GitHub
devtools::install_github("trissyboats/oystermapR")

Requires R >= 4.1.0. Core dependencies: dplyr, terra, sf, cli, rlang.

Quick start

library(oystermapR)

# Load your survey data
df <- read.csv("my_survey.csv")

# Run the prediction
result <- predict_oyster(
  data           = df,
  species        = "ostrea_edulis",
  output_geotiff = "oyster_suitability.tif",
  verbose        = TRUE
)

# Inspect high-suitability sites
subset(result, suitability_class == "High")

# Export matching QGIS colour ramp style
export_qml_style("oyster_suitability.tif")

The output dataframe contains suitability (0–1), suitability_class (High / Moderate / Low / Very Low / Excluded), data_completeness (fraction of variables scored), n_layers_scored (integer count of variables that contributed at each point), and per-variable component scores (score_depth, score_temperature, etc.).

Input data

Your CSV needs at minimum lat, lon, and date. Any environmental columns present are scored automatically — missing variables are skipped and their weights redistributed. Column names are matched case-insensitively against the aliases below.

A sample dataset is included at inst/extdata/sample_survey.csv.

Ocean acidification scoring

pH and aragonite saturation state (Ω_arag) are scored for all 17 species. If your data contains ph and alkalinity columns, predict_oyster() automatically computes omega_aragonite before scoring — no preprocessing required.

# Aragonite auto-calculated from pH + alkalinity columns
result <- predict_oyster(df, "ostrea_edulis")

# Manual calculation (no external packages)
df$omega_aragonite <- calculate_aragonite(
  pH          = df$ph,
  alkalinity  = df$alkalinity,    # µmol/kg
  temperature = df$temperature,   # °C
  salinity    = df$salinity       # PSU
)

calculate_aragonite() implements the full carbonate chemistry using Lueker et al. (2000) K1/K2, Dickson (1990) KB, Millero (1995) KW, Mucci (1983) Ksp_arag, Uppstrom (1974) [B_T], and Riley & Tongudai (1967) [Ca²⁺]. No external package dependencies.

Diagnostics and QA

Variable impact

variable_impact() returns a ranked table showing which environmental variables contribute most to the suitability score across a dataset. Use it to identify limiting factors and data gaps.

impact <- variable_impact(result, "ostrea_edulis")
#>   variable          rank  norm_weight_pct  mean_score  mean_contribution  n_points  pct_coverage
#>   temperature          1            24.1        0.82               0.198       842          99.5
#>   depth                2            18.7        0.91               0.170       842         100.0
#>   dissolved_oxygen     5            12.4        0.68               0.084       511          60.7
#>   ...

Sort by "mean_score" to find bottlenecks (variables scoring poorly across the dataset), or by "pct_coverage" to find sparsely observed layers.

Data layer coverage

n_layers_scored in the result dataframe records how many variables contributed to the suitability score at each location. Low values indicate points where the score rests on few observations. GeoTIFF band 5 maps this spatially in QGIS.

Evidence weighting

By default the engine weights each variable by 1 / rank and normalises those weights over whichever variables are present in your data. That renormalisation makes absence indistinguishable from irrelevance: a survey measuring four of a species' seventeen scored variables gets a confident number on [0, 1] with nothing on it to say the rest was never observed.

# Reliability is measured, never assumed. Two independent estimates of the same
# quantity -- line against line, pass against pass.
rel <- c(
  depth              = reliability_from_replicates(depth_pass1, depth_pass2),
  slope              = reliability_from_replicates(slope_pass1, slope_pass2),
  substrate_hardness = reliability_from_replicates(hard_pass1,  hard_pass2)
)

ev <- predict_oyster(survey, "ostrea_edulis",
                     weighting   = "evidence",
                     reliability = rel)

attr(ev, "evidence_weights")     # the derivation, layer by layer
attr(ev, "evidence_coverage")    # what fraction of the evidence base you have
cat(weight_table_md(ev), sep = "\n")   # both, as markdown for a methods section

w = prior x reliability x spread. Reliability multiplies the biological prior and never replaces it — weighting on measured reliability alone promotes whatever the instrument happens to measure well, which is circular. A layer with no measured reliability gets no weight; a layer that scores the same everywhere gets no weight either, because it carries no information about where within this survey conditions differ. Coverage is reported separately and never multiplied into the per-cell score.

Locations below confidence_floor come back as "Insufficient data" rather than being handed a class they cannot support. suitability_rank and suitability_class_rank are kept alongside so the two weightings can be compared directly.

Fine-scale area analysis

area_summary() converts point-based scores into habitat area estimates at sub-hectare resolution — designed for restoration ecology where the difference between 200 m² and 800 m² matters.

s <- area_summary(result, cell_size_m = 10, viable_area_m2 = 100)

# Per-class breakdown in m²
s$class_summary

# Total suitable area
s$total["suitable_area_m2"]
s$total["pct_suitable"]

# Contiguous patches — flag those meeting OSPAR minimum viable area
s$patches[s$patches$viable, ]

Cell size is auto-estimated from median nearest-neighbour point spacing when cell_size_m = NULL. The viable_area_m2 threshold (default 100 m² = 0.01 ha) aligns with OSPAR oyster reef restoration guidance.

Tolerance visualisation

plot_tolerance() draws the mathematical scoring curve for any species/variable combination directly from the tolerance specification — no dataset needed. Useful for QA, stakeholder communication, and report figures.

# Scoring curve with colour-coded zones (optimal = green, excluded = red)
plot_tolerance("ostrea_edulis", "temperature")

# All four seasons overlaid
plot_tolerance("ostrea_edulis", "temperature", season = "all")

# Ocean acidification variables
plot_tolerance("ostrea_edulis", "ph")
plot_tolerance("magallana_gigas", "omega_aragonite")

# Compare DO tolerance: European flat vs brackish slipper oyster
plot_tolerance("ostrea_edulis",       "dissolved_oxygen")
plot_tolerance("crassostrea_iredalei","dissolved_oxygen")

For data-driven partial dependence curves from an actual survey, use sensitivity_analysis() instead.

Sensor readers

# Nortek Signature 500 ADCP
adcp <- read_nortek_adcp("adcp_export.csv")

# Nortek Aquadopp (moored or vessel-mounted)
adcp2 <- read_nortek_aquadopp("aquadopp_export.csv")

# Teledyne RDI binary PD0 (requires oce)
adcp3 <- read_rdi_adcp("transect.pd0")

# Ping 3DSS / BioBase bathymetric raster
bathy <- read_sonar_tif("bathymetry.tif")

# Bathymetric sounding XYZ point cloud
bathy2 <- read_soundings_xyz("soundings.xyz")

# Any standard CSV export (auto column matching)
ctd <- read_generic_csv("ctd_cast.csv")

# Merge multiple sensor sources onto a common spatial grid
survey <- merge_sensor_data(adcp = adcp, bathy = bathy2, ctd = ctd)

Validation and diagnostics

# Validate against known presence/absence records
val <- validate_against_records(result, records)
val$auc      # ROC-AUC
val$tss      # True Skill Statistic
val$f1       # F1 score
val$brier    # Brier score

# Spatial block cross-validation (avoids inflated AUC from autocorrelation)
cv <- spatial_block_cv(result, records, n_blocks = 5)

# Permutation variable importance (AUC drop per variable)
imp <- permutation_importance(result, records)

# Partial dependence curve for a single variable (data-driven)
sensitivity_analysis(result, records, variable = "temperature")

# Dataset-level weighted contribution summary (no records needed)
variable_impact(result, "ostrea_edulis")

# Multi-species comparison at the same grid
compare_species(survey_clean, species = c("ostrea_edulis", "magallana_gigas"))

Bayesian tolerance updating

Update species tolerance parameters from your own field observations:

fit <- update_species_tolerances(
  records     = field_data,
  species     = "ostrea_edulis",
  update_vars = c("temperature", "salinity", "depth")
)

# Subsequent predict_oyster() calls use the updated parameters automatically
result2 <- predict_oyster(df, "ostrea_edulis")

# Persist across sessions
save_tolerance_update("ostrea_edulis")

Additional modules

QGIS workflow

# Export five-band GeoTIFF + contour lines
predict_oyster(df, "ostrea_edulis", output_geotiff = "suitability.tif")

# Export matching colour ramp style file
export_qml_style("suitability.tif")

Load suitability.tif in QGIS, then drag the .qml file onto the layer to apply the standard oystermapR colour scheme instantly.

GeoTIFF bands:

Bands 3, 4, 6 and 7 are the ones that say how much to trust band 1. Band 1 will interpolate a confident-looking score across water nobody surveyed.

export_qml_style() writes two files: the band 1 colour ramp, and a companion <name>_insufficient.qml that styles band 7 as a grey overlay. Load the raster twice and drop one style on each to show unclassified ground as its own class.

Citation

If you use oystermapR in published work, please cite:

Tucker T. (2026). oystermapR: Predict and Map Oyster Growth Suitability from Environmental Data. R package version 1.6.0. https://github.com/trissyboats/oystermapR

License

GPL-3 © T Tucker

Free for research, education, and non-commercial use. Commercial entities wishing to embed oystermapR in a proprietary product should contact [email protected] to discuss a commercial licence.