The purpose of this script is to evaluate the spatial and temporal overlap of the different surveys used in the SEM.

TASK LIST:

Inventory datasets:

Annual SEM: - Fish: FMWT, Bay Study, DJFMP Beach seine (ISS) - Zoop/nutrients: EMP - Benthic: EMP Benthic - Flows: Dayflow

Monthly SEM: - Fish: FMWT, Bay Study, DJFMP Beach seine (ISS) - Zoop/ nutrients: EMP - Benthic: EMP Benthic - Flows: Dayflow

Plot them:

Plot all monthly stations on one map
Evaluate overlap with region designations
Create new region map, based on station overlap
Export polygon layer and station designations
Use new polygon layer to get dayflow estimates

Load libraries

# tidy data
library(tidyverse)
library(ggthemes)
library(patchwork)

library(DT)
library(leaflet)
library(sf)
library(janitor)

#devtools::install_github("yutannihilation/ggsflabel")
library(ggsflabel)

#devtools::install_github("tomroh/leaflegend")
library(leaflegend)
require(knitr)

Fish surveys

Bay Study

bay_annual <- read_csv("data/stations/stations_fish_BayStudy_annual.csv") %>% 
  mutate(Category="Fish", Survey="Bay Study",
         Station=as.character(Station))

bay_month <- read.csv("data/stations/stations_fish_BayStudy_monthly.csv") %>% 
  mutate(Category="Fish", Survey="Bay Study",
         Station=as.character(Station)) %>% 
  select(Station, Latitude, Longitude, Survey, Category) # rearrange columns

Fall Midwater Trawl

FMWT_annual <- read_csv("data/stations/stations_fish_FMWT_annual.csv")%>% 
  mutate(Category="Fish", Survey="FMWT",
         Station=as.character(Station))

Delta Juvenile Fish Monitoring Program

seine_annual <- read_csv("data/stations/stations_fish_DJFMP_annual.csv")%>% 
  mutate(Category="Fish", Survey="DJFMP",
         Station=as.character(Station))

seine_month <- read_csv("data/stations/stations_fish_DJFMP_monthly.csv")%>% 
  mutate(Category="Fish", Survey="DJFMP",
         Station=as.character(Station))

Summer Townet

STN_annual <- read_csv("data/stations/stations_fish_STN_annual.csv")%>% 
  mutate(Category="Fish", Survey="STN",
         Station=as.character(Station))

Zoop surveys

zoop_annual <- read_csv("data/stations/stations_zoop_annual.csv")%>% 
  mutate(Category="Zoop", Survey="EMP Zoop")

zoop_month <- read_csv("data/stations/stations_zoop_month.csv")%>% 
  mutate(Category="Zoop", Survey="EMP Zoop")

Water quality

Same stations as EMP zooplankton.

Nutrient surveys

nutrients_annual <- read_csv("data/stations/stations_nutrients_annual.csv")%>% 
  mutate(Category="Nutrients", Survey="EMP Nutrients")

nutrients_month <- read_csv("data/stations/stations_nutrients_monthly.csv")%>% 
  mutate(Category="Nutrients", Survey="EMP Nutrients")

Benthic surveys

benthic_month <- read_csv("data/stations/stations_benthic_monthly.csv")%>% 
  mutate(Category="Clams", Survey="EMP Benthic") %>% 
  select(Station:Survey)

benthic_annual <- read_csv("data/stations/stations_benthic_annual.csv")%>% 
  mutate(Category="Clams", Survey="EMP Benthic") %>% 
  select(Station:Survey)

Pesticide data

TBD.

# library(dataRetrieval)
# ?dataRetrieval

Stitch together

Monthly

station_list_month <- bind_rows(bay_month, seine_month, zoop_month, benthic_month, nutrients_month)%>%
  mutate(Survey=factor(Survey),
         Station=factor(Station))

station_list_month %>% 
  write_csv("data/stations/station_list_for_monthly_analysis.csv")

Annual

station_list_annual <- bind_rows(bay_annual, FMWT_annual, STN_annual, 
                                 seine_annual, zoop_annual, benthic_annual,
                                 nutrients_annual)%>%
  mutate(Survey=factor(Survey),
         Station=factor(Station))

station_list_annual %>% 
  write_csv("data/stations/station_list_for_annual_analysis.csv")

Regional designations

Code from deltamapr documentation. https://github.com/InteragencyEcologicalProgram/deltamapr#regions

library(deltamapr)
library(ggplot2)
library(sf)

# query region
deltamapr::R_EDSM_Regions_1617P1

## Simple feature collection with 4 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 558333.4 ymin: 4182233 xmax: 650325.5 ymax: 4273755
## Projected CRS: NAD83 / UTM zone 10N
##     Region                       geometry
## 1 Far West POLYGON ((562664.4 4227189,...
## 2    North POLYGON ((627665.6 4226993,...
## 3    South POLYGON ((640148.7 4220717,...
## 4     West POLYGON ((615648.7 4210824,...

st_crs(R_EDSM_Regions_1617P1) <- 26910
ggplot(R_EDSM_Regions_1617P1)+
   geom_sf(aes(fill=Region))+
   theme_bw()

deltamapr::R_EDSM_Subregions_Mahardja

## Simple feature collection with 41 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 534372.4 ymin: 4140336 xmax: 650325.5 ymax: 4273755
## Projected CRS: NAD83 / UTM zone 10N
## First 10 features:
##      Region                             SubRegion       SQM
## 1     North                     Upper Yolo Bypass 471689278
## 2  Far West                      Upper Napa River  21496636
## 3     North            Sacramento River near Ryde  92220125
## 4     South                 Upper Mokelumne River 213170012
## 5     North       Sacramento River near Rio Vista  44594324
## 6     South San Joaquin River at Twitchell Island  41103028
## 7     South     San Joaquin River at Prisoners Pt  31760557
## 8     South                 Disappointment Slough 182085789
## 9     South                          Franks Tract  94761329
## 10    South                           Holland Cut  28298743
##                          geometry
## 1  POLYGON ((618378.8 4244972,...
## 2  POLYGON ((560777 4227437, 5...
## 3  POLYGON ((630208.1 4233426,...
## 4  POLYGON ((639477.9 4220463,...
## 5  POLYGON ((617164.8 4226544,...
## 6  POLYGON ((622278.3 4220219,...
## 7  POLYGON ((632383.9 4213297,...
## 8  POLYGON ((650317 4202723, 6...
## 9  POLYGON ((624582.5 4215447,...
## 10 POLYGON ((627223.3 4211334,...

st_crs(R_EDSM_Subregions_Mahardja) <- 26910

ggplot(R_EDSM_Subregions_Mahardja)+
  geom_sf(aes(fill=Region))+
  geom_sf_label(aes(label=SubRegion), size=2) +
  theme_bw()+
  ggthemes::scale_fill_colorblind()

Monthly map: leaflet

Extract subregions of interest and combine

Regions <- R_EDSM_Subregions_Mahardja %>%
  filter(!SubRegion %in% c("South Bay", "San Francisco Bay", "Upper Napa River", "Upper Yolo Bypass",
                           "Cache Slough and Lindsey Slough", "Liberty Island", "Upper Sacramento River Ship Channel",
                           "Grant Line Canal and Old River", "Disappointment Slough", "Rock Slough and Discovery Bay",
                           "Steamboat and Miner Slough", "Franks Tract"))

# plot with single dataframe for simplicity

pal2 <- colorFactor(
  palette = c("green", "purple", "pink", "black", "white", "yellow"),
  domain = station_list_month$Survey
)

leaflet() %>% 
  addTiles() %>% 
  addPolygons(data=Regions %>% 
              st_transform(crs=4326),
              color = ~pal(SubRegion),
              label = ~SubRegion,
              group = "polygons",
              weight=2,
              opacity=1) %>% 
    addCircleMarkers(data = station_list_month,  
                         lat = ~Latitude,
                         lng = ~Longitude,
                          color = ~pal2(Survey),
                          radius = 2,
                          weight = 2,
                          opacity= 1,
                          label = ~Station) %>% 
  # help here-- need the legend to show one label for each survey type
    addLegend(data=station_list_month,
                          "bottomleft",
                          #values = ~Survey,
                          colors = ~pal2(Survey),
                          title = "Survey stations",
                          labels = ~Survey,
                          opacity = 1,
                          group = "circles", 
                          labFormat = labelFormat(style = list(
                  "font-family" = "serif",
                  "font-size" = "8px"))) %>%
    addLayersControl(baseGroups = c("polygons"))

stations as simple features

station_month_tbl <- tibble(lat=station_list_month$Latitude, lon=station_list_month$Longitude)

station_month_points <- st_as_sf(station_list_month, coords = c("Longitude", "Latitude"), 
                 crs= 4326, agr = "constant") %>% 
                  st_transform(crs=4326)

station_annual_tbl <- tibble(lat=station_list_annual$Latitude, lon=station_list_annual$Longitude)

station_annual_points <- st_as_sf(station_list_annual, coords = c("Longitude", "Latitude"), 
                 crs= 4326, agr = "constant") %>% 
                  st_transform(crs=4326)

regions as simple features

URL: https://r-spatial.github.io/sf/reference/geos_combine.html

Pilot code for joining regions with st_union function Compare different regional designations for San Pablo Bay as an example

# san pablo bay only
sanpablo1 <- Regions %>% 
  filter(SubRegion == "San Pablo Bay") %>% 
  mutate(Designation = "San Pablo Bay")
  
sanpablo1_union <- st_union(sanpablo1$geometry, by_feature= FALSE)
plot(sanpablo1_union)

# san pablo bay, lower napa, and carquinez strait
sanpablo2 <- Regions %>% 
  filter(SubRegion == "San Pablo Bay" |
        SubRegion == "Lower Napa River" |
        SubRegion == "Carquinez Strait") %>% 
  mutate(Designation = "San Pablo Bay")

sanpablo2_union <- st_union(sanpablo2$geometry, by_feature= FALSE)
plot(sanpablo2_union)

scenario 1, 2, 3, 4

region_scenarios <- Regions %>% 
      mutate(Scenario1 = case_when(
                                  # far west
                                  SubRegion == "San Pablo Bay" | 
                                  SubRegion == "Carquinez Strait" | 
                                  SubRegion == "Lower Napa River" ~ "Far West",
                                  # west
                                  SubRegion == "West Suisun Bay" |
                                  SubRegion == "Mid Suisun Bay" |
                                  SubRegion == "Grizzly Bay" | 
                                  SubRegion == "Suisun Marsh" |
                                  SubRegion == "Honker Bay" |
                                  SubRegion == "Confluence" ~ "West",
                                  # north
                                  SubRegion == "Lower Sacramento River" ~ "North",
                                  Region == "North" ~ "North",
                                  # south
                                  SubRegion == "Lower San Joaquin River" ~ "South",
                                  Region == "South" ~ "South",
                                  TRUE ~ "Other")) %>% 
    mutate(Scenario2 = case_when(
                                  # far west
                                  SubRegion == "San Pablo Bay" |
                                  SubRegion == "Lower Napa River" ~ "Far West",
                                  # north
                                  Region == "North" ~ "North",
                                  # south
                                  Region == "South" ~ "South",
                                  # west
                                  SubRegion == "Lower Sacramento River" ~ "West",
                                  TRUE ~ "West")) %>% 
  mutate(Scenario3 = case_when (
                                  # far west
                                  SubRegion == "San Pablo Bay" |
                                  SubRegion == "Lower Napa River" ~ "Far West",
                                  # north
                                  Scenario1 == "North" ~ "North",
                                  # south
                                  Scenario1 == "South" ~ "South",
                                  # west
                                  TRUE ~ "West")) %>% 
  mutate(Scenario4 = case_when(
                                  # far west
                                  Scenario1 == "Far West" ~ "Far West",
                                  # north
                                  Scenario2 == "North" ~ "North",
                                  # south
                                  Scenario2 == "South" ~ "South",
                                  # west
                                  TRUE ~ "West"))%>%
  st_transform(crs=4326)


# use st_union code to combine these designated polygons for each scenario

# scenario 1 - two approaches

plot1 <- ggplot(region_scenarios)+
  geom_sf(aes(fill=Scenario1))+
  scale_fill_viridis_d(name="Region")+
  theme_bw()+
  labs(title="Scenario 1")


# scenario 2

plot2 <- ggplot(region_scenarios)+
  geom_sf(aes(fill=Scenario2))+ 
  labs(title="Scenario 2")+
  scale_fill_viridis_d(name="Region")+
  theme_bw()+
  labs(title="Scenario 2")


# scenario 3
plot3 <- ggplot(region_scenarios)+
  geom_sf(aes(fill=Scenario3))+ 
  labs(title="Scenario 3")+
  scale_fill_viridis_d(name="Region")+
  theme_bw()+
  labs(title="Scenario 3")


# scenario 4
plot4 <- ggplot(region_scenarios)+
  geom_sf(aes(fill=Scenario4))+ 
  labs(title="Scenario 4")+
  scale_fill_viridis_d(name="Region")+
  theme_bw()+
  labs(title="Scenario 4")


plot1 + plot2 + plot3 + plot4 + plot_layout(nrow=2, guides="collect")

station densities

Monthly

station_month_points_filtered<-station_month_points%>%
  st_filter(region_scenarios)

# scenario 1 with stations
plot1b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario1))+
  geom_sf(data=station_month_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(1,2,15,18,4))+
  labs(title="Scenario 1: Monthly", subtitle="1990s to present")

# scenario 2 with stations
plot2b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario2))+
  geom_sf(data=station_month_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(1,2,15,18,4))+
  labs(title="Scenario 2: Monthly", subtitle="1990s to present")

# scenario 3 with stations
plot3b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario3))+
  geom_sf(data=station_month_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(1,2,15,18,4))+
  labs(title="Scenario 3: Monthly", subtitle="1990s to present")

# scenario 4 with stations
plot4b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario4))+
  geom_sf(data=station_month_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(1,2,15,18,4))+
  labs(title="Scenario 4: Monthly", subtitle="1990s to present")

plot1b + plot2b + plot3b + plot4b + plot_layout(nrow=2, guides="collect")

Ugly code to identify which regions they fall in

scenario_month_counts<-station_month_points_filtered%>%
  st_join(region_scenarios%>%
            select(-Region, -SQM)%>%
            pivot_longer(cols=contains("Scenario"), names_to = "Scenario", values_to="Region", names_prefix="Scenario")
  )%>%
  st_drop_geometry()%>%
  group_by(Scenario, Region, Survey)%>%
  summarise(Count=n(), .groups="drop")%>%  
  pivot_wider(names_from="Survey", values_from="Count")%>% 
  mutate(Analysis="Monthly") %>% 
  select(Analysis, Scenario, Region, `Bay Study`: `EMP Zoop`)

kable(scenario_month_counts)

Analysis	Scenario	Region	Bay Study	EMP Benthic	EMP Nutrients	EMP Zoop
Monthly	1	Far West	12	2	1	1
Monthly	1	North	6	2	1	2
Monthly	1	South	5	2	2	4
Monthly	1	West	8	2	3	6
Monthly	2	Far West	10	2	1	1
Monthly	2	North	4	1	NA	1
Monthly	2	South	3	2	2	3
Monthly	2	West	14	3	4	8
Monthly	3	Far West	10	2	1	1
Monthly	3	North	6	2	1	2
Monthly	3	South	5	2	2	4
Monthly	3	West	10	2	3	6
Monthly	4	Far West	12	2	1	1
Monthly	4	North	4	1	NA	1
Monthly	4	South	3	2	2	3
Monthly	4	West	12	3	4	8

scenario_month_counts %>% 
  write_csv("data/stations/region_scenarios_and_station_counts_monthly.csv")

Annual

station_annual_points_filtered<-station_annual_points%>%
  st_filter(region_scenarios)

# scenario 1 with stations
plot1b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario1))+
  geom_sf(data=station_annual_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(20,1,2,15,18,4,0))+
  labs(title="Scenario 1: Annual", subtitle="1980s to present")

# scenario 2 with stations
plot2b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario2))+
  geom_sf(data=station_annual_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(20,1,2,15,18,4,0))+
  labs(title="Scenario 2: Annual", subtitle="1980s to present")

# scenario 3 with stations
plot3b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario3))+
  geom_sf(data=station_annual_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(20,1,2,15,18,4,0))+
  labs(title="Scenario 3: Annual", subtitle="1980s to present")

# scenario 4 with stations
plot4b <- ggplot()+
  geom_sf(data=region_scenarios, aes(fill=Scenario4))+
  geom_sf(data=station_annual_points_filtered, aes(color=Survey, shape=Survey))+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c(20,1,2,15,18,4,0))+
  labs(title="Scenario 4: Annual", subtitle="1980s to present")

plot1b + plot2b + plot3b + plot4b + plot_layout(nrow=2, guides="collect")

Ugly code to identify which regions they fall in

scenario_annual_counts<-station_annual_points_filtered%>%
  st_join(region_scenarios%>%
            select(-Region, -SQM)%>%
            pivot_longer(cols=contains("Scenario"), names_to = "Scenario", values_to="Region", names_prefix="Scenario")
  )%>%
  st_drop_geometry()%>%
  group_by(Scenario, Region, Survey)%>%
  summarise(Count=n(), .groups="drop")%>%  
  pivot_wider(names_from="Survey", values_from="Count")%>% 
  mutate(Analysis="Annual") %>% 
  select(Analysis, Scenario, Region, `Bay Study`: `EMP Zoop`)

kable(scenario_annual_counts)

Analysis	Scenario	Region	Bay Study	FMWT	STN	DJFMP	EMP Benthic	EMP Nutrients	EMP Zoop
Annual	1	Far West	1	8	1	NA	NA	NA	NA
Annual	1	North	1	10	4	5	1	1	2
Annual	1	South	1	17	7	2	1	2	3
Annual	1	West	7	36	13	NA	1	3	5
Annual	2	North	NA	6	2	4	NA	NA	1
Annual	2	South	NA	13	6	NA	1	2	2
Annual	2	West	10	52	17	3	2	4	7
Annual	3	North	1	10	4	5	1	1	2
Annual	3	South	1	17	7	2	1	2	3
Annual	3	West	8	44	14	NA	1	3	5
Annual	4	Far West	1	8	1	NA	NA	NA	NA
Annual	4	North	NA	6	2	4	NA	NA	1
Annual	4	South	NA	13	6	NA	1	2	2
Annual	4	West	9	44	16	3	2	4	7

scenario_annual_counts %>% write_csv("data/stations/region_scenarios_and_station_counts_annual.csv")

stitch together

scenario_counts_all <- bind_rows(scenario_month_counts, scenario_annual_counts)
scenario_counts_all %>% write_csv("data/stations/region_scenarios_and_station_counts_all.csv")

Filtering the overall region to the area of interest, then saving it for use in data processing scripts

region_scenarios_cropped<-region_scenarios%>%
  filter(!SubRegion%in%c("Lower Sacramento River Ship Channel", "Upper Sacramento River", "Middle Sacramento River", 
                         "Sacramento River near Ryde", "Georgiana Slough", "Upper Mokelumne River", "Lower Mokelumne River", 
                         "San Joaquin River near Stockton", "Upper San Joaquin River", "Mildred Island", 
                         "Middle River", "Victoria Canal"))%>%
  group_by(Scenario3)%>%
  summarise()%>% 
  st_make_valid()%>%
  rename(Region=Scenario3)

Plot the spatial extent

ggplot()+
  geom_sf(data=deltamapr::WW_Delta, color="slategray1", fill="slategray1")+
  geom_sf(data=region_scenarios_cropped%>%st_transform(crs=st_crs(deltamapr::WW_Delta)), aes(fill=Region), alpha=0.4)+
  theme_bw()

Plot the stations within the annual or monthly regional extents

plot1<-ggplot()+
  geom_sf(data=region_scenarios_cropped, aes(fill=Region))+
  geom_sf(data=station_month_points_filtered%>%mutate(Survey=factor(Survey, levels=levels(station_annual_points_filtered$Survey))), aes(color=Survey, shape=Survey))+
  coord_sf(xlim=st_bbox(station_month_points_filtered)[c(1,3)], ylim=st_bbox(station_month_points_filtered)[c(2,4)])+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey", drop=FALSE)+
  scale_shape_manual(name="Survey", values=c("Bay Study"=1, "DJFMP"=2, "EMP Benthic"=15, 
                                             "EMP Nutrients"=18, "EMP Zoop"=4, "FMWT"=20),
                     drop=FALSE)+
  labs(title="Scenario 3: Monthly cropped", subtitle="1990s to present")

plot2<- ggplot()+
  geom_sf(data=region_scenarios_cropped, aes(fill=Region))+
  geom_sf(data=station_annual_points_filtered, aes(color=Survey, shape=Survey))+
  coord_sf(xlim=st_bbox(station_month_points_filtered)[c(1,3)], ylim=st_bbox(station_month_points_filtered)[c(2,4)])+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c("Bay Study"=1, "DJFMP"=2, "EMP Benthic"=15, 
                                             "EMP Nutrients"=18, "EMP Zoop"=4, "FMWT"=20))+
  labs(title="Scenario 3: Annual cropped", subtitle="1980s to present")

plot1 + plot2 + plot_layout(nrow=2, guides="collect")

Now how many stations are in each region with the cropped extents?

scenario_counts_cropped<-station_annual_points_filtered%>% 
  mutate(Analysis="Annual") %>%
  bind_rows(station_month_points_filtered%>%
              mutate(Analysis="Monthly"))%>%
  st_join(region_scenarios_cropped)%>%
  filter(!is.na(Region))%>%
  st_drop_geometry()%>%
  group_by(Analysis, Region, Survey)%>%
  summarise(Count=n(), .groups="drop")%>%  
  pivot_wider(names_from="Survey", values_from="Count")%>% 
  select(Analysis, Region, `Bay Study`: `EMP Zoop`)

kable(scenario_counts_cropped)

Analysis	Region	Bay Study	DJFMP	EMP Benthic	EMP Nutrients	EMP Zoop
Annual	North	1	5	1	1	2
Annual	South	1	2	1	2	3
Annual	West	8	NA	1	3	5
Monthly	Far West	10	NA	2	1	1
Monthly	North	6	5	2	1	2
Monthly	South	5	4	2	2	4
Monthly	West	10	NA	2	3	6

Output the regions

# For use in R
saveRDS(object=region_scenarios_cropped, file="data/data_in/regions.Rds")

# Shapefile for use anywhere
write_sf(region_scenarios_cropped%>%st_transform(crs=26910), "data/data_in/regions_shapefile/regions.shp")

#HELP PLZ - code is breaking at line 657 # Map figure for ppt

# insert ggplot code
p<-ggplot()+
  geom_sf(data=region_scenarios_cropped, aes(fill=Region))+
  geom_sf(data=station_month_points_filtered%>%mutate(Survey=factor(Survey, levels=levels(station_annual_points_filtered$Survey))), aes(color=Survey, shape=Survey))+
  coord_sf(xlim=st_bbox(station_month_points_filtered)[c(1,3)], ylim=st_bbox(station_month_points_filtered)[c(2,4)])+
  scale_fill_grey(name="Region", start = 0.4)+
  scale_color_viridis_d(name="Survey")+
  scale_shape_manual(name="Survey", values=c("Bay Study"=1, "DJFMP"=2, "EMP Benthic"=15, 
                                             "EMP Nutrients"=18, "EMP Zoop"=4))+
  labs(title="Monthly dataset", subtitle="1990s to present")

ggsave(p, file="data/data_in/regions.png",
       device="png", width=6, height=5, units="in")

Survey extents

Denise-Colombano

11/5/2021

Load libraries

Fish surveys

Bay Study

Fall Midwater Trawl

Delta Juvenile Fish Monitoring Program

Summer Townet

Zoop surveys

Water quality

Nutrient surveys

Benthic surveys

Pesticide data

Stitch together

Regional designations

Monthly map: leaflet

Extract subregions of interest and combine

stations as simple features

regions as simple features

scenario 1, 2, 3, 4

station densities

Monthly

Annual

stitch together

Filtering the overall region to the area of interest, then saving it for use in data processing scripts