DFO-NOAA-Pacific / gfsynopsis-noaa

Code for extending the gfsynopsis report idea from DFO to Bering Sea / Gulf of Alaska / West Coast of USA
https://dfo-noaa-pacific.github.io/gfsynopsis-noaa/
GNU Affero General Public License v3.0
1 stars 0 forks source link

Data Visualization #10

Open chantelwetzel-noaa opened 2 days ago

chantelwetzel-noaa commented 2 days ago

Plotting or table functions need to be developed to visualize the following data products:

library(ggplot2)
library(tidyverse)

# bring in sdmTMB results of the the predicted index 
# From: predict.sdmTMB {sdmTMB}
df = readRDS(rds output from an sdmTMB model)
# get just predictions for the prediction grid
p = data.frame(df[[3]]) 
# get CPUE from delta-gamma
p$cpue = plogis(p$est1) * exp(p$est2) 
# scale to max of 1.0 just for plotting
# could also just leave as cpue
p = p %>% mutate(`Biomass index` = cpue/max(p$cpue) )

# start mapping, set up base map
states <- map_data("state")
west_coast <- subset(states, region %in% c("california", "oregon", "washington"))

# some fidgety stuff to plot years along the x-axis. 
axis.mod = 3
max.year = 2022
p$yrs = axis.mod*(p$Year-max.year)
p$lon2 = p$lon + p$yrs
xmin = floor( min(p$lon2))

ggplot(data=west_coast, aes(x = long, y = lat), 
                      fill = grey(0.9), color = "black") +
  geom_polygon(aes(x = long, y = lat, group=group), fill = grey(0.9), color = "black") +
  geom_point(data=p, aes(lon2, lat, color=`Biomass index`), size=0.01) +
  # might need to adjust here based on the scale of the data, esp 'trans' term.
  scale_color_viridis_c(option = 'turbo', direction = 1, trans='log', breaks=c(0,0.001,0.01,0.1,1)) +
  coord_fixed(xlim=c(xmin,-117),ylim=c(32,48),ratio=1.3) +
  # might need to adjust x-axes to plot all data
  scale_x_continuous(breaks = seq(-177,-120,axis.mod ), minor_breaks = -117:-177, labels = 2003:max.year) +
  xlab("Year") +
  ylab("Latitude") +
  theme_bw() + 
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.title = element_text(size = 8),
        legend.text = element_text(size = 8),
        legend.key.size = unit(1,'lines'))