Prioritizing Potential Finfish Species
With the U.N. predicting that the global human population will reach approximately 9.7 billion people by the year 2050, access to basic resources (particularly food and water) will continue to become increasing pressing issues (United Nations 2024). The global assessment of the state of food security and nutrition in 2022 highlights the severity of this issue at current population levels (~8.2 billion):
“It is estimated that between 691 and 783 million people in the world faced hunger in 2022. Considering the midrange (about 735 million), 122 million more people faced hunger in 2022 than in 2019, before the pandemic (Food and Agriculture Organization of the United Nations 2023).”
This highlights the increasingly urgent need to address the lack of global food security. As a result, marine aquaculture has become the center of research; with many notable characteristics that make it a sustainable alternative to meet the growing demands for protein sources (Hall et al. 2011; Nash et al. 2021). Gentry et al. mapped potential locations that would be suitable for marine aquaculture around the world, and they found that global seafood demand could be met using less than 0.015% of the global ocean area (Gentry et al. 2017). Additionally, an analysis of systems that have adopted an integrated multi-trophic aquaculture (IMTA) model show promise as a sustainable approach to open-ocean marine aquaculture (Troell et al. 2010).
Gentry et al. identified potential locations based on human and environmental constraints. This included excluding areas; that would negatively impact the growth of bivalves (low phytoplanktonic food availability) and finfish (low dissolve oxygen), were too expensive (too deep >200m) to anchor farms, were already allocated to other uses, or were in high-density shipping areas (Gentry et al. 2017).
Based on this analysis, we aim to identify locations along the west coast of the U.S. that would be optimal for the establishment of a marine aquaculture farm for both oysters and various finfish species based on two environmental parameters that impact optimal growth:
Research has shown that the following conditions are required for optimal oyster growth:
Potential finfish species were selected based on the current commercial fishes managed by the NOAA West Coast operation (NOAA Fisheries: West Coast 2024). Highly migratory species (such as sharks) were removed from consideration resulting in a list of 26 potential finfish species. We then used FishBase to compile a table of minimum temperature, maximum temperature, minimum depth, and maximum depth for each species (FishBase Consortium 2024). If a temperature range was not explicitly listed, the preferred minimum temperature and preferred maximum temperature (estimated from a model) were recorded (Froese 2020).
First, load the necessary packages for this analysis:
#this clears out the environment
rm(list = ls())
#load necessary data packages
library(tidyverse)
library(here)
library(janitor)
library(stringr)
library(terra)
library(sf)
library(calecopal)
library(tmap)
library(ggplot2)
library(knitr)
library(tibble)
library(kableExtra)
library(htmltools)All functions used in this workflow are outlined and defined in this section, and grouped into three subcategories:
This function can be used to generate a Kable from a tibble or dataframe that is formatted to match the general formatting in this document (Zhu 2024).
The input is:
data: a tibble or dataframe that contains the values to be displayed in the table
col_names: a named vector with the names to be used for the column headers
caption_text: the title to be used for the table
The output is:
kable_output: a kable generated from the data input with the col_names as headers
this kable is saved with the name “data + _kable”
this output is also saved to the tables subfolder in the figures folder
#write a function to generate a standard kable format
standard_kable <- function(data, col_names, caption_text) {
#extract the name of the data to save as the kable filename
filename <- deparse(substitute(data))
#write the standard formatting for the kable outputs
kable_output <- data %>%
kable("html",
col.names = col_names,
caption = htmltools::tags$div(style = "text-align: center; font-size: 20px;",
htmltools::tags$strong(caption_text))) %>%
kable_styling(full_width = FALSE, font_size = 14) %>%
row_spec(row = 0, bold = TRUE) %>%
kable_classic(html_font = "Times New Roman")
#define the file path
kable_file <- file.path(tables_folder, paste0(filename, "_kable.html"))
#save the kable
save_kable(kable_output, file = kable_file)
#print the kable
return(kable_output)
}This function can be used to generate a map tmap from a SpatRaster (only compatible with terra rasters) (Tennekes 2018; Hijmans 2024). In this workflow, this function is used to generate a map with the “optimal” and “not optimal” land classifications for each EEZ region for each potential finfish species.
Colors:
The input is:
The output is:
map: a tmap of the rast_obj generated for the region_name
the static tmap is saved with the name “map_title + _map.png”
the interactive tmap is saved with the name “map_title + _map.html”
these outputs are also saved to the maps subfolder in the figures folder
#assign a variable name for the optimal color
optimal_color <- "tomato3"
#assign a variable name for the not optimal color
not_optimal_color <- "lightgrey"
#write a function to generate a map for the species raster object
generate_maps <- function(rast_list, eez_data, map_title = NULL) {
#create an empty list to hold the interactive and static maps
map_list <- list()
#extract the first raster from the list to be used as the base raster
base_raster <- rast_list[[1]]
#save the raster layer name as base_name
base_name <- names(rast_list)[1]
#if there is a map_title entry, use that for the title, if not, use the base_name
title <- if (!is.null(map_title)) map_title else base_name
#create a base map to be used from the first raster in the list
combined_map <- tm_shape(base_raster) +
tm_raster(palette = c(not_optimal_color, optimal_color),
alpha = 0.5,
labels = c("Not Optimal", "Optimal"),
title = "Suitability") +
tm_shape(eez_data) +
tm_fill(col = "color",
popup.vars = c("EEZ Region" = "region"),
alpha = 0.5,
title = "EEZ Regions") +
tm_layout(main.title = title,
main.title.position = "center",
main.title.size = 1.5,
main.title.fontfamily = "Times New Roman",
main.title.fontface = "bold",
legend.outside = TRUE)
#for every object in the raster list create a map by iterating through the remaining rasters in the list
if (length(rast_list) > 1) {
for (i in 2:length(rast_list)) {
#extract the raster from the list
raster_layer <- rast_list[[i]]
#add the new layer to the map
combined_map <- combined_map +
tm_shape(raster_layer) +
tm_raster(palette = c(not_optimal_color, optimal_color),
alpha = 0.5,
legend.show = FALSE)
}
}
#make the map name using map_title and define file path for the static .png file
map_name_png <- file.path(maps_folder, paste(title, "map.png", sep = "_"))
#save the static tmap as PNG
tmap_save(combined_map, filename = map_name_png)
#create an interactive title
title_int <- htmltools::tags$h1(title,
style = "text-align: center; font-family: Times New Roman; font-size: 24px; font-weight: bold;")
#combine the interactive title and map into an interactive map object
map_int <- htmltools::browsable(htmltools::tagList(title_int,
tmap_leaflet(combined_map)))
#make the HTML file name and define file path
map_name_html <- file.path(maps_folder, paste0(title, "_map.html"))
#save the interactive map as an HTML file
htmltools::save_html(map_int, file = map_name_html)
#add the interactive map to the map_list
map_list$interactive_map <- map_int
#add the static map to the map_list
map_list$static_map <- map_name_png
#return the file path for the combined interactive map and the static file
return(map_list)
}This function can be used to generate a heatmap from the final region rank dataframe generated from the rank_regions function. In this workflow, this function is used to create a heatmap which visualizes the regions that are the most suitable for aquaculture farms for both oysters and the 26 potential finfish species.
Colors:
The input is:
The output is:
heatmap: a heatmap
the heatmap is saved with the name “filename + _heatmap.png”
this output is also saved to the heatmaps subfolder in the figures folder
#assign a variable name for the optimal color
optimal_color <- "tomato3"
#write a function to generate a heatmap for the potential species and their ranking per region
generate_heatmap <- function(rank_df) {
#extract the name of the data to save as the filename
filename <- deparse(substitute(rank_df))
#create a heatmap from the rank_df
heatmap <- ggplot(rank_df, aes(x = species, y = region, fill = percent_cover)) +
geom_tile(color = "black", size = 0.5) +
geom_text(aes(label = round(percent_cover, 1)), #show the percent cover to 1 decimal
color = "black", size = 2) +
scale_fill_gradient2(midpoint = 0, #make the midpoint 0 so we can make the 0 percent_cover values the "azure" color
mid = "azure",
high = optimal_color,
name = "Percent\nCover",
breaks = seq(0, 5, by = 1),
limits = c(0, 5)) +
scale_y_discrete(labels = c("washington" = "Washington",
"oregon" = "Oregon",
"northern_california" = "Northern California",
"central_california" = "Central California",
"southern_california" = "Southern California")) +
labs(title = "Species vs Region Rankings",
subtitle = "Heatmap of Percent Optimal Cover",
x = "Species",
y = "Region") +
theme(plot.title = element_text(family = "Times New Roman",
face = "bold",
size = 16,
hjust = 0.5),
plot.subtitle = element_text(family = "Times New Roman",
size = 14,
hjust = 0.5),
axis.title = element_text(family = "Times New Roman",
face = "bold",
size = 14),
axis.text.x = element_text(family = "Times New Roman",
size = 10,
angle = 45,
hjust = 1,
vjust = 1),
axis.text.y = element_text(family = "Times New Roman",
size = 10),
legend.title = element_text(family = "Times New Roman",
face = "bold",
size = 10),
legend.text = element_text(family = "Times New Roman",
size = 10),
legend.position = "right",
legend.box.background = element_rect(color = "black", size = 1))
#define the file path
heatmap_file <- file.path(heatmap_folder, paste0(filename, "_heatmap.png"))
#save the heatmap
ggsave(heatmap, file = heatmap_file)
#return the heatmap
return(heatmap)
}This function can be used to create a tibble with all of the geometric properties (resolution, extent, origin, and coordinate reference system) for a list of spatial objects with the class “SpatRaster” or “sf object” (both terra and sf compatible). This is useful for comparing the spatial properties of the spatial objects in a workflow to ensure compatibility before processing begins.
The input is:
The output is:
#write a function to create a tibble with all the spatial properties for a list of spatial objects
extract_spatial_properties <- function(spatial_file_list) {
#create an empty tibble to hold the output
spatial_properties <- tibble(name = character(),
resolution = character(),
extent = character(),
origin = character(),
crs = character())
#extract the name of the list
list_name <- deparse(substitute(spatial_file_list))
#write a for loop to extract the properties into the tibble for review
for (file_name in names(spatial_file_list)) {
#extract the current spatial object
obj <- spatial_file_list[[file_name]]
#since we have sf objects and SpatRasters, we must handle them differently:
##if the object is a SpatRaster
if (inherits(obj, "SpatRaster")){
#if the object is a raster stack (has more than one layer):
if (nlyr(obj) > 1) {
##for each layer of the raster stack
for (layer_name in names(obj)) {
#extract the layer
layer <- obj[[layer_name]]
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string with 5 decimals
resolution <- paste(round(res(obj), 5), collapse = ", ")
##extract the extent (raster size) as a string with 2 decimals
extent <- paste(round(ext(obj), 2), collapse = ", ")
##extract the origin (0,0) as a string with 7 decimals
origin <- paste(round(origin(obj), 7), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(obj)
#extract the EPSG ID part from the CRS string to make the CRS easier to read
crs_id <- sub('.*ID\\["EPSG",(\\d+)\\].*', 'EPSG:\\1', as.character(crs))
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = layer_name,
resolution = resolution,
extent = extent,
origin = origin,
crs = crs_id)
}
#if the object is a single raster:
} else {
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string with 5 decimals
resolution <- paste(round(res(obj), 5), collapse = ", ")
##extract the extent (raster size) as a string with 2 decimals
extent <- paste(round(ext(obj), 2), collapse = ", ")
##extract the origin (0,0) as a string with 7 decimals
origin <- paste(round(origin(obj), 7), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(obj)
#extract the EPSG ID part from the CRS string to make the CRS easier to read
crs_id <- sub('.*ID\\["EPSG",(\\d+)\\].*', 'EPSG:\\1', as.character(crs))
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = file_name,
resolution = resolution,
extent = extent,
origin = origin,
crs = crs_id)
}
##if the object is a sf object
} else if (inherits(obj, "sf")) {
#extract the spatial properties for this layer:
##extract the EPSG ID part from the CRS string to make the CRS easier to read
crs_id <- gsub('.*ID\\[\\\"(EPSG)\\\",(\\d+)\\].*', '\\1:\\2', as.character(crs))
#add the new information to the tibble for raster properties
spatial_properties <- spatial_properties %>%
add_row(name = file_name,
resolution = NA,
extent = NA,
origin = NA,
crs = crs_id)
}
}
#remove duplicate rows based on the "name" column that are generated for sf objects (inherently sf and dataframes)
spatial_properties <- spatial_properties %>%
distinct(name, .keep_all = TRUE)
#create a tibble name to save
tibble_name <- paste0(list_name, "_tibble")
#rename the tibble with the list name and save to the global environment
return(assign(tibble_name, spatial_properties, envir = .GlobalEnv))
}This function can be used to verify that all the spatial objects in a list (both terra and sf compatible) have a pre-defined coordinate reference system (crs). This is necessary to ensure before any map algebra can be performed.
The input is:
The output is:
#write a function to verify that all the spatial objects have the same crs
check_crs <- function(obj, target_crs) {
#if the object is a "SpatRaster"
if (inherits(obj, "SpatRaster")) {
#returns a logical (TRUE/FALSE) if the crs matches the target_crs (TRUE) or not (FALSE)
return(terra::crs(obj) == target_crs)
#if the object is an "sf object"
} else if (inherits(obj, "sf")) {
#returns a logical (TRUE/FALSE) if the crs matches the target_crs (TRUE) or not (FALSE)
return(st_crs(obj)$wkt == target_crs) ##need to use wkt since terra and sf print crs differently
}
}This function can be used to check that the geometric properties (resolution, extent, origin, and crs) of two spatial objects (compatible with terra and sf) match. This function is a modified version of the extract_spatial_properties function.
The input is:
The output is:
#write a modified version of the extract_spatial_properties function to check the spatial properties
check_spatial_properties <- function(spatial_file_list, ref_crs) {
#create an empty tibble to hold the output
spatial_properties_check <- tibble(name = character(),
resolution_match = character(),
extent_match = character(),
origin_match = character(),
crs_match = character())
#extract the name of the list
list_name <- deparse(substitute(spatial_file_list))
#we will use the first file in the list as the reference object for verification
ref_obj <- spatial_file_list[[1]]
#now identify the spatial properties of the first object
##identify the reference resolution
ref_resolution <- paste(res(ref_obj), collapse = ", ")
##identify the reference extent
ref_extent <- paste(ext(ref_obj), collapse = ", ")
##identify the reference origin
ref_origin <- paste(origin(ref_obj), collapse = ", ")
##identify the reference crs
ref_crs <- ref_crs
#write a for loop to extract the properties into the tibble for review
for (file_name in names(spatial_file_list)) {
#extract the current spatial object
obj <- spatial_file_list[[file_name]]
#since we have sf objects and SpatRasters, we must handle them differently:
##if the object is a SpatRaster
if (inherits(obj, "SpatRaster")){
#if the object is a raster stack (has more than one layer):
if (nlyr(obj) > 1) {
##for each layer of the raster stack
for (layer_name in names(obj)) {
#extract the layer
layer <- obj[[layer_name]]
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string
resolution <- paste(res(obj), collapse = ", ")
##extract the extent (raster size) as a string
extent <- paste(ext(obj), collapse = ", ")
##extract the origin (0,0) as a string
origin <- paste(origin(obj), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(obj)
#check if resolution, extent, origin, and CRS match the reference raster
resolution_match <- ifelse(resolution == ref_resolution, "MATCH", "DOES NOT MATCH")
extent_match <- ifelse(extent == ref_extent, "MATCH", "DOES NOT MATCH")
origin_match <- ifelse(origin == ref_origin, "MATCH", "DOES NOT MATCH")
crs_match <- ifelse(crs == ref_crs, "MATCH", "DOES NOT MATCH")
#add the new information to the tibble for raster properties
spatial_properties_check <- spatial_properties_check %>%
add_row(name = layer_name,
resolution_match = resolution_match,
extent_match = extent_match,
origin_match = origin_match,
crs_match = crs_match)
}
#if the object is a single raster:
} else {
#extract the spatial properties for this layer:
##extract the resolution (cell size) as a string
resolution <- paste(res(obj), collapse = ", ")
##extract the extent (raster size) as a string
extent <- paste(ext(obj), collapse = ", ")
##extract the origin (0,0) as a string
origin <- paste(origin(obj), collapse = ", ")
##extract the crs from the current object
crs <- terra::crs(obj)
#check if resolution, extent, origin, and CRS match the reference raster
resolution_match <- ifelse(resolution == ref_resolution, "MATCH", "DOES NOT MATCH")
extent_match <- ifelse(extent == ref_extent, "MATCH", "DOES NOT MATCH")
origin_match <- ifelse(origin == ref_origin, "MATCH", "DOES NOT MATCH")
crs_match <- ifelse(crs == ref_crs, "MATCH", "DOES NOT MATCH")
#add the new information to the tibble for raster properties
spatial_properties_check <- spatial_properties_check %>%
add_row(name = file_name,
resolution_match = resolution_match,
extent_match = extent_match,
origin_match = origin_match,
crs_match = crs_match)
}
##if the object is a sf object
} else if (inherits(obj, "sf")) {
#extract the spatial properties for this layer:
#extract the crs from the current object
crs <- sf::st_crs(obj)$wkt
#check that the CRS matches the ref_crs
crs_match <- ifelse(crs == ref_crs, "MATCH", "DOES NOT MATCH")
#add the new information to the tibble for raster properties
spatial_properties_check <- spatial_properties_check %>%
add_row(name = file_name,
resolution_match = NA,
extent_match = NA,
origin_match = NA,
crs_match = crs_match)
}
}
#remove duplicate rows based on the "name" column that are generated for sf objects (inherently sf and dataframes)
spatial_properties_check <- spatial_properties_check %>%
distinct(name, .keep_all = TRUE)
#create a tibble name to save
tibble_name <- paste0(list_name, "_tibble")
#rename the tibble with the list name and save to the global environment
return(assign(tibble_name, spatial_properties_check, envir = .GlobalEnv))
}This function can be used to reclassify the geometries in a SpatRaster (only compatible with terra) based on an optimal parameter range. It generates a new SpatRaster that has three classifications:
The input is:
The output is:
#write a function to reclassify the geometries in a raster based on the optimal parameter range
suitable_locations <- function(raster_obj, min_parameter, max_paramter) {
#duplicate the raster that will be reclassified
mask <- raster_obj
#create a reclassification matrix to reclassify the raster based on the optimal parameter range
rcl <- matrix(c(-Inf, min_parameter, 1, #below the minimum value ("below_range")
min_parameter, max_paramter, 2, #in the optimal range ("optimal")
max_paramter, Inf, 3), #above the maximum value ("above_range")
ncol = 3, byrow = TRUE)
#use the reclassification matrix to reclassify the duplicated raster with the rcl
mask <- classify(mask, rcl = rcl)
#assign the group to a classification (not optimal or optimal)
levels(mask) <- data.frame(id = c(1, 2, 3),
cats = c("below_range", "optimal", "above_range"))
#rename the mask to have the name of the raster_obj + _mask
mask_name <- paste(names(raster_obj), "mask", sep = "_")
names(mask) <- mask_name
#return reclassified raster object
return(mask)
}This function is used to generate a new raster object from two optimized SpatRaster (only compatible with terra). This function is formated to work with SpatRasters that were generated from the previously defined suitable_locations function. It generates a new SpatRaster that has two classifications based on the optimal parameter range of both SpatRasters (raster_1 and raster_2):
The input is:
The output is:
#write a function to generate a new raster object from two optimized rasters (this function is formatted to work with rasters generated from the suitable_locations locations defined above)
optimal_locations <- function(raster_1, raster_2){
#duplicate raster_1 to serve as the combined raster
combined <- raster_1
#first assign any locations with missing data from either raster as "not optimal" group 1
combined[raster_1 == 2 | raster_2 == 2] <- 1
#assign locations where BOTH parameters are "optimal" to group 2
combined[raster_1 == 2 & raster_2 == 2] <- 2
#assign locations that are not identified as optimal in both rasters as "not optimal" to group 1
combined[!(combined[] == 2)] <- 1
#assign the categories a name (optimal or not optimal)
levels(combined) <- data.frame(id = c(1, 2),
cats = c("not_optimal", "optimal"))
#return reclassified raster object
return(combined)
}This function is used to rasterize the individual geometries in a sf object based on the name of a column in the dataframe. In this workflow, this was used to rasterize the region geometries from the EEZ dataframe.
The input is:
The output is:
#write a function to rasterize individual geometries based on the column that has the corresponding name
rasterize_geometries <- function(df, template_raster, name_column) {
#create a list of names the name column
name_list <- unique(df[[name_column]])
#create an empty list to store the individual layers
rast_list <- list()
#create a duplicated empty raster from the template raster to be the raster stack
rast_stack <- rast(template_raster, nlyrs = 0)
#create an empty list to hold the bboxes
bbox_list <- list()
#write a for loop to extract the geometry for each id and save it as a raster in the list
for(unique_name in name_list) {
#filter the df to the geometry that corresponds to the current name
geometry <- df[df[[name_column]] == unique_name, ]
#create a duplicated empty raster from the template raster to be the raster layer
rast_layer <- template_raster
#rasterize the extracted geometry
geom_rast <- rasterize(geometry, rast_layer)
#rename the raster with the name
names(geom_rast) <- unique_name
#add the new raster to the list
rast_list[[unique_name]] <- geom_rast
#create a boundary box from the geometry
bbox <- st_bbox(geometry)
#rename the bbox with the unique_name + _bbox
bbox_name <- paste(names(unique_name), "bbox", sep = "_")
names(bbox) <- bbox_name
#add the new bbox to the list
bbox_list[[unique_name]] <- bbox
#add the new raster as a layer to the raster stack
rast_stack <- c(rast_stack, geom_rast)
}
#return a list with the individual layers, the raster stack, and the bbox_list
return(list(rast_list = rast_list,
rast_stack = rast_stack,
bbox_list = bbox_list))
}This function is used to create a mask of a SpatRaster (only compatible with terra) for each layer in a SpatRaster stack. This function is a modified version of the rasterize_geometries function. In this workflow, this function is used to create a masked SpatRaster for each region in the EEZ with the input raster (raster_obj_input) being the SpatRaster generated for each potential species.
The input is:
The output is:
#write a function to create a list with a raster stack, rasters, and a map for each region
extract_regions <- function(raster_obj_input, rast_stack_input) {
#create an empty raster to store the output of the for loop
rast_list <- list()
#create a duplicate of the raster object to be the masked layer
rast_obj <- raster_obj_input
#create a duplicated empty raster from the raster object to be the raster stack
rast_stack <- rast(rast_stack_input, nlyrs = 0)
for (region_name in names(rast_stack_input)) {
#extract the region layer from the raster stack
region_layer <- rast_stack_input[[region_name]]
#use the region_layer to mask the rast_obj
region_layer_mask <- terra::mask(rast_obj, region_layer)
#rename the raster with the region_name
names(region_layer_mask) <- region_name
#add the new raster to the list
rast_list[[region_name]] <- region_layer_mask
#add the new raster as a layer to the raster stack
rast_stack <- c(rast_stack, region_layer_mask)
}
#return a list with the individual layers and the raster stack
return(list(rast_list = rast_list,
rast_stack = rast_stack))
}This function is used to create a tibble that ranks the regional SpatRasters (only compatible with terra) generated for each of the species using the previously defined functions. The regions are ranked based on the percent of the total area that is classified as optimal for cultivation. In this workflow, this function ranks the regions from 1 (best EEZ region for cultivation) to 5 (worst EEZ region for cultivation) based on a calculated optimal percent_cover.
The input is:
The output is:
#write a function to extract the create a tibble with each species and the rank for their optimal regions
rank_regions <- function(rast_stack){
#extract the name of the data to save as the plots filename
filename <- deparse(substitute(rast_stack))
#create a tibble to hold the output of the for loop
region_rank <- tibble(species = character(),
region = character(),
optimal_area = numeric(),
total_area = numeric(),
percent_cover = numeric(),
rank = numeric())
#for every layer in the rast_stack
for (layer_name in names(rast_stack)){
#extract the region layer from the raster stack
layer <- rast_stack[[layer_name]]
#calculate the area (km^2) for each cell
cell_area <- cellSize(layer, unit = "km")
#calculate the area of all the cells with the "not_optimal" (group 1) classification
not_optimal_area_km2 <- sum(values((layer == 1) * cell_area), na.rm = TRUE)
#calculate the area of all the cells with the "optimal" (group 2) classification
optimal_area_km2 <- sum(values((layer == 2) * cell_area), na.rm = TRUE)
#calculate the total area in the region
total_area_km2 <- optimal_area_km2 + not_optimal_area_km2
#calculate the percent cover of optimal area
percent_cover <- (optimal_area_km2 /total_area_km2) *100
#add the new information to the tibble
region_rank <- region_rank %>%
add_row(species = NA,
region = layer_name,
optimal_area = optimal_area_km2,
total_area = total_area_km2,
percent_cover = percent_cover,
rank = NA)
}
#rank the regions with the highest optimal percent cover being ranked 1
region_rank <- region_rank %>%
mutate(rank = ifelse(percent_cover == 0, NA,
rank(-percent_cover, ties.method = "max"))) #if the percent cover is zero assign the rank NA, and if there is a tie, assign it the max rank to be conservative
#create a tibble name to save
tibble_name <- paste0(filename, "_tibble")
#rename the tibble with the filename and save a table with the rankings for each region and each species to the global environment
return(assign(tibble_name, region_rank, envir = .GlobalEnv))
}This section sets up the filepaths and files for later analysis:
Sea Surface Temperature (SST):
To characterize the average sea surface temperature inside the region of interest, we obtained average annual sea surface temperature data for the years 2008 to 2012.
This data was generated by the NOAA Coral Reef Watch Daily 5km Satellite Sea Surface Temperature (SST) (v3.1) (NOAA Coral Reef Watch (2024)).
Bathymetry:
Economic Exclusion Zone (EEZ):
To define the area of U.S. waters that could potential be regions of interest for a marine aquaculture operation, we obtained a shapefile outlining the economic exclusion zones (EEZ) along the west coast of the U.S.
This data was developed by the Flanders Marine Institute (Flanders Marine Institute (VLIZ) (2024)).
#define the pathway to access the average annual sea surface temperature (sst) files
sst_files <- list.files(here::here("data"),
pattern = "average_annual_sst.*\\.tif$",
full.names = TRUE)
#define the pathway to access the depth file
depth_file <- here::here("data", "depth.tif")
#define the pathway to access the economic exclusion zone (eez) files
eez_file <- here::here("data", "wc_regions_clean.shp")
#define the pathway to access the csv with the potential finfish species data
species_file <- here::here("data", "potential_finfish_species.csv")#define the pathway to save to the maps subfolder in the figures folder
maps_folder <- here::here("figures", "maps")
#define the pathway to save to the tables subfolder in the figures folder
tables_folder <- here::here("figures", "tables")
#define the pathway to save to the heatmap subfolder in the figures folder
heatmap_folder <- here::here("figures", "heatmap")Now that we have defined the file paths, we can read in the necessary spatial data:
#read in the sea surface temperature (sst) rasters
sst_rasters <- terra::rast(sst_files)
#write a for loop to rename each layer to have the format sst_year
for (i in seq_along(names(sst_rasters))) {
#extract the year for the corresponding file, by extracting the last four digits
year <- str_extract(basename(sst_files[i]), "\\d{4}") #extracting the last four digits
#format the name of the current layer to be sst_year and save it as the new name for the layer
names(sst_rasters)[i] <- paste("sst", year, sep = "_")
}#read in the depth raster
depth_raster <- terra::rast(depth_file)
#rename the layer
names(depth_raster) <- "depth"#read in the economic exclusion zone (eez) raster
eez <- sf::st_read(eez_file)
#create a vector with region_id's that correspond to the regions based on the geographic order
region_ids <- c("Washington" = 1,
"Oregon" = 2,
"Northern California" = 3,
"Central California" = 4,
"Southern California" = 5)
#create a region_id column in the eez dataset that assigns the regions a number from north (washington) to south (southern_california)
eez <- eez %>%
mutate(region_id = region_ids[rgn])
#create a version of the dataframe with the capitalized region names still
eez_caps <- eez %>%
clean_names() %>%
rename(region = rgn) %>%
select(region_id, region, area_km2, geometry)
#tidy up the data
eez <- eez %>%
clean_names() %>% #clean up the column names
mutate(region = str_replace_all(rgn, " ", "_") %>%
str_to_lower()) %>% #create a region column with the name in snake_case
select(region, area_km2, geometry) #only keep the necessary columnsNext, we need to examine the geometric properties (resolution, extent, origin, and coordinate reference system (crs)) for each of the spatial objects that we are working with. We will need to transform any spatial objects that have mismatched properties.
We will start by identifying the geometric properties of each of the spatial objects with the extract_spatial_properties function and save these properties in a table with the standard_kable function:
#create a summary list with all the files
spatial_file_list <- list(sst_rasters = sst_rasters,
depth_raster = depth_raster,
eez = eez)
#identify the spatial properties of all the spatial files in the list
spatial_properties <- extract_spatial_properties(spatial_file_list)
#save the column names for the kable
col_names <- c("Name" = "Name",
"Resolution" = "Resolution",
"Extent" = "Extent",
"Origin" = "Origin",
"CRS" = "CRS")
#generate the kable from the tibble output
standard_kable(spatial_properties, col_names, caption = "Geometric Properties Identification Results")| Name | Resolution | Extent | Origin | CRS |
|---|---|---|---|---|
| sst_2008 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2009 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2010 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2011 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| sst_2012 | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| depth_raster | 0.00417, 0.00417 | ext(-132, -114, 29, 50) | 0, 0 | EPSG:4326 |
| eez | NA | NA | NA | EPSG:4326 |
It appears that the coordiante reference system for the bathymetry data (depth_raster) and the economic exclusion zone data (eez) does not match that of the sea surface temperature data (sst_rasters).
Therefore, we will transform the bathymetry data (depth_raster) and the economic exclusion zone data (eez) to match the coordinate reference system (crs = EPSG:9122) of sea surface temperature data (sst_rasters).
#define the target CRS (EPSG:9122) to be used in the transformations
target_crs <- crs(sst_rasters)
#transform the depth_raster to the target crs
depth_raster <- terra::project(depth_raster, target_crs)
#transform the eez to the target crs
eez <- st_transform(eez, crs = target_crs)
#create a summary list with all the files
verify_crs_list <- list(sst_rasters = sst_rasters,
depth_raster = depth_raster,
eez = eez)Now verify that all of the spatial files have the same coordinate reference system with the check_crs function before continuing with analysis:
#check that all the files have the same crs (target_crs)
all_crs_match <- sapply(verify_crs_list, check_crs, target_crs = target_crs)
#verify that all the files have the same crs
if (all(all_crs_match)) {
#if all the files have the same crs (all TRUE) then print
print("All spatial files have the target CRS")
} else {
#if all the files do NOT have the same crs (at least one FALSE) then print
print("The following files do not match the target CRS:")
print(names(spatial_file_list)[!all_crs_match])
}[1] "All spatial files have the target CRS"Since all the spatial files have the same coordinate reference system, we can move to the analysis process.
This section identifies locations that would be suitable for only bivalve farming:
To start, we will ensure that the sea surface temperature (SST) and the bathymetry data are compatible to process together.
First, we need to calculate the average sea surface temperature by taking the mean across all of the rasters (from the years 2008 to 2012), and then convert this temperature from Kelvin to Celsius.
#find the mean sea surface temperature and store it as a single layer raster
sst_mean <- terra::mean(sst_rasters, na.rm = TRUE)
#now convert the temperature to Celcius by subtracting 273.15 to convert from Kelvin
sst_mean_celsius <- sst_mean - 273.15
#rename the layer for the celcisus raster
names(sst_mean_celsius) <- "mean_sst_c"Now we need to verify that the rasters are compatible to be combined. We will start by identifying the geometric properties of with the extract_spatial_properties function and save these properties in a table with the standard_kable function:
#create a summary list with all the files
sst_bath_file_list <- list(sst_mean_celsius = sst_mean_celsius,
depth_raster = depth_raster)
#identify the spatial properties of all the files in the list
properties_sst_bath <- extract_spatial_properties(sst_bath_file_list)
#generate the kable from the tibble output
standard_kable(properties_sst_bath, col_names, caption = "Geometric Properties Identification Results")| Name | Resolution | Extent | Origin | CRS |
|---|---|---|---|---|
| sst_mean_celsius | 0.04166, 0.04166 | ext(-131.98, -114.99, 29.99, 49.99) | -7.6e-06, -3.8e-06 | EPSG:9122 |
| depth_raster | 0.00417, 0.00417 | ext(-132, -114, 29, 50) | 0, 0 | EPSG:9122 |
It appears that the rasters have the same coordinate reference systems (since we transformed them to match above), but differ in their resolution, extent, and origin.
To ensure compatible geometric properties we will transform the data so both raster have the same resolution, extent, origin, and coordinate reference system. Since the bathymetry data has a larger extent than the sea surface temperature, we can crop this raster to match the extent and origin of the temperature raster.
#crop the bathymetry raster to match the extent/origin of the sst raster
depth_raster_crop <- terra::crop(depth_raster, sst_mean_celsius)Since the bathymetry data has a finer resolution than the sea surface temperature (SST) data, we will coarsen the bathymetry data by resampling the data using the nearest neighbor approach.
#resample the cropped bathymetry data using the nearest neighbor approach
depth_resampled <- terra::resample(depth_raster_crop, sst_mean_celsius, method = "near")Now we will verify that the transformations were successful and the resulting rasters have compatible geometric properties.
Now we need to verify that the rasters are compatible to be combined by ensuring the geometric properties have been successfully transformed with the check_spatial_properties function and save these properties in a table with the standard_kable function:
#create a summary list with all the files
sst_bath_transform_file_list <- list(sst_mean_celsius = sst_mean_celsius,
depth_resampled = depth_resampled)
#identify the spatial properties of all the files in the list
properties_sst_bath_transform <- check_spatial_properties(sst_bath_transform_file_list, ref_crs = target_crs)
#generate the kable from the tibble output
standard_kable(properties_sst_bath_transform, col_names, caption = "Geometric Properties Identification Results")| Name | Resolution | Extent | Origin | CRS |
|---|---|---|---|---|
| sst_mean_celsius | MATCH | MATCH | MATCH | MATCH |
| depth_resampled | MATCH | MATCH | MATCH | MATCH |
Since the transformations were successful, we can move forward with identifying suitable locations for oyster cultivation along the west coast.
Before we can identify potential finfish species that will be suitable options for marine aquaculture, we need to determine which locations are within the scope of the temperature and depth parameters for oysters.
To start, we will classify locations as “optimal” or “not optimal” for oyster cultivation based on the average sea surface temperature calculated above. This is done by creating a mask of the average sea surface temperature data (mean_sst_celcius) and assigning appropriate categorical levels using the suitable_locations function:
#define the minimum temperature
min_temp <- 11
#define the maximum temperature
max_temp <- 30
#generate a masked version of the temperature raster to identify optimal and not optimal locations
sst_mask <- suitable_locations(sst_mean_celsius, min_temp, max_temp)Now we classified locations as “optimal”, “not optimal”, or “land” for oyster cultivation based on the bathymetry data processed above. This is done by creating a mask of the bathymetry data (depth_resampled) and assigning appropriate categorical levels using the suitable_locations function:
#define the "minimum" depth (meters below sea level)
min_depth <- -70
#define the "maximum" depth (sea surface)
max_depth <- 0
#generate a masked version of the bathymetry raster to identify optimal and not optimal locations
depth_mask <- suitable_locations(depth_resampled, min_depth, max_depth)
#we will now rename the "not optimal" group 3 to be land (elevations above sea level)
levels(depth_mask) <- data.frame(id = c(1, 2, 3),
cats = c("not_optimal", "optimal", "land"))Once the physical condition rasters have been reclassified, we can combine the reclassified rasters (sst_mask and depth_mask) into a single raster that contains the optimal locations for oyster cultivation based on both parameters using the optimal_locations function:
#create a raster for the optimal geographic range for oyster cultivation based on the temperature and depth parameters
oyster_range <- optimal_locations(sst_mask, depth_mask)After generating a raster with the optimal location for oyster cultivation, we can overlay the economic exclusion zones (eez) data. Before combining this information, we must rasterize the eez data to make it compatible with the other SpatRaster objects. We will create an individual raster object for each region (Southern California, Central California, Northern California, Oregon, and Washington) of the economic exclusion zones as well as a raster stack which contains each region as a layer using the rasterize_geometries function:
#rasterize the regions and create a raster stack with all of the regions
eez_results <- rasterize_geometries(eez, oyster_range, "region")
#save the list of region rasters
eez_regions <- eez_results$rast_list
#save the raster stack
eez_stack <- eez_results$rast_stack
#save the list of boundary boxes for each region
eez_bbox_list <- eez_results$bbox_listNow we can visualize the suitable locations for oyster cultivation in each region of the economic exclusion zone (Washington, Oregon, Northern California, Central California, Southern California) along the west coast of the U.S.
Let’s start by making a color palette to visualize each of the regions within the west coast EEZ using these colors.
#generate a palette from the calecopal bigsur color palette
bigsur_palette <- cal_palette(name = "bigsur")
#create a vector of selected hex codes from the bigsur_palette
region_palette <- c("#E4DECE", "#ECBD95", "#9BB1BB", "#346575", "#0B4221")
#assign the colors in the palette to region names
region_names <- c("Washington",
"Oregon",
"Northern California",
"Central California",
"Southern California")
#assign the colors in the palette to the list of names in region_names
region_palette_named <- setNames(region_palette, region_names)
#create a copy of the eez dataframe with a color that corresponds to each region in the region_palette_named
eez_color <- eez_caps %>%
mutate(color = region_palette_named[region])And then we can use the generate_maps function to create an interactive map to visualize the optimal locations (optimal and not optimal) for oyster cultivation within each EEZ region:
#use the generate maps function to create an interactive map of the suitable oyster cultivation locations
oyster_map_results <- generate_maps(oyster_range,
eez_data = eez_color,
map_title = "Suitable Oyster Cultivation Locations and EEZ Regions")
#extract the interactive map from the oyster_map_results
interactive_map <- oyster_map_results$interactive_map
#print the interactive map
interactive_mapThis section identifies locations that would be suitable for BOTH bivalve and finfish aquaculture:
Start by loading in the csv with the metadata for the 26 potential finfish species identified in the background section:
#read in the csv for the potential species
potential_species <- read_csv(species_file)
#subset the data to only the species common name and scientific name
potential_species_tbl <- potential_species %>%
select(1:2) %>%
mutate(number = paste0(row_number())) %>% #add a number based on the row number in the df
select(3,1,2)
#select only the first 13 rows of the table
potential_species_tbl_1 <- potential_species_tbl %>%
slice(1:13)
#select only the last 13 rows of the table
potential_species_tbl_2 <- potential_species_tbl %>%
slice(14:26)
#combine the two tables
potential_species_comb <- cbind(potential_species_tbl_1, potential_species_tbl_2)
#tidy up the data
potential_species <- potential_species %>%
clean_names() %>%
mutate(common_name = str_replace_all(common_name, " ", "_") %>% str_to_lower(),
scientific_name = str_replace_all(scientific_name, " ", "_") %>% str_to_lower())Then print out a table with a list of the potential species being considered and save this output in a table with the standard_kable function:
#save the column names for the kable
col_names_p <- c("Number" = "Number",
"Common Name" = "Common Name",
"Latin Name" = "Latin Name",
"Number" = "Number",
"Common Name" = "Common Name",
"Latin Name" = "Latin Name")
#generate the kable from the tibble output
standard_kable(potential_species_comb, col_names_p, caption = "Potential Finfish Species") %>%
column_spec(column = c(3, 6), italic = TRUE) #reformat for the resized table output | Number | Common Name | Latin Name | Number | Common Name | Latin Name |
|---|---|---|---|---|---|
| 1 | Arrowtooth Flounder | Atheresthes stomias | 14 | Pacific Ocean Perch | Sebastes alutus |
| 2 | Bocaccio | Sebastes paucispinis | 15 | Pacific Sardine | Sardinops sagax |
| 3 | Canary Rockfish | Sebastes pinniger | 16 | Pacific Whiting | Merluccius productus |
| 4 | Chinook Salmon | Oncorhynchus tshawytscha | 17 | Petrale Sole | Eopsetta jordani |
| 5 | Chum Salmon | Oncorhynchus keta | 18 | Pink Salmon | Oncorhynchus gorbuscha |
| 6 | Coho Salmon | Oncorhynchus kisutch | 19 | Rex Sole | Glyptocephalus zachirus |
| 7 | Dover Sole | Microstomus pacificus | 20 | Pacific Rock Sole | Lepidopsetta billineta |
| 8 | English Sole | Parophrys vetulus | 21 | Northern Rock Sole | Lepidopsetta polyxystra |
| 9 | Flathead Sole | Hippoglossoides elassodon | 22 | Sablefish | Anoplopoma fimbria |
| 10 | Lingcod | Ophiodon elongatus | 23 | Shortspine Thornyhead | Sebastolobus alascanus |
| 11 | Northern Anchovy | Engraulis mordax | 24 | Sockeye Salmon | Oncorhynchus nerka |
| 12 | Pacific Cod | Gadus macrocephalus | 25 | Widow Rockfish | Sebastes entomelas |
| 13 | Pacific Halibut | Hippoglossus stenolepis | 26 | Yellowtail Rockfish | Sebastes flavidus |
Now we will use the functions defined above to evaluate the compatibility of each potential finfish species in each EEZ region. For each species, we will generate a SpatRaster for each EEZ region whose area has been classified as “optimal” or “not optimal” for the cultivation of that species and oysters. We will then stack these SpatRasters to generate a SpatRaster stack for each species containing a layer for each region. Additionally we will generate a tibble which ranks each EEZ region for each potential finfish species based on the percent of optimal area in the region. Finally, we will also create a map for each region layer generated for each potential species. The functions are used in the below order:
Output: This workflow will generate three outputs:
#create an empty list to hold the output of the for loop
species_rast_list <- list()
#create a dataframe column with the scientific names from the potential_species_tbl
map_name_df <- potential_species_tbl[, 3]
#convert the map_name_df to a list
map_name_list <- as.list(map_name_df)
#bind the map_name_list on to the potential_species dataframe
potential_species_names <- cbind(potential_species, map_name_df)
#rename the new column to be tidyverse friendly
potential_species_names <- potential_species_names %>%
rename(species_name_caps = 9)
#create a tibble to hold the output of the for loop
region_rank <- tibble(species = character(),
region = character(),
optimal_area = numeric(),
total_area = numeric(),
percent_cover = numeric(),
rank = numeric())
#create a list to hold the species maps
species_map <- list()
#for every species in the potential species dataframe
for(i in 1:nrow(potential_species_names)) {
#extract the properly formatted species name from the map_name_list
species_name_caps <- potential_species_names$species_name_caps[i]
#save the species name as the scientific name
species_name <- potential_species_names$scientific_name[i]
#extract the minimum temperature for this species
min_temp <- potential_species_names$min_temp[i]
#extract the maximum temperature for this species
max_temp <- potential_species_names$max_temp[i]
#since the thermal range of tolerance was predicted from a model for some species, we have classified this in other columns (preferred_min_temp and preferred_max_temp) and left the min_temp and max_temp column empty
##so if the min_temp is NA we will use the preferred_min_temp instead
min_temp <- ifelse(is.na(min_temp), #if the min_temp is NA
potential_species_names$preferred_min_temp[i], #use the preferred_min_temp
min_temp) #if is.na(min_temp) = FALSE then use the min_temp already extracted
##or if the max_temp is NA we will use the preferred_max_temp instead
max_temp <- ifelse(is.na(max_temp), #if the max_temp is NA
potential_species_names$preferred_max_temp[i], #use the preferred_max_temp
max_temp) #if is.na(max_temp) = FALSE then use the max_temp already extracted
#the min_depth and max_depth in the dataframe are recorded as negative numbers to represent meters below sea level to align with the format in the bathymetry data which has a vertical datum of mean sea level
##however the function wants the min_parameter value which would correspond to the max_depth, so extract the maximum depth for this species
min_depth <- potential_species_names$max_depth[i]
##the function wants the max_parameter value which would correspond to the min_depth, so extract the minimum depth for this species
max_depth <- potential_species_names$min_depth[i]
#FUNCTION ONE A: generate a masked version of the sea surface temperature data using the suitable_locations function
temp_rast <- suitable_locations(sst_mean_celsius, min_temp, max_temp)
#FUNCTION ONE B:generate a masked version of the bathymetry using the suitable_locations function
depth_rast <- suitable_locations(depth_resampled, min_depth, max_depth)
#FUNCTION TWO A: determine the optimal range for this species using the optimal_locations function
optimal_rast <- optimal_locations(temp_rast, depth_rast)
#FUNCTION TWO B: determine the optimal range for this species AND the oysters using the optimal_locations function
potential_locations <- optimal_locations(oyster_range, optimal_rast)
#FUNCTION THREE: generate a raster stack with each layer in the stack being a region in the economic exclusion zone
results_list <- extract_regions(potential_locations, eez_stack)
#extract the species raster stack from the results list
species_stack <- results_list$rast_stack
#extract the species raster list from the results list
species_list <- results_list$rast_list
#FUNCTION FOUR: rank the regions for each species
species_rank <- rank_regions(species_stack)
#add the species name to the tibble
species_rank <- species_rank %>%
mutate(species = species_name)
#add the new species_rank tibble to the full tibble
region_rank <- bind_rows(region_rank, species_rank)
#FUNCTION FIVE: generate a static and an interactive map
map_results <- generate_maps(rast_list = species_list,
eez_data = eez_color,
map_title = species_name_caps)
#add the map to the species_map list
species_map[[species_name]] <- map_results$interactive_map
#save the species results list to the species rast list
species_rast_list[[species_name]] <- results_list
}Now let’s visualize some of the results from our analysis.
This section summarizes the results of the above analysis through a variety of visualizations:
First, we can create a heatmap with the generate_heatmap function to determine which regions are suitable most for aquaculture farms for both oysters and the 26 potential finfish species. Here we can examine the optimal percent cover calculated for each species in each region with locations with ZERO optimal area in light blue and locations with HIGH percent cover of optimal area shown in a gradient of red.
#make sure the regions are in order from south to north
region_rank$region <- factor(region_rank$region,
levels = c("southern_california",
"central_california",
"northern_california",
"oregon",
"washington"))
#create a heatmap of the potential species and their ranking per region
generate_heatmap(region_rank)
Then we can determine which regions will be the most optimal for creating farms for oyster and a finfish species and save these properties in a table with the standard_kable function:
#subset the region_rank dataframe to only include the regions that are the best for cultivation
top_regions <- region_rank %>%
group_by(species) %>%
filter(rank == 1) %>%
ungroup()
#create a list of column names to be used in the kable
col_names_r <- c("Species" = "Species",
"Region" = "Region",
"Optimal Area" = "Optimal Area",
"Total Area" = "Total Area",
"Percent Cover" = "Percent Cover",
"Rank" = "Rank")
#generate the kable from the tibble output
standard_kable(top_regions, col_names_r, caption = "Top Ranked Regions for Each Species")| Species | Region | Optimal Area | Total Area | Percent Cover | Rank |
|---|---|---|---|---|---|
| sebastes_pinniger | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| oncorhynchus_tshawytscha | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| oncorhynchus_keta | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| oncorhynchus_kisutch | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| hippoglossoides_elassodon | washington | 2842.345 | 67478.12 | 4.212248 | 1 |
| engraulis_mordax | washington | 3209.864 | 67478.12 | 4.756897 | 1 |
| sardinops_sagax | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| merluccius_productus | washington | 3180.263 | 67478.12 | 4.713028 | 1 |
| oncorhynchus_gorbuscha | washington | 3224.738 | 67478.12 | 4.778938 | 1 |
| sebastes_flavidus | washington | 2842.345 | 67478.12 | 4.212248 | 1 |
Finally, we can examine one of the regional maps generated for a potential finfish species. We will visualize the maps generated for the Sebastes pinniger (Canary Rockfish) as an example of the maps created for all of the potential finfish species. We will color the regions similarly with the optimal locations (optimal and not optimal) used to visualize areas for cultivation within each EEZ region (Washington, Oregon, Northern California, Central California, Southern California):
#examine the maps generated for the sebastes_pinniger species
tmap_mode("view")
#extract the species map from the map list
sebastes_pinniger_map <- species_map[["sebastes_pinniger"]]
#print the map
sebastes_pinniger_mapThese results show that for all of the potential finfish species chosen, the Washington EEZ is the most suitable for paired oyster aquaculture. However, further parameters for each species should be examined more closely for potential areas with the Washington region before moving forward with development.