R Notebook

Set up

First make sure to install and load all packages needed. This chunk will also load functions related to the analysis.

# Load libraries and functions #

ipak <- function(pkg){
  new.pkg <- pkg[!(pkg %in% installed.packages()[, "Package"])]
  if (length(new.pkg)) 
    install.packages(new.pkg, dependencies = TRUE)
  sapply(pkg, require, character.only = TRUE)
}


#### Library ####
packages <- c(
  "readxl", # Read dataframe in excel format
  "dplyr", # Data manipulation
  "tidyr", # Data manipulation
  "ggplot2", #Nice grpahs and spatial analysis
  # "cowplot",
  # "forcasts",
  # "rgdal",
  # "RColorBrewer",
  "knitr",
  "kableExtra",
  "data.table", # Fast read data
  # "ggrepel",
  # "gridExtra",
  # "ggmap",
  # "rgeos",
  "sf",
  # "sp",
  # "rgdal", #Spatial analysis
  "tools", #Spatial analysis
  # "png", # For reading plots in chunk codes
  # "grid" # For reading plots in chunk codes
  "stringr" # For filtering by string
)

ipak(packages)
##__________________FUNCTIONS _________________________##

# Standarization of plots
ggtheme_plot <- function() {
  theme(
    plot.title = element_text(size = rel(1), hjust = 0, face = "bold"),
    panel.background = element_blank(),
    strip.background = element_blank(),
    # strip.text       = element_text(size = base_size, face = "italic"),
    panel.border     = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_blank(),
    axis.ticks = element_blank(),
    axis.text.x = element_text(size = 10,
                               angle = 0,
                               face = "plain"),
    axis.text.y = element_text(size = 12),
    axis.title = element_text(size = 14),
    legend.key = element_rect(colour = NA, fill = NA),
    legend.position  = "top",
    legend.title = element_text(size = 16),
    legend.text = element_text(size = 14),
    strip.text.x = element_text(size = 18, colour = "darkgrey"),
    strip.text = element_text(size = 18)
  )
}

# Standarization of maps
ggtheme_map <- function(base_size = 9, Region = "NA") {
  
  theme(text             = element_text(#family = "Helvetica",
    color = "gray30", size = base_size),
    plot.title       = element_text(size = rel(1.25), hjust = 0, face = "bold"),
    panel.background = element_blank(),
    panel.border     = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major = element_line(color = "transparent"),
    strip.background =element_rect(fill = "transparent"),
    strip.text.x = element_text(size = 18, colour = "black",face= "bold.italic", angle = 0),
    axis.line        = element_blank(), # element_line(colour = "grey30", size = .5))
    axis.ticks       = element_blank(),
    axis.text        = element_blank(),
    axis.title       = element_blank(),
    # legend.key       = element_rect(colour = NA, fill = NA, size = 4),
    legend.position = "bottom"#,
    # legend.key.width =unit(6,"line")
  )
}

Overall Map

Loading DBEM data

Data Wrangling

• We need to make some corrections to the shapefile to project everything with North America centered. Basically solve for the fact that Alaska is in the other side of the globe because the Earth is round, shocking!

• We also have to include a column referencing the country’s EEZ for the analysis

# LAT LONG Corrections for map projection

DBEM_Corrected <- DBEM_Results %>% 
  filter(INDEX %in% All_EEZs$INDEX) %>% 
  filter(Longitude > 0) %>% 
  mutate(Longitude = Longitude - 360)

DBEM_Corrected <- DBEM_Results %>% 
  filter(INDEX %in% All_EEZs$INDEX) %>% 
  filter(Longitude < 0) %>% 
  bind_rows(DBEM_Corrected)

# head(DBEM_Corrected)

# Just checking #####

# DBEM_Corrected %>%
#   group_by(INDEX, Longitude, Latitude) %>%
#   summarise(n =n()) %>%
#   ggplot(.,
#          aes(x = Longitude,
#              y = Latitude,
#              fill = n
#          ),
#          size = 0.05,
#          alpha = 0.5) +
#   geom_tile() #+
# coord_map(projection = "mercator")


### Get ridd of any Pacific Species in NAFO

NAFO_E <- DBEM_Corrected %>% 
  filter(Longitude >= -100)

Pacific <- DBEM_Corrected %>% 
  filter(Treaty != "NAFO")

DBEM_Corrected <- NAFO_E %>% 
  bind_rows(Pacific)

# unique(DBEM_Corrected$Species)

# Remove salmons...

Salmons <- c("Oncorhynchus keta","Oncorhynchus kisutch","Oncorhynchus nerka","Oncorhynchus tshawytscha","Oncorhynchus gorbuscha","Salmo salar")

Chapter_Data <- DBEM_Corrected %>%
  # filter(!Species %in% Salmons) %>% # Reomve or only analysis with salmons
  left_join(All_EEZs,
            by ="INDEX") %>% 
  select(-EEZID)

# Show table with species 

# names(Chapter_Data)
# head(Chapter_Data)

Chapter_Data %>% 
  group_by(Species, Treaty) %>% 
  summarise(n=n()) %>% 
  ungroup() %>% 
  left_join(exploited_species_list,
            by = "Species") %>% 
  arrange(Treaty) %>% 
  select(Species,CommonName,Treaty) %>% 
  kable()

Loading Spatial data

# Countries included 

Countries <- c("Canada","United States","Alaska")


# World shapefile
path_world <- "~/Documents/R Code/Spatial_Analysis/Data/TM_WORLD_BORDERS-0.3"
file_name <- "TM_WORLD_BORDERS"

World_Land_sf <- st_read(dsn = path_world, 
                         layer = file_path_sans_ext(file_name)
) %>% 
  st_transform(crs = 3832)

North_America_Land <- World_Land_sf %>% 
  filter(NAME %in% Countries)

# Checkpoint 
# unique(North_America_Land$NAME) [1] Canada        Greenland     United States

# Map it? Yup!
# ggplot(North_America_Land) +
# geom_sf()

# EEZs shapefile ####

path_eez<- ("~/Documents/R Code/Spatial_Analysis/Data/World_EEZ_v10_2018")
#The File
fnam_eez <- "eez_boundaries_v10.shp" # Just lines
# fnam_eez <- "eez_v10.shp" # Actual shapefiles (takes 4-ever, for ever ever? fooor eever eeever!?)

eez_world_sf <- st_read(dsn = path_eez, 
layer = file_path_sans_ext(fnam_eez)
) %>% 
  st_transform(crs = 4326) #3832 centered

# head(eez_world_sf)

#Checkpoint  
# unique(eez_world_sf$Territory1) 

# North_America_EEZ <- eez_world_sf %>% 
#   filter(Territory1 %in% Countries) %>% 
#   st_transform(crs = 4326)

#Checkpoint
# unique(North_America_EEZ$Territory1)
# Map it? 
# ggplot(North_America_EEZ) +
#   geom_sf()
# Nope, Alaska is not showing up. Neither East Greenland

# Greenland Shapefile Bug ####
# First lets get all of greenalnd in the map... 
# 
# eez_world_sf %>%
#   filter(str_detect(Territory1,"Canada")) %>%
#   # mutate(Basin = ifelse( x_1 >= -100, "Atlantic","Pacific")) %>%
#   head()
#   ggplot() +
#   # geom_sf(aes(fill = Basin))
#   geom_sf()
# 
# 
# 
# # The chunk missing is from Iceland
# Iceland <-  eez_world_sf %>% 
#   # group_by(Territory1) %>% 
#   # summarise(n()) %>% 
#   # filter(str_detect(Territory1,"Iceland")) # %>% 
#   filter(Line_ID ==135)
#   # ggplot() +
#   # geom_sf()
# 
# eez_world_sf %>%
#   filter(str_detect(Territory1,"Greenland") |
#          Line_ID == 135) %>% # Include this in final map plot
#   ggplot() +
#   geom_sf() #+
#   # geom_sf(data = Iceland)
# 
# # Alaska Shapefile Bug ####
# eez_world_sf %>% 
#   filter(Territory1 %in% Countries) %>% 
#   st_transform(crs = 3832) %>% 
#   ggplot() +
#   geom_sf()
#   
# # The chunk missing is from Russia
# eez_world_sf %>%
#   filter(str_detect(Territory1,"Alaska") |
#            MRGID_Ter2 == 8684) %>% # This includes Russia side (Include in final plot)
#   st_transform(crs = 3832) %>%
#   # View()
#   ggplot() +
#   geom_sf(aes(col=Line_type))
# 
# 
# # South US (Atlantic) Shapefile Bug ####
# 
# # The chunk missing is from ... Mexico and Cuba?... and Bahamas. Nope, just Mexico
# eez_world_sf %>%
#   # filter(str_detect(Territory1"United States")) %>% 
#   filter(Territory1 %in% c ("United States", "Mexico", "Cuba","Bahamas"))
# 
# 
# ### Mexico First (Taht's it! )
# 
# eez_world_sf %>%
#   # filter(str_detect(Territory1"United States")) %>% 
#   filter(Territory1 %in% c ("United States") |
#          MRGID_Ter2 == 2204) %>%  # Inlcude this for Mexico's southern boundary
#   # st_transform(crs = 3832) %>%
#   # View()
#   ggplot() +
#   geom_sf(aes(col=Line_type))


#### The Actual Map ####

# We need to inlcude alaska in here
Countries <- c("Canada","United States", "Alaska")

# Here we go again
North_America_EEZ <- eez_world_sf %>% 
  filter(Territory1 %in% Countries |
           Line_ID == 135 | # # Line missing from Greenland (Iceland)
         MRGID_Ter2 == 8684 |
           MRGID_Ter2 == 2204) %>%  # Line missing from Alaska (Russia)
  st_transform(crs = 3832)

# Finall check! 
# ggplot(North_America_EEZ) +
#   geom_sf()

# Yup! 

# Map both shapefiles to see how it looks

# ggplot() +
#   geom_sf(data = North_America_Land) +
#   geom_sf(data = North_America_EEZ)

# Annnnd.... Is... Done!!!!!!!

# But not really because those are lines and we actually need polygons...
#  Its way easier, they are all there... ¬ ¬ It does take for ever, thou so projection issues should be deal with in the line version. But maybe the line fizes cuz they are polygons? We will see...

North_America_EEZ <- eez_world_sf %>% 
  filter(Territory1 %in% Countries) %>% 
  st_transform(crs = 3832) #3832


# Let's try this... I think It will take 4-ever

# Shapefile_Map <- ggplot() +
#   geom_sf(data = North_America_Land) +
#   geom_sf(data = North_America_EEZ)
# 
#   # 6.26pm and 18:34 on carmelia
#   ggsave("Shapefile_Map.png",
#          plot = last_plot(),
#          width = 12,
#          height = 10,
#          units = "in"
#          )


### Explorting new dataset, hopefully is faster... ####

path_eez<- ("~/Documents/R CODE/Spatial_Analysis/Data/Intersect_EEZ_IHO_v3")
#The File
# fnam_eez <- "eez_boundaries_v10.shp" # Just lines
fnam_eez <- "EEZ_IHO_v3.shp" # Actual shapefiles (takes 4-ever, for ever ever? fooor eever eeever!?)

new_eez_sf <- st_read(dsn = path_eez, 
layer = file_path_sans_ext(fnam_eez)
) #%>% 
  # st_transform(crs = 4326) #3832 centered

# See map
# ggplot(new_eez_sf) +
#   geom_sf()

# head(new_eez_sf)

# new_eez_sf %>% 
#   filter(Territory1 == "Canada",
#     Longitude > -100) %>%
  # tail()
  # filter(Territory1 == "Canada",
         # IHO_Sea %in% c("North Atlantic Ocean")) %>% 
  # head()
  # ggplot() +
  # geom_sf()


# Canada Pacific (Half arctic sea) <- Longitude < -100

Canada_shp <- new_eez_sf %>% 
  filter(Territory1 == "Canada") %>% 
  mutate(Basin = ifelse(Longitude > -100, "Can E","Can W")) %>% 
  select(Territory1,
         Basin,
         geometry)

# head(Canada_shp) yes! 


# US Pacific (Half arctic sea) <- Longitude < -100

US_shp <- new_eez_sf %>% 
  filter(Territory1 == "United States") %>% 
  mutate(Basin = ifelse(Longitude > -100, "USA E","USA W")) %>% 
  select(Territory1,
         Basin,
         geometry)

# Alaska

Alaska_shp <- new_eez_sf %>% 
  filter(Territory1 == "Alaska") %>% 
    mutate(Basin = "Alaska") %>% 
  select(Territory1,
         Basin,
         geometry)


# Greenland

# Greenland_shp <- 
  # new_eez_sf %>% 
  # filter(Territory1 == "Greenland") %>% 
  #   mutate(Basin = "Alaska") %>% 
  # select(Territory1,
  #        Basin,
  #        geometry)
  
  North_America_EEZ_sp <- Alaska_shp %>% 
    rbind(US_shp,Canada_shp)
  

DBEM data analysis for NAFO

Areas <- c(
  "Canada (Arctic)",
  "Canada (East Coast)",
  "Canada (Pacific)",
  "USA (East Coast)",
  "USA (West Coast)",
  "USA (Alaska, Arctic)",
  "USA (Alaska, Subarctic)"
)

RCP = c("GFDL26F1","IPSL26F1","MPI26F1") #Low RCP

Selected_Areas <- EEZIDs_List %>% 
  filter(Name %in% Areas)

# DBEM_Corrected %>% 
#   filter(Treaty == "NAFO") %>% 
#   View()


Overall <- DBEM_Corrected %>%
  left_join(EEZ_CellID,
            by = "INDEX") %>% 
  # head() %>%
  filter(EEZID %in% Selected_Areas$EEZID) %>% 
  mutate( # uncomment to get overall by nation
    Nation = ifelse(EEZID >= 958, "Alaska",
                    ifelse(EEZID == 925,"Can W",
                           ifelse(EEZID == 851, "USA E",
                                  ifelse( EEZID== 848, "USA W","Can E"))
                    )),
    RCP = ifelse(Model %in% RCP,"Low","High"),
    Basin = ifelse(Nation %in% c("USA E","Can E"),"Atlantic","Pacific")
  ) %>% 
  # View()
  group_by(Model,Nation,RCP,Basin,Species) %>%
  summarise_if(is.numeric,sum,na.rm=T) %>%
  group_by(Nation,RCP,Basin,Species) %>% 
  summarise_at(vars(`2005`:`2099`),mean,na.rm=T) %>%  #Average of models
  gather("Year","MCP",`2005`:`2099`)

# Average of each Species MCP today

Today_Mean <- Overall %>% 
  filter(Year >= 2005 & Year <= 2014) %>% 
  group_by(Nation,RCP,Species) %>% 
  summarise(Today_MCP = mean(MCP))

# View(Today_Mean)

#### Change in proportions ####

# Average of each Species MCP today in both nations
Today_Total_MCP <- Overall %>% 
  filter(Year >= 2005 & Year <= 2014) %>% 
  group_by(RCP,Nation,Basin, Species) %>% 
  summarise(Mean_MCP = mean(MCP)) %>% 
  group_by(RCP,Basin,Species) %>% 
  summarise(Total_MCP = sum(Mean_MCP))


Today_Proportion <- Overall %>%
  # mutate(Basin = ifelse(Nation %in% c("USA E","Can E"),"Atlantic","Pacific")) %>% 
  filter(Year >= 2005 & Year <= 2014) %>% 
  group_by(Nation,RCP,Basin,Species) %>% 
  summarise(Today_MCP = mean(MCP)) %>% 
  left_join(Today_Total_MCP,
            by = c("RCP","Basin","Species")
  ) %>% 
  mutate(
    Prop_Today = (Today_MCP/Total_MCP)*100
  ) %>% 
  select(Nation,RCP,Basin,Prop_Today,Species)

### Check point

# Today_Proportion %>% 
#   group_by(RCP,Basin) %>% 
#   summarise(sum(Prop_Today))

### Proportion change ##

Future_Totals <- Overall %>% 
  filter(Year >= 2046) %>% 
  mutate(#Basin = ifelse(Nation %in% c("USA E","Can E"),"Atlantic","Pacific"),
         Period = (ifelse(Year >= 2046 & Year <= 2055,"Mid Century","End Century"))
) %>% 
  # filter(Year >= 2046 & Year <= 2055) %>% 
  group_by(RCP,Nation,Basin,Period,Species) %>% 
  summarise(Mean_MCP = mean(MCP)) %>% 
  group_by(RCP,Basin,Period,Species) %>% 
  summarise(Total_MCP_Fut = sum(Mean_MCP))

Result_Table <- Overall %>%
  mutate(#Basin = ifelse(Nation %in% c("USA E","Can E"),"Atlantic","Pacific"),
         Period = ifelse(Year >= 2046 & Year <= 2055,"Mid Century",
                          ifelse(Year >= 2090 & Year <= 2099,"End Century","Other_Years"))
                   ) %>% 
  filter(Period != "Other_Years") %>%
  group_by(Nation,RCP,Basin,Period,Species) %>% 
  summarise(Fut_MCP = mean(MCP)) %>% 
  # View()
  left_join(Future_Totals,
            by = c("RCP","Basin","Period","Species")
  ) %>% 
  mutate(
    Prop_Fut = (Fut_MCP/Total_MCP_Fut)*100
  ) %>% 
  left_join(Today_Proportion,
            by = c("Nation","RCP","Basin","Species")) %>% 
  mutate(Change = Prop_Fut/Prop_Today,
         Percentage_Mid = round(-(1-Change)*100)
  ) %>% 
  View()
  # select(Species,Nation:Period,Percentage_Mid) %>%
  group_by(Nation,RCP,Period) %>%
  summarise(
    Mean_All = mean(Percentage_Mid, na.rm = T),
    sd_All = sd(Percentage_Mid, na.rm = T)
  ) %>%
  View()


Result_Map <- Result_Table  %>% 
  filter(RCP == "High",
         Period == "Mid Century") %>% 
  ungroup() %>% 
  select(
    Basin= Nation,
    Mean_All
    )

  # COnvert DBEM data to a shapefile (points)
  Mid_Century_Plot <- North_America_EEZ_sp %>% 
    left_join(Result_Map,
              by ="Basin") %>% 
    st_transform(crs = 3832)
  
  # head(Mid_Century_Plot)
  # tail(Mid_Century_Plot)
  # unique(Mid_Century_Plot$)
  
  
  #### Map checking #### Works!!!!
  
  Seq <- seq(-12,12,2)
  
  Plot <- ggplot(Mid_Century_Plot) +
    geom_sf(aes(fill = Mean_All,
                colour=Mean_All)) +
    scale_fill_gradient2(
      limits=c(-12,12),
      breaks = Seq ) +
    scale_colour_gradient2(
      limits=c(-12,12),
      breaks = Seq) +
    # geom_sf(data = North_America_Land) +
    # geom_sf(data = North_America_EEZ, colour = "black") +
    ggtheme_map() +
    theme(legend.key = element_rect(size = 4),
      legend.key.width =unit(6,"line"));Plot
    
  # ggsave("Mid_Century.png", #03:24
  #        plot = Plot,
  #        width = 12,
  #        height = 10,
  #        units = "in"
  #        )
  

Chinook Map

# Run function to get all Data
suppressMessages(
  Chinook_Data <- bind_rows(
    lapply(Acronyms$Sheet,get_chinook_data)
  )
) 

NAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercionNAs introduced by coercion