scripts to explore the gslea EA.data and create figures and maps

EA_GSL.R

devtools::install_github("duplisea/gslea", build_vignettes = TRUE)
packages = c('tidyverse','magrittr', 'lubridate', 'data.table','gslea','formattable', 'sf', 'sp', 'tmap', 'rgdal','rgeos' ,'raster','maptools','mapdata','gridExtra','shapefiles','installr','crosstalk','flexdashboard','DT','rmapshaper' )
invisible(lapply(packages, function(x) {if (!require(x, character.only = T)) {install.packages(x);require(x)}}))

# Load data
load("./EA.data.rda")
load("./variable.description.rda")
load("./field.description.rda")
AEbox_boundary <- read_csv("./AEbox.boundary.csv")

# Load map ... I made this shape file from the AE bounding boxes and can provide it
EAR_map <- st_read("./EAR_map.shp")

# Add variable descriptions to the data set
EA.data <- left_join(EA.data, variable.description, by = "variable")
glimpse(EA.data)

# split climatic and non climatic data
EA.clim <- EA.data %>% dplyr::filter(type == "climatic")
EA.data <- EA.data %>% dplyr::filter(type != "climatic")

# Change variable classes
EA.data <- as_tibble(EA.data)
EA.data$value <- as.numeric(EA.data$value)
glimpse(EA.data)

# add common grouping variable
EAR_map %<>% mutate(EAR = as.character(PID))

# change variable class
EA.data$EAR <- as.character(EA.data$EAR)

#simplify geometry
EAR_map_lowfi <- ms_simplify(EAR_map)

# merge data with spatial data
EAR <- left_join(EA.data, EAR_map_lowfi)
#EAR <- left_join(phys, EAR_map_lowfi)

# Rename EAR variable 
level.key <- c("1" = "northwest_estuary", 
               "2" = "northeast",
               "3" = "centre",
               "4" = "mecatina",
               "5" = "magdalen_shallows",
               "6" = "northumberland_strait",
               "7" = "laurentian_hermitage",
               "50" = "baie_des_chaleurs",
               "10" = "estuary",
               "-1" = "climatic") 

# recode
EAR$EAR <- recode(EAR$EAR, !!!level.key)
EAR_std$EAR <- recode(EAR_std$EAR, !!!level.key)
unique(EAR_std$EAR)


# transform to simple features (sf) and set projection to wgs84
EAR_sf <- st_as_sf(EAR)
st_crs(EAR_sf) <- 4326

# Make standardised data df
EA.data.wide <- as_tibble(EA.data %>% 
                            dplyr::filter(!is.na(value)) %>%
                            dplyr::select(year, EAR, variable, value) %>%
                            pivot_wider(names_from = "variable", values_from = "value"))

# make standardized values (i.e. zscore, aka anomalies) - caution here as some variables are already standardised
EAR_std <- EA.data.wide %>% 
pivot_longer(cols = 3:112, names_to = "variable", values_to = "value")  %>% # with addition of new variables the columns need to be revised
  dplyr::filter(!is.na(value)) %>% 
  group_by(year, EAR, variable) %>% 
  dplyr::summarise(value = mean(value,na.rm=T))
EAR_std <- as.data.frame(EAR_std)
EAR_std %<>% 
  mutate(year = as.character(year), EAR = as.character(EAR)) %>% 
  group_by(EAR, variable) %>% 
  mutate(m_val = mean(value,na.rm=T), sd_val = sd(value, na.rm = T), value_std = (value-m_val)/sd_val)
EAR_std <- left_join(EAR_std, variable.description)
EAR_std <- as.data.frame(EAR_std)

# some example figures using tmap

## SST in july maps by year
tm_shape(EAR_sf %>% 
           dplyr::filter(variable == "SST.month7")) +
  tm_polygons("value", palette = "plasma", style = "fixed", breaks = c(0,3.5,7,9.5,13,16.5,20)) +
  tm_layout(frame = FALSE, bg.color = "transparent") + 
  tm_layout(legend.text.size = 1,
            legend.show =  TRUE,
            legend.outside = TRUE,
            legend.outside.position = "right",
            panel.label.size = 1,
            panel.label.bg.color = "white",
            panel.show = T,
            scale = 1, 
            outer.margins = c(0,0,0,0),
            inner.margins = c(0,0,0,0),
            between.margin = 0) +
  tm_facets(by = "year", free.coords = FALSE, ncol = 4)
  
  # plots of the taylor gulf stream north wall over time
  EA.clim %>%
           dplyr::filter(variable %in% 
                           c("T.GSNW.Q1","T.GSNW.Q2","T.GSNW.Q3","T.GSNW.Q4"),
                         year > 1899, year < 2020) %>% 
           ggplot(aes(year, value)) +
           geom_point() +
           geom_smooth(span = 0.05) +
           scale_x_continuous(breaks=seq(1950, 2020, by = 5)) +
           geom_hline(yintercept = 0)