Directory and doc rules

knitr::opts_chunk$set(
  echo = TRUE,         # Display code chunks
  eval = TRUE,         # Evaluate code chunks
  warning = FALSE,     # Hide warnings
  message = FALSE,     # Hide messages
  fig.width = 20,       # Set plot width in inches
  fig.height = 9,      # Set plot height in inches
  fig.align = "center" # Align plots to the center
)

Load packages

#install.packages(c("ggplot2", "sf", "viridis", "rnaturalearth", "rnaturalearthdata"))
library(tidyr)
#library(tidyverse)
library(ggplot2)
library(vegan)
library(sf)
library(viridis)
library(rnaturalearth)
library(rnaturalearthdata)

Load data

Weight = mg
Length, width, height = mm
p450, SOD = activity/ (mg/protein)
Condition factor, economic factor = unitless

getwd()
## [1] "/Users/cmantegna/Documents/WDFWmussels/code"
#data has all sites, coordinates, p450, sod, condition factor, economic factor data
data<- read.csv("/Users/cmantegna/Documents/WDFWmussels/data/biomarkerfull.csv")

#alldata has the site names, biomarkers, condition factor, average thickness and analyte data - each row is an individual sample
alldata<- read.csv("/Users/cmantegna/Documents/WDFWmussels/data/alldata.csv")

fix zero’s in the data frame and alldata frame

# Data contains 0's and must be adjusted in this order to preserve all usable data.

#sod
#replace any SOD values at or below 0 with half of the lower detection limit of .005 (.005*.5). Lower detection limit determined by assay protocol by the manufacturer, Cayman.
data$SOD[data$SOD <= 0] <- 0.0025
alldata$SOD[alldata$SOD <= 0] <- 0.0025
#p450
#remove any p450 values that are 0 - those are true 0's not non-detectable. I am replacing with na so I don't lose the entire row of data, including the SOD values.
data$p450[data$p450 <= 0] <- NA
alldata$p450[alldata$p450 <= 0] <- NA

#write.csv(alldata, "/Users/cmantegna/Documents/WDFWmussels/data/alldata.csv")

Map of Washington State

world <- ne_states(country = "united states of america", returnclass = "sf")
washington_map <- world[world$name == "Washington", ]

p450 map- WA state

pmap<- ggplot() + 
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point(data = alldata, aes(x = longitude, y = latitude, color = p450), size = 2) +
  scale_color_viridis(option = "C", name = "p450") +
  theme_minimal() +
  labs(title = "Washington State - p450 Data")

print(pmap)

ggsave(plot=pmap, filename="/Users/cmantegna/Documents/WDFWmussels/output/p450map.png", width=15, height=8)

p450 map - Puget Sound

#zoom into puget sound region & note the legend, lighter colors are higher values
xlim <- c(-124, -122)  # longitude bounds
ylim <- c(47, 49)  # latitude bounds

sodpugetsound<- ggplot() + 
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point(data = alldata, aes(x = longitude, y = latitude, color = p450), size = 3) +
  scale_color_viridis(option = "C", name = "SOD") +
  coord_sf(xlim = xlim, ylim = ylim, expand = FALSE)+
  theme_minimal() +
  labs(title = "Washington State - SOD Data")

print(sodpugetsound)

p450 map by regions

fix the coordinates and maps, not a priority 4/21

library(ggplot2)
library(dplyr)
library(sf)  

# Prepare data for plotting by adding a region column
data <- alldata %>%
  mutate(
    region = case_when(
      latitude >= 47.78184 & latitude <= 48.8208 ~ "North",
      latitude >= 47.38566 & latitude <= 47.7296079 ~ "Central",
      latitude >= 47.052362 & latitude <= 47.3539722 ~ "South"
    )
  ) %>%
  filter(!is.na(region) & !is.na(p450))

# Plotting p450 values by region over the Washington state map
ggplot(data = data, aes(x = longitude, y = latitude, color = p450)) +
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point() +
  scale_color_gradient(low = "blue", high = "red") +
  facet_wrap(~ region) +
  labs(title = "p450 Levels by Region in Washington State",
       x = "Longitude",
       y = "Latitude",
       color = "p450 Levels") +
  theme_minimal()

SOD map- WA state

smap<- ggplot() + 
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point(data = alldata, aes(x = longitude, y = latitude, color = SOD), size = 3) +
  scale_color_viridis(option = "C", name = "SOD") +
  theme_minimal() +
  labs(title = "Washington State - SOD Data")

print(smap)

ggsave(plot=smap, filename="/Users/cmantegna/Documents/WDFWmussels/output/sodmap.png", width=15, height=8)

#SOD map - Puget Sound

#zoom into puget sound region & note the legend, lighter colors are higher values
xlim <- c(-124, -122)  # longitude bounds
ylim <- c(47, 49)  # latitude bounds

sodpugetsound<- ggplot() + 
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point(data = alldata, aes(x = longitude, y = latitude, color = SOD), size = 3) +
  scale_color_viridis(option = "C", name = "SOD") +
  coord_sf(xlim = xlim, ylim = ylim, expand = FALSE)+
  theme_minimal() +
  labs(title = "Washington State - SOD Data")

print(sodpugetsound)

SOD map by regions

fix the coordinates and maps, not a priority 4/21

library(ggplot2)
library(dplyr)
library(sf)  

# Prepare data for plotting by adding a region column
data <- alldata %>%
  mutate(
    region = case_when(
      latitude >= 47.78184 & latitude <= 48.8208 ~ "North",
      latitude >= 47.38566 & latitude <= 47.7296079 ~ "Central",
      latitude >= 47.052362 & latitude <= 47.3539722 ~ "South"
    )
  ) %>%
  filter(!is.na(region) & !is.na(SOD))

# Plotting p450 values by region over the Washington state map
ggplot(data = data, aes(x = longitude, y = latitude, color = SOD)) +
  geom_sf(data = washington_map, fill = "lightgrey", color = "white") +
  geom_point() +
  scale_color_gradient(low = "blue", high = "red") +
  facet_wrap(~ region) +
  labs(title = "p450 Levels by Region in Washington State",
       x = "Longitude",
       y = "Latitude",
       color = "SOD Activity") +
  theme_minimal()