This blog explores the environmental suitability of oyster aquaculture along the U.S. West Coast. Using sea surface temperature (SST) data, bathymetry rasters, and Exclusive Economic Zone (EEZ) boundaries, we identify areas most favorable for oyster cultivation. By combining geospatial modeling with cartographic visualization, we highlight how environmental data science can inform sustainable aquaculture planning.
Highlights of the Analysis
Raster Reclassification: Converted SST and depth rasters into binary suitability maps based on oyster thresholds. Map Algebra: Combined multiple criteria (temperature × depth) to produce a composite suitability raster. Zonal Statistics: Summed suitable areas within EEZ boundaries to quantify aquaculture potential by region. Cartographic Visualization: Produced clear maps with legends, units, and orientation cues to communicate results effectively.
Dataset Descriptions
Sea Surface Temperature (2008–2012) NOAA Coral Reef Watch 5‑km SST data, averaged across five years and converted from Kelvin to Celsius. Accessed Nov. 2025.
Bathymetry GEBCO gridded bathymetry dataset, reprojected and resampled to align with SST rasters. Accessed Nov. 2025.
Exclusive Economic Zones (EEZs) Marine Regions shapefiles defining U.S. West Coast EEZ boundaries. Accessed Nov. 2025.
Species Requirements Oyster depth and temperature preferences from SeaLifeBase. Accessed Nov. 2025.
Findings
Oyster aquaculture suitability is concentrated in specific EEZs along the West Coast.
Combining SST and depth thresholds provides a robust method for identifying viable cultivation zones.
Zonal statistics allow policymakers to compare aquaculture potential across regions.
Conclusion
This analysis demonstrates how spatial data science can inform sustainable aquaculture planning. By integrating environmental thresholds with geospatial datasets, we can identify where oysters are most likely to thrive. Future work could expand this framework to other species — such as grapevines requiring specific elevations or trees adapted to certain climates — showcasing the versatility of raster reclassification and zonal statistics in environmental modeling.
References
SeaLifeBase. (2025). Species depth and temperature requirements. https://www.sealifebase.ca/search.php [Accessed Nov. 2025]
# Import libraries for geospatial and data analysissuppressPackageStartupMessages({library(tidyverse)library(stars)library(terra)library(sf)library(tmap)library(knitr)library(kableExtra)})
Warning: package 'tidyverse' was built under R version 4.5.2
Warning: package 'ggplot2' was built under R version 4.5.2
Warning: package 'stars' was built under R version 4.5.2
Warning: package 'sf' was built under R version 4.5.2
Warning: package 'terra' was built under R version 4.5.2
Warning: package 'knitr' was built under R version 4.5.2
Code
# Function to map EEZ suitability for a speciesmap_species_suitability <-function(sst_min, sst_max, depth_min, depth_max, species_name) {# Load in data regions_clean <-st_read("data/wc_regions_clean.shp", quiet =TRUE) depth <-rast("data/depth.tif") sst_stack <-c(rast("data/average_annual_sst_2008.tif"),rast("data/average_annual_sst_2009.tif"),rast("data/average_annual_sst_2010.tif"),rast("data/average_annual_sst_2011.tif"),rast("data/average_annual_sst_2012.tif") )# Process the SST and convert to celsius sst_mean <-mean(sst_stack) sst_celsius <- sst_mean -273.15# Align depth to SST depth_proj <-project(depth, crs(sst_celsius)) depth_crop <-crop(depth_proj, sst_celsius) depth_resample <-resample(depth_crop, sst_celsius, method ="near")# Create suitability classification matrixes sst_suit <-classify(sst_celsius,rbind(c(-Inf, sst_min, 0),c(sst_min, sst_max, 1),c(sst_max, Inf, 0))) depth_suit <-classify(depth_resample,rbind(c(-Inf, depth_min, 0),c(depth_min, depth_max, 1),c(depth_max, Inf, 0))) suitability <- sst_suit * depth_suit# Create EEZ zonal stats# Convert cleaned EEZ regions to vector eez_v <-vect(regions_clean)# Reproject the EEZ vector to match CRS of sst_celsius eez_v <-project(eez_v, crs(sst_celsius))# Add unique ID for each EEZ zone for later join eez_v$eez_id <-seq_len(nrow(eez_v))# Rasterize the EEZ polygon so each cell in raster is assigned an ID that is within eez_r <-rasterize(eez_v, sst_celsius, field ="eez_id")# Calculate the area of each raster in square kilometers cell_area <-cellSize(sst_celsius, unit ="km")# Multiply the suitability raster by cell area to get area of suitable habitat per cell area_suitable <- suitability * cell_area# Perform zonal math by getting the sum of the suitable area in each EEZ zone area_by_eez <-zonal(area_suitable, eez_r, fun ="sum", na.rm =TRUE)# Convert above results into a dataframe and rename column area_by_eez <-as.data.frame(area_by_eez) |>rename(suitable_km2 = mean)# Convert EEZ vector back to an sf obejct and join stats by EEZ ID eez_sf_out <-st_as_sf(eez_v) |>left_join(area_by_eez, by ="eez_id")# Plot maptm_shape(eez_sf_out) +tm_polygons(fill ="suitable_km2",fill.scale =tm_scale(values ="viridis"),fill.legend =tm_legend(title ="Suitable area (km²)")) +tm_title(paste(species_name, "Suitability Across West Coast EEZs"),position =c("center", "top")) +tm_scalebar(position =c("left", "bottom")) +tm_compass(position =c("right", "top"))}
Code
# Use the function for finding suitable area for Geoducks# map_species_suitability(sst_min, sst_max, depth_min, depth_max, species_name)map_species_suitability(sst_min =0, sst_max =80, depth_min =-100, depth_max =0, species_name ="Geoduck")
[plot mode] fit legend/component: Some legend items or map compoments do not
fit well, and are therefore rescaled.
ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
Analysis 1: Mean Sea Surface Temperature
Import and Prepare the Data
Code
# Read in Exclusive Economic Zone (EEZ) shapefileregions_clean <-st_read("data/wc_regions_clean.shp")
Reading layer `wc_regions_clean' from data source
`C:\Users\joshu\Documents\MEDS\Career\Awoo56709.github.io\posts\Aquaculture\data\wc_regions_clean.shp'
using driver `ESRI Shapefile'
Simple feature collection with 5 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -129.1635 ymin: 30.542 xmax: -117.097 ymax: 49.00031
Geodetic CRS: WGS 84
Code
# Read in depth rasterdepth <-rast("data/depth.tif")# Read in Sea Surface Temperature rasterssst_2008 <-rast("data/average_annual_sst_2008.tif")sst_2009 <-rast("data/average_annual_sst_2009.tif")sst_2010 <-rast("data/average_annual_sst_2010.tif")sst_2011 <-rast("data/average_annual_sst_2011.tif")sst_2012 <-rast("data/average_annual_sst_2012.tif")
Process the Data
We assume WGS84 (EPSG:4326) as the common CRS. If CRS mismatches occur, analysis would be invalid.
# Compute the mean of SST across yearssst_mean <-mean(sst_stack)
Code
# Convert Kelvin to Celsiussst_celsius <- sst_mean -273.15# Replace NA values with 0sst_celsius[is.na(sst_celsius)] <-0# Add warning if the range is wrongif (max(values(sst_celsius), na.rm=TRUE) >40) warning("SST unusually high — check conversion.")
Code
# Convert SST raster to stars object for a more robust plotsst_stars <-st_as_stars(sst_celsius)tm_shape(sst_stars) +tm_raster(fill.scale =tm_scale(values ="viridis"),fill.legend =tm_legend(title ="Temperature (C)")) +tm_title("Mean Sea Surface Temperature (°C, 2008–2012)", position =c("center", "top")) +tm_scalebar(position =c("left", "bottom")) +tm_compass(position =c("right", "bottom"))
[tm_raster()] Arguments `fill.scale` and `fill.legend` unknown.
[plot mode] fit legend/component: Some legend items or map compoments do not
fit well, and are therefore rescaled.
ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.
Reflection on Mean Sea Surface Temperature Raster
The mean sea surface temperature (SST) raster from the years 2008 - 2012 provides insight into the temperature variation along the west coast of the united states from California to Washington state. Because temperature data is continuous, a color scale was used to illustrate temperature variability. The above illustrate is not varied because bins of mean temperature across five years of sea surface temperature.
Analysis 2: Combined Suitability Raster
Code
# Align depth raster to the SST raster# Reproject CRSdepth_proj <-project(depth, crs(sst_celsius)) # Crop to SST extentdepth_crop <-crop(depth_proj, sst_celsius)# Resample resolutiondepth_resample <-resample(depth_crop, sst_celsius, method ="near")
Resampling with the nearest neighbor preserves categorical depth bins.
Code
# Verify alignments of CRS, resolution and extentcrs(sst_celsius) ==crs(depth_resample)
[1] TRUE
Code
if (!identical(res(sst_celsius), res(depth_resample))) stop("Resolution mismatch.")
Code
# Create matrix to reclassify SST suitability for oysters based on 11–30 Csst_suit <-classify(sst_celsius,rbind(c(-Inf, 11, 0), # Unsuitable temp below 11 Cc(11, 30, 1), # Suitable temp at 11 - 30 Cc(30, Inf, 0))) # Unsuitable temp above 30 C# Create matrix to reclassify depth suitability for oysters based on 0–70 mdepth_suit <-classify(depth_resample,rbind(c(-Inf, -70, 0), # Unsuitable elevation above sea levelc(-70, 0, 1), # Suitable between -70 - 0 m sea levelc(0, Inf, 0))) # Unsuitable elevation above sea level# Test to verify our matrix of binary outputtestthat::expect_true(all(unique(values(sst_suit)) %in%c(0,1)))
We assume depths are negative, so suitable oyster depths are −70 to 0 m. If positive, thresholds would be 0–70.
Code
# Combine suitability (map algebra: both must be 1)suitability <- sst_suit * depth_suit# Convert raster to dataframesuit_df <-as.data.frame(suitability, xy =TRUE, na.rm =TRUE)# Plot with ggplotggplot(suit_df) +geom_raster(aes(x = x, y = y, fill =factor(mean))) +scale_fill_manual(values =c("0"="purple", "1"="gold"),name ="Oyster suitability",labels =c("Unsuitable", "Suitable")) +labs(title ="Combined Suitability (Oysters)") +theme_minimal()
Reflection on Combined Suitability Raster
The bathymetry raster illustrates ocean depth across the study area. With the given information, suitability classification was determined by depth preferences of oysters and the categories mapped “suitable” and “unsuitable” was based on binary outcomes of the suitable matrix sea depth. This binary map showcases environments suitable for oyster aquaculture.
Analysis 3: Suitability by EEZ
Code
# Convert EEZ polygons to terra vectoreez_v <-vect(regions_clean)# Align CRSeez_v <-project(eez_v, crs(sst_celsius))# Rasterize EEZ IDs to each polygoneez_v$eez_id <-seq_len(nrow(eez_v))eez_r <-rasterize(eez_v, sst_celsius, field ="eez_id")# Calculate cell area in km^2cell_area <-cellSize(sst_celsius, unit ="km")# Multiply suitable area raster by cell to get suitable area per cellarea_suitable <- suitability * cell_area# Zonal sum of suitable area per EEZarea_by_eez <-zonal(area_suitable, eez_r, fun ="sum", na.rm =TRUE)# Convert to data frame and rename columns correctlyarea_by_eez <-as.data.frame(area_by_eez) |> dplyr::rename(suitable_km2 = mean)
Summing suitable cell areas by EEZ allows us to rank zones by aquaculture potential
Code
# Join zonal statistics to EEZ polygonseez_sf_out <-st_as_sf(eez_v) |>left_join(area_by_eez, by ="eez_id")
Code
# Plot map of suitable area by EEZeez_sf <-st_as_sf(eez_v) |>left_join(area_by_eez, by ="eez_id")ggplot(eez_sf) +geom_sf(aes(fill = suitable_km2), color ="black", size =0.2) +scale_fill_viridis_c(option ="C", name ="Suitable area (km^2)") +labs(title ="Oyster Suitability Across West Coast EEZs") +theme_minimal() +theme(plot.title.position ="plot", # move title above plot panelplot.title =element_text(margin =margin(b =10)) # add spacing below title )
Reflection on Oyster Suitability Across West Coast EEZs
This map aggregates suitability into Exclusive Economic Zones (EEZs), showing total suitability across the zones and illustrates aquaculture potential. We see that based on the table below that EEZ ID magnitude of suitability from largest to smallest are EEZ IDs: 3,4,5,1 and 2 illustrating the southern west coast tend to be an optimal environment for Oyster aquaculture.
Code
# Create table of suitable areas by EEZarea_by_eez %>%arrange(desc(suitable_km2)) %>% knitr::kable(digits =2, caption ="Suitable Area by EEZ (km^2)") |> kableExtra::kable_classic(full_width =FALSE, position ="center")
