Code
# Load required packages
library(tidyverse)
library(sf)
library(terra)
library(here)
library(tmap)
library(kableExtra)
library(testthat)November 15, 2024
Marine aquaculture represents a critical opportunity for sustainable food production as global seafood demand continues to rise. Unlike land-based meat production, marine aquaculture can provide protein with lower environmental costs when strategically located in suitable environments (Hall et al. 2011). Recent research has demonstrated that global seafood demand could theoretically be met using less than 0.015% of the global ocean area (Gentry et al. 2017), highlighting the importance of identifying optimal locations for development.
The US West Coast, with its diverse marine ecosystems and extensive Exclusive Economic Zones (EEZs), presents significant potential for marine aquaculture expansion. However, successful aquaculture operations depend largely on environmental conditions, particularly sea surface temperature and water depth, which determine species survival and growth rates.
This analysis evaluates the potential for marine aquaculture along the US West Coast by identifying EEZs with suitable environmental conditions for two commercially important species: oysters (Crassostrea spp.) and black abalone (Haliotis cracherodii). Using geospatial analysis of satellite-derived sea surface temperature and bathymetry data, we address three key questions:
By developing a generalizable workflow, this analysis can be applied to assess aquaculture potential for additional species, supporting evidence-based decision-making for sustainable marine resource development.
This analysis uses several R packages for geospatial data processing and visualization:
tidyverse — Data manipulation and visualizationsf — Spatial vector data handlingterra — Raster data processingtmap — Thematic mappingkableExtra — Table formattingWe use maritime boundary data for the US West Coast EEZs, which define the regions where the United States has rights to explore and exploit marine resources (Flanders Marine Institute 2024). The shapefile includes boundaries for California (Northern, Central, and Southern), Oregon, and Washington.
The EEZ data uses EPSG:4326 (WGS 84), a geographic coordinate reference system.
Average annual sea surface temperature data from 2008-2012 provides a multi-year baseline for identifying thermally suitable aquaculture locations (NOAA Coral Reef Watch 2024). The data comes from NOAA’s 5km Daily Global Satellite Sea Surface Temperature product, aggregated to annual averages and provided in degrees Kelvin.
# Load SST rasters for 2008-2012
# Define a function to read rast
read_tif <- function(year) {
rast(file.path("data", paste0("average_annual_sst_", year, ".tif")))
}
# Apply function
sst_list <- lapply(2008:2012, read_tif)
# Combine into raster stack
sst_stack <- rast(sst_list)
# Inspect properties
names(sst_stack)
st_crs(sst_stack["average_annual_sst_2008"])
summary(values(sst_stack))
global(sst_stack, fun = "isNA")We’ve loaded 5 years of annual SST data and combined them into a raster stack for further processing.
Ocean depth data comes from the General Bathymetric Chart of the Oceans (GEBCO) 2022 Grid (GEBCO Compilation Group 2022), providing global bathymetric measurements. Depth is a critical parameter for aquaculture site selection as different species have specific depth requirements based on their feeding behavior, oxygen requirements, and substrate preferences.
The bathymetry data also uses EPSG:4326, matching our SST and EEZ data.
Before proceeding with spatial analysis, we must ensure all datasets share a common coordinate reference system (CRS). Spatial operations like overlays and intersections will fail or produce incorrect results if the CRS doesn’t match.
# Check sst crs
crs(sst_stack, describe = TRUE)$code
# Set epsg explicitly from metadata
crs(sst_stack) <- "EPSG:4326"
# Check again
crs(sst_stack, describe = TRUE)$code
# Check extent makes sense for this CRS
ext(sst_stack)
# Verify sst and depth crs are the same
crs(sst_stack) == crs(depth)
# Check if wc_df crs matches
crs(sst_stack, describe = TRUE)$code == st_crs(wc_df)$epsgAll three datasets (SST, depth, and EEZ boundaries) now use EPSG:4326 (WGS 84), ensuring compatibility for subsequent spatial operations.
To create a stable baseline for identifying suitable temperatures, we calculate the mean SST across all five years (2008-2012). This multi-year average reduces the influence of short-term temperature anomalies and provides a more representative measure of typical thermal conditions.
The spatial pattern shows a clear north-south temperature gradient along the West Coast, with cooler waters in the north (Washington/Oregon) and warmer waters in Southern California.
The NOAA SST data is provided in degrees Kelvin. For easier interpretation and to match species requirement data (which uses Celsius), we convert to degrees Celsius by subtracting 273.15.
# Convert from Kelvin to Celsius (subtract 273.15)
sst_mean_c <- sst_mean - 273.15
# Extract values as vector and calculate summary statistics
sst_values <- as.vector(values(sst_mean_c))
sst_values <- sst_values[!is.na(sst_values)]
# Create summary statistics table
temp_summary <- data.frame(
Statistic = c("Minimum", "1st Quartile", "Median", "Mean", "3rd Quartile", "Maximum"),
Temperature_C = round(c(
min(sst_values),
quantile(sst_values, 0.25),
median(sst_values),
mean(sst_values),
quantile(sst_values, 0.75),
max(sst_values)
), 2)
)
kable(temp_summary,
caption = "Summary Statistics for Mean SST (2008-2012) in °C",
col.names = c("Statistic", "Temperature (°C)"),
align = c("l", "r")) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE,
position = "left")| Statistic | Temperature (°C) |
|---|---|
| Minimum | 4.98 |
| 1st Quartile | 12.38 |
| Median | 14.07 |
| Mean | 14.16 |
| 3rd Quartile | 16.06 |
| Maximum | 32.89 |
The temperature range along the West Coast spans approximately 8-21°C, which falls within the thermal tolerance of many marine aquaculture species.
To perform raster algebra (combining temperature and depth layers), we need the rasters to have identical spatial properties: extent, resolution, and CRS. The depth and SST rasters currently have different resolutions and extents, so we crop and resample the depth data to match the SST grid.
# Crop depth to SST extent using crop()
depth_crop <- crop(depth, sst_mean_c)
# Resample depth to match SST resolution using resample() with method = "near"
depth_resampled <- resample(depth_crop, sst_mean_c, method = "near")
# Verify alignment
res(depth_resampled) == res(sst_mean_c)
ext(depth_resampled) == ext(sst_mean_c)
crs(depth_resampled) == crs(sst_mean_c)We use the “nearest neighbor” resampling method, which is appropriate for preserving discrete depth values without interpolation artifacts. The depth and SST rasters are now perfectly aligned.
Oysters are among the most widely cultivated marine organisms globally, valued for their hardiness and ability to filter feed in coastal waters. Based on data from SeaLifeBase (Palomares, M.L.D. and D. Pauly 2024), oysters require:
These requirements reflect oyster’s tolerance for relatively warm waters and their preference for shallow to moderate depths where food availability is high and substrate is accessible.
We use binary reclassification to identify areas meeting both temperature and depth criteria. For each environmental layer, cells within the suitable range are assigned a value of 1, while unsuitable cells receive 0. Multiplying the two layers together produces a final suitability map where only locations meeting both criteria receive a value of 1.
# Reclassify SST for oysters using classify()
# Create rcl matrix: <11°C = 0, 11-30°C = 1, >30°C = 0
sst_rcl_mat_oyster <- matrix(c(-Inf, 11, 0,
11, 30, 1,
30, Inf, 0),
ncol = 3,
byrow = TRUE)
sst_rcl_oyster <- classify(sst_mean_c, sst_rcl_mat_oyster)
# Reclassify depth for oysters using classify()
# Create rcl matrix: <-70m = 0, -70m to 0m = 1, >0m = 0
depth_rcl_mat_oyster <- matrix(c(-Inf, -70, 0,
-70, 0, 1,
0, Inf, 0),
ncol = 3,
byrow = TRUE)
depth_rcl_oyster <- classify(depth_resampled, depth_rcl_mat_oyster)
# Combine suitability layers
suitable_zones_oyster <- sst_rcl_oyster * depth_rcl_oyster
# Visualize with plot()
plot(suitable_zones_oyster,
main = "Oyster Suitability\n(1 = Suitable, 0 = Unsuitable)",
cex.main = 0.9,
col = c("gray90", "darkgreen"))The map reveals extensive suitable habitat for oysters along the West Coast, particularly in California’s coastal waters where temperatures and depths align with oyster requirements.
Black abalone is a marine gastropod native to the Pacific coast, historically important both ecologically and commercially. However, the species has experienced significant population declines and is now listed as endangered. Understanding suitable habitat is crucial for both aquaculture development and conservation planning. Based on SeaLifeBase (Palomares, M.L.D. and D. Pauly 2024), black abalone requires:
These narrower requirements reflect black abalone’s preference for cool, shallow rocky intertidal and subtidal zones typical of the California coastline.
We apply the same binary reclassification approach used for oysters, but with black abalone’s more restrictive temperature and depth ranges.
# Reclassify SST for black abalone using classify()
# Suitable: 12.2-18.6°C = 1, Unsuitable: <12.2°C or >18.6°C = 0
sst_rcl_mat_abalone <- matrix(c(-Inf, 12.2, 0,
12.2, 18.6, 1,
18.6, Inf, 0),
ncol = 3,
byrow = TRUE)
# Apply reclassification matrix to mean SST raster
sst_rcl_abalone <- classify(sst_mean_c, sst_rcl_mat_abalone)
# Reclassify depth for black abalone using classify()
# Suitable: 0-6m below sea level = -6 to 0 in raster values
depth_rcl_mat_abalone <- matrix(c(-Inf, -6, 0,
-6, 0, 1,
0, Inf, 0),
ncol = 3,
byrow = TRUE)
# Apply reclassification matrix to resampled depth raster
depth_rcl_abalone <- classify(depth_resampled, depth_rcl_mat_abalone)
# Combine suitability layers
suitable_zones_abalone <- sst_rcl_abalone * depth_rcl_abalone
# Visualize with plot()
plot(suitable_zones_abalone,
main = "Black Abalone Suitability\n(1 = Suitable, 0 = Unsuitable)",
cex.main = 0.9,
col = c("gray90", "darkblue"))The map shows much more limited suitable habitat for black abalone compared to oysters, with suitable areas concentrated in central California where temperatures remain within the narrow thermal window and very shallow depths occur.
To calculate suitable area within each EEZ region, we need to convert the EEZ vector boundaries to raster format. This allows us to use zonal statistics to sum suitable areas within each region.
Each EEZ region is now represented in raster format with cells assigned their corresponding region ID, enabling region-specific area calculations.
We use zonal statistics to sum the suitable area within each EEZ region. The process involves:
# Mask suitable areas to EEZ using mask()
eez_suitable_oyster <- mask(suitable_zones_oyster, eez_rast)
# Calculate cell area in km² using cellSize()
cell_area_km2 <- cellSize(eez_suitable_oyster, unit = "km")
# Multiply suitable locations by cell area to get suitable area per cell
suitable_area_oyster <- eez_suitable_oyster * cell_area_km2
# Sum area by EEZ region using zonal()
area_by_region_oyster <- zonal(suitable_area_oyster, eez_rast, fun = "sum", na.rm = TRUE)
# Rename area column for clarity
names(area_by_region_oyster)[2] <- "suitable_area_km2"
# Join results with EEZ sf object to get region names
area_by_region_oyster_joined <- eez %>%
left_join(area_by_region_oyster, by = "rgn_id") %>%
arrange(desc(suitable_area_km2))
oyster_results <- area_by_region_oyster_joined %>%
st_drop_geometry() %>%
select(rgn, suitable_area_km2) %>%
head(10)
# Display top 10 regions in table with kable()
kable(oyster_results,
caption = "Top 10 EEZ Regions for Oyster Aquaculture",
col.names = c("Region", "Suitable Area (km²)"),
digits = 2) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Region | Suitable Area (km²) |
|---|---|
| Central California | 4923.15 |
| Southern California | 4096.54 |
| Washington | 3224.74 |
| Oregon | 1533.09 |
| Northern California | 438.15 |
Key Finding: Central California leads with over 4,000 km² of suitable habitat, followed by Southern California and Washington, together accounting for more than 90% of all suitable oyster habitat on the West Coast.
# Test that area calculations are reasonable
test_that("Oyster area calculations are valid", {
# Check total area > 0
expect_true(sum(area_by_region_oyster$suitable_area_km2, na.rm = TRUE) > 0)
# Check all values >= 0
expect_true(all(area_by_region_oyster$suitable_area_km2 >= 0, na.rm = TRUE))
# Check that we have results for multiple regions
expect_true(nrow(area_by_region_oyster) > 0)
})The map below visualizes how suitable area is distributed across West Coast EEZs, with darker blues indicating greater total suitable area.
# Create map using tmap
tm_shape(area_by_region_oyster_joined) +
tm_polygons(
fill = "suitable_area_km2",
fill.scale = tm_scale_continuous(values = "brewer.yl_gn_bu"),
fill.legend = tm_legend(title = "Suitable Area (km²)")
) +
tm_borders(col = "gray40", lwd = 0.5) +
tm_title("Oyster Aquaculture Suitability by EEZ Region") +
tm_compass(type = "4star", position = c("left", "bottom")) +
tm_scalebar(position = c("left", "bottom")) +
tm_layout(legend.outside = TRUE,
legend.outside.position = "right",
frame = FALSE)
The map confirms that California EEZs dominate oyster aquaculture potential, with Central and Southern California showing the highest concentration of suitable habitat.
We repeat the same zonal statistics workflow for black abalone, calculating the total suitable area within each EEZ region.
# Mask suitable areas to EEZ using mask()
eez_suitable_abalone <- mask(suitable_zones_abalone, eez_rast)
# Calculate cell area in km² using cellSize()
cell_area_km2_abalone <- cellSize(eez_suitable_abalone, unit = "km")
# Multiply suitable locations by cell area to get suitable area per cell
suitable_area_abalone <- eez_suitable_abalone * cell_area_km2_abalone
# Sum area by EEZ region using zonal()
area_by_region_abalone <- zonal(suitable_area_abalone, eez_rast, fun = "sum", na.rm = TRUE)
# Rename area column for clarity
names(area_by_region_abalone)[2] <- "suitable_area_km2"
# Join results with EEZ sf object to get region names
area_by_region_abalone_joined <- eez %>%
left_join(area_by_region_abalone, by = "rgn_id") %>%
arrange(desc(suitable_area_km2))
abalone_results <- area_by_region_abalone_joined %>%
st_drop_geometry() %>%
select(rgn, suitable_area_km2) %>%
head(10)
# Display top 10 regions in table with kable()
kable(abalone_results,
caption = "Top 10 EEZ Regions for Black Abalone Aquaculture",
col.names = c("Region", "Suitable Area (km²)"),
digits = 2) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Region | Suitable Area (km²) |
|---|---|
| Central California | 186.24 |
| Washington | 103.47 |
| Southern California | 89.01 |
| Northern California | 16.22 |
| Oregon | 14.90 |
Key Finding: Black abalone has dramatically less suitable habitat than oysters—only 199.88 km² total across the entire West Coast. Central California accounts for 76% (152.37 km²) of all suitable habitat, while Oregon and Northern California have essentially no suitable areas.
# Test that area calculations are reasonable
test_that("Black Abalone area calculations are valid", {
# Check total area > 0
expect_true(sum(area_by_region_abalone$suitable_area_km2, na.rm = TRUE) > 0)
# Check all values >= 0
expect_true(all(area_by_region_abalone$suitable_area_km2 >= 0, na.rm = TRUE))
# Check that we have results for multiple regions
expect_true(nrow(area_by_region_abalone) > 0)
})# Create map using tmap
tm_shape(area_by_region_abalone_joined) +
tm_polygons(
fill = "suitable_area_km2",
fill.scale = tm_scale_continuous(values = "brewer.yl_gn_bu"),
fill.legend = tm_legend(title = "Suitable Area (km²)")
) +
tm_borders(col = "gray40", lwd = 0.5) +
tm_title("Black Abalone Aquaculture Suitability by EEZ Region") +
tm_compass(type = "4star",
position = c("left", "bottom")) +
tm_scalebar(position = c("left", "bottom")) +
tm_layout(legend.outside = TRUE,
legend.outside.position = "right",
frame = FALSE)
The map reveals stark geographic concentration, with black abalone suitable habitat almost entirely restricted to Central California EEZ, reflecting the species’ narrow thermal tolerance and very shallow depth requirements.
The analysis workflow developed for oysters and black abalone can be generalized into a single function that accepts species-specific parameters. This enables rapid assessment of aquaculture suitability for any marine species with known temperature and depth requirements.
The function takes five parameters:
min_temp, max_temp — Temperature range in degrees Celsiusmin_depth, max_depth — Depth range in meters (negative values below sea level)species_name — Species name for labeling outputscalculate_suitability <- function(min_temp, max_temp, min_depth, max_depth, species_name) {
# Validate inputs
if (min_temp >= max_temp) {
stop("min_temp must be less than max_temp")
}
if (min_depth >= max_depth) {
stop("min_depth must be less than max_depth")
}
if (!is.character(species_name) || species_name == "") {
stop("species_name must be a non-empty string")
}
# Create reclassification matrix for sea surface temperature
sst_rcl_matrix <- matrix(c(-Inf, min_temp, 0,
min_temp, max_temp, 1,
max_temp, Inf, 0),
ncol = 3,
byrow = TRUE)
# Apply reclassification matrix to mean SST raster
sst_rcl <- classify(sst_mean_c, sst_rcl_matrix)
# Create reclassification matrix for depth
depth_rcl_matrix <- matrix(c(-Inf, min_depth, 0,
min_depth, max_depth, 1,
max_depth, Inf, 0),
ncol = 3,
byrow = TRUE)
# Apply reclassification matrix to resampled depth raster
depth_rcl <- classify(depth_resampled, depth_rcl_matrix)
# Combine temperature and depth suitability using multiplication
suitable_zones <- sst_rcl * depth_rcl
# Mask suitable zones to only include areas within EEZ boundaries
eez_suitable <- mask(suitable_zones, eez_rast)
# Calculate the area of each raster cell in square kilometers
cell_area <- cellSize(eez_suitable, unit = "km")
# Multiply suitable locations (1s) by cell area (unsuitable cells with have 0 area)
suitable_area <- eez_suitable * cell_area
# Sum suitable area within each EEZ region
area_by_region <- zonal(suitable_area, eez_rast, fun = "sum", na.rm = TRUE)
# Rename the area column to be more descriptive
names(area_by_region)[2] <- "suitable_area_km2"
# Join area calculations with EEZ spatial data to add region names
results_joined <- eez %>%
left_join(area_by_region, by = "rgn_id") %>%
arrange(desc(suitable_area_km2))
# Create a summary table showing top 10 regions
results_table <- results_joined %>%
st_drop_geometry() %>%
select(rgn, suitable_area_km2) %>%
head(10)
# Create choropleth map showing suitable area by EEZ region
map <- tm_shape(results_joined) +
tm_polygons(
fill = "suitable_area_km2",
fill.scale = tm_scale_continuous(values = "brewer.yl_gn_bu"),
fill.legend = tm_legend(title = "Suitable Area (km²)")
) +
tm_borders(col = "gray40", lwd = 0.5) +
tm_title(paste(species_name, "Aquaculture Suitability by EEZ Region")) +
tm_compass(type = "4star", position = c("right", "top")) +
tm_scalebar(stack = "vertical", position = c("left", "bottom")) +
tm_layout(legend.outside = TRUE,
legend.outside.position = "right",
frame = FALSE)
return(list(
results = results_table,
full_data = results_joined,
map = map
))
}We can verify the function produces identical results to our step-by-step analysis by applying it to oysters with the same parameters.
# Call function with oyster parameters
oyster_output <- calculate_suitability(
min_temp = 11,
max_temp = 30,
min_depth = -70,
max_depth = 0,
species_name = "Oyster"
)
# Display results table
kable(oyster_output$results,
caption = "Top 10 EEZ Regions for Oyster Aquaculture (Function Output)",
col.names = c("Region", "Suitable Area (km²)"),
digits = 2) %>%
kable_styling(bootstrap_options = c("striped", "hover"))| Region | Suitable Area (km²) |
|---|---|
| Central California | 4923.15 |
| Southern California | 4096.54 |
| Washington | 3224.74 |
| Oregon | 1533.09 |
| Northern California | 438.15 |

The function successfully replicates our earlier oyster analysis, confirming Central California’s dominance.

| Region | Suitable Area (km²) |
|---|---|
| Central California | 186.24 |
| Washington | 103.47 |
| Southern California | 89.01 |
| Northern California | 16.22 |
| Oregon | 14.90 |
The function also correctly identifies black abalone’s highly restricted suitable habitat, concentrated almost entirely in Central California.
Oysters demonstrate substantial aquaculture potential across the US West Coast, with over 11,000 km² of suitable habitat identified across all five EEZ regions:
These three top regions (Central CA, Southern CA, and Washington) together account for over 90% of all suitable oyster habitat on the West Coast. The broad temperature tolerance (11-30°C) and moderate depth range (0-70m) enable oysters to thrive across diverse coastal environments.
Black abalone presents a stark contrast, with only 199.88 km² of suitable habitat total—approximately 50 times less than oysters:
This dramatic restriction reflects black abalone’s narrow thermal tolerance (12.2-18.6°C) and very shallow depth requirement (0-6m below sea level). The species is geographically constrained to cool, shallow coastal waters predominantly found in Central California.
The 50-fold difference in suitable habitat between oysters and black abalone highlights how species-specific environmental requirements shape aquaculture potential:
This analysis demonstrates that Central California EEZ represents the most versatile region for marine aquaculture development on the West Coast, supporting substantial habitat for both broad-tolerance species (oysters) and narrow-tolerance species (black abalone).
For aquaculture development priorities:
The generalizable function developed here enables rapid assessment of additional candidate species, facilitating evidence-based site selection and multi-species aquaculture planning.
This analysis was completed for EDS 223: Geospatial Analysis and Remote Sensing at the Bren School of Environmental Science & Management, UC Santa Barbara (November 2024). The workflow demonstrates the application of raster processing, spatial statistics, and reproducible geospatial analysis to real-world marine resource management questions.
@online{miller2024,
author = {Miller, Emily},
title = {Prioritizing {Marine} {Aquaculture} {Locations} on the {US}
{West} {Coast}},
date = {2024-11-15},
url = {https://rellimylime.github.io/posts/eds223-final2/},
langid = {en}
}