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.
| Key | Species | Common name | Region |
|---|---|---|---|
ostrea_edulis |
Ostrea edulis | European Flat Oyster | NE Atlantic / Mediterranean |
magallana_gigas |
Magallana gigas | Pacific Oyster | Global aquaculture |
crassostrea_angulata |
Crassostrea angulata | Portuguese Oyster | Iberian Peninsula |
ostrea_stentina |
Ostrea stentina | Denticulate Flat Oyster | Mediterranean / W Africa |
ostrea_lurida |
Ostrea lurida | Olympia Oyster | Pacific NW (N America) |
crassostrea_virginica |
Crassostrea virginica | Eastern Oyster | NW Atlantic |
saccostrea_glomerata |
Saccostrea glomerata | Sydney Rock Oyster | E Australia |
magallana_sikamea |
Magallana sikamea | Kumamoto Oyster | Japan / Pacific NW |
magallana_ariakensis |
Magallana ariakensis | Suminoe Oyster | E Asia |
crassostrea_hongkongensis |
Crassostrea hongkongensis | Hong Kong Oyster | S China Sea |
crassostrea_nippona |
Crassostrea nippona | Iwagaki Oyster | Japan / Korea |
crassostrea_belcheri |
Crassostrea belcheri | Tropical Rock Oyster | SE Asia |
ostrea_chilensis |
Ostrea chilensis | Chilean Oyster | S Chile / New Zealand |
ostrea_denselamellosa |
Ostrea denselamellosa | Korean Flat Oyster | E Asia |
ostrea_angasi |
Ostrea angasi | Australian Flat Oyster | S Australia |
saccostrea_cucullata |
Saccostrea cucullata | Hooded Rock Oyster | Indian Ocean / W Pacific |
crassostrea_iredalei |
Crassostrea iredalei | Slipper Oyster | SE Asia (brackish) |
# 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.
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.).
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.
| Variable | Recognised column names |
|---|---|
| Temperature (°C) | temperature, temp,
temp_c |
| Salinity (PSU) | salinity, sal, salinity_psu,
psu, sal_psu |
| Dissolved oxygen (mg/L) | dissolved_oxygen, do, do_mgl,
oxygen, dissolved_o2,
do_mg_l |
| pH | ph, ph_nbs, ph_total,
seawater_ph |
| Total alkalinity (µmol/kg) | alkalinity, total_alkalinity,
ta, talk, alk_umol_kg |
| Aragonite saturation (Ω) | omega_aragonite, omega_arag,
aragonite_saturation, arag_sat,
omega |
| Depth (m) | depth, depth_m |
| Current velocity (m/s) | current_velocity, velocity,
current |
| Shear stress (N/m²) | shear_stress, tau,
bed_shear |
| Chlorophyll-a (µg/L) | chlorophyll_a, chla,
chlorophyll |
| Turbidity (NTU) | turbidity, ntu, turb |
| Slope (degrees) | slope, slope_deg |
| Substrate hardness | substrate_hardness, hardness |
A sample dataset is included at
inst/extdata/sample_survey.csv.
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.
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.
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.
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 sectionw = 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.
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.
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.
# 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)# 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"))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")| Function | Purpose |
|---|---|
variable_impact() |
Ranked variable contribution table for QA and survey planning |
area_summary() |
Fine-scale area estimates in m² with contiguous patch analysis |
plot_tolerance() |
Tolerance scoring curve from species parameters (no data required) |
calculate_aragonite() |
In-house aragonite saturation (Ω_arag) from pH and alkalinity |
score_wave_exposure() |
JONSWAP wave height from fetch and wind speed |
score_sediment_stability() |
Shields parameter sediment mobility analysis |
score_larval_connectivity() |
Hybrid Gaussian kernel + OpenDrift connectivity matrix |
score_predation_risk() |
Starfish, crab, and snail predation pressure |
score_hab_risk() |
Harmful algal bloom risk (PSP/ASP/DSP/AZP) |
score_anthropogenic_disturbance() |
Bottom trawling, anchor damage, dredging |
add_shellfish_classification() |
UK/EU harvesting area classification (A/B/C) |
compare_species() |
Side-by-side suitability across multiple species |
composite_seasonal() |
Merge summer/winter surveys into a composite score |
generate_report() |
Export a formatted PDF or HTML report |
evidence_weights() |
Weights from published rank x measured reliability x realised spread |
evidence_coverage() |
Share of the species’ evidence base the survey actually observed |
evidence_suitability() |
Re-score on evidence weights; curves untouched |
reliability_from_replicates() |
Reliability from real repeat coverage |
reliability_from_perturbation() |
Reliability from a measured error model |
confidence_score() |
Per-location confidence from weighted components |
weight_table_md() |
Weight derivation and coverage ledger as markdown |
# 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:
| Band | Name | Description |
|---|---|---|
| 1 | suitability |
Continuous suitability score [0, 1] |
| 2 | excluded_mask |
1 = hard-excluded (failed exclusion threshold) |
| 3 | n_observations |
Number of survey points per raster cell |
| 4 | dist_to_nearest_m |
Distance to nearest survey point (m) |
| 5 | n_layers_scored |
Number of environmental variables that contributed to the score |
| 6 | confidence |
Per-location confidence [0, 1] under evidence weighting;
NA under rank weighting |
| 7 | insufficient_mask |
1 = the model declined to classify this location |
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.
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
GPL-3 © T Tucker
Free for research, education, and non-commercial use. Commercial entities wishing to embed oystermapR in a proprietary product should contact tristantucker48@gmail.com to discuss a commercial licence.