tests/hw/06_mitigation.R

In ericmkeen/shipstrike: Utilities for spatially- and temporally-explicit predictions of whale-ship collisions

################################################################################
# Mitigation scenarios -- HUMPBACK WHALES
################################################################################
################################################################################

library(dplyr)
library(ggplot2)
library(devtools)
library(knitr)
document()

################################################################################
################################################################################
################################################################################
################################################################################
# (1) Speed reductions

################################################################################
# 1a	LNG-specific speed reduction (7 – 9 kn)

# LNG Canada ===================================================================

# Encounter rate
data(lng_canada)
encounter_rate(vessels = lng_canada,
               whales = humpback_params(),
               outcome_dir = 'tests/hw/mitigation/1a/lng_canada/',
               month_batches = list(all = 1:12),
               new_speeds = runif(100, 7, 9),
               runs  = 100, iterations = 100, toplot = TRUE)

# Outcomes
bootstraps <- readRDS('tests/hw/pwhale_seasonal_boots.RData')
data(lng_canada)
(vessels <- unique(lng_canada$type))
data(p_surface)

data(p_collision) ; (avoidance <- p_collision)
data(p_lethality) ; (lethality <- p_lethality)
(avoidance <- avoidance[c(2,4),])
(avoidance$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))
(lethality <- lethality[c(2,4),])
(lethality$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))

lng_canada_mod <- lng_canada
lng_canada_mod$speed <- runif(nrow(lng_canada_mod), 7, 9)

outcome_predict(traffic = lng_canada_mod,
                scale_factors=NULL,
                whale = bootstraps,
                seasonal = NULL,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/hw/mitigation/1a/lng_canada/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

# Cedar LNG ====================================================================

# Encounter rate
data(cedar_lng)
encounter_rate(vessels = cedar_lng,
               whales = humpback_params(),
               outcome_dir = 'tests/hw/mitigation/1a/cedar_lng/',
               month_batches = list(all = 1:12),
               new_speeds = runif(100, 7, 9),
               runs  = 100, iterations = 100, toplot = TRUE)

# Outcomes
bootstraps <- readRDS('tests/hw/pwhale_seasonal_boots.RData')
data(cedar_lng)
(vessels <- unique(cedar_lng$type))
data(p_surface)
data(p_collision) ; (avoidance <- p_collision)
data(p_lethality) ; (lethality <- p_lethality)
(avoidance <- avoidance[c(2,4),])
(avoidance$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))
(lethality <- lethality[c(2,4),])
(lethality$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))

cedar_lng_mod <- cedar_lng
cedar_lng_mod$speed <- runif(nrow(cedar_lng_mod), 7, 9)

outcome_predict(traffic = cedar_lng_mod,
                scale_factors=NULL,
                whale = bootstraps,
                seasonal = NULL,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/hw/mitigation/1a/cedar_lng/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

################################################################################
# 1b	Speed reduction for all traffic > 100m (7 – 9 kn for LNG, < 10 kn for all other)

# Encounter rate: AIS 2030
data(ais_2019)
encounter_rate(vessels = ais_2019,
               whales = humpback_params(),
               outcome_dir = 'tests/hw/mitigation/1b/ais_2030/',
               month_batches = list(all = 1:12),
               speed_restriction = 9,
               lengths_restricted = 180,
               runs  = 100, iterations = 100, toplot = TRUE)

# Predict outcomes on AIS 2030 with < 9 kn
bootstraps <- readRDS('tests/hw/pwhale_seasonal_boots.RData')
data(p_surface)
data(ais_2019)
scale_factors <- readRDS('tests/ais/vessel_trends.RData')
data(p_collision)
data(p_lethality)

ais_2019_mod <- ais_2019
bads <- which(ais_2019_mod$length > 180 & ais_2019_mod$speed > 9)
ais_2019_mod$speed[bads] <- 9

outcome_predict(traffic = ais_2019_mod,
                scale_factors=NULL,
                whale = bootstraps,
                seasonal = NULL,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/hw/mitigation/1b/ais_2030/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

################################################################################
################################################################################
################################################################################
################################################################################
# (2) Daytime-only transits

################################################################################
# 2a	LNG-only

# alter LNG traffic to day, then re-run predict outcomes
# do so for LNG Canada
# do so for Cedar LNG

# LNG Canada ===================================================================

bootstraps <- readRDS('tests/hw/pwhale_seasonal_boots.RData')
data(p_surface)
data(lng_canada)
(vessels <- unique(lng_canada$type))

lng_canada_mod <- lng_canada # Modify diel here
lng_canada_mod$diel <- 'day'

data(p_collision) ; (avoidance <- p_collision)
data(p_lethality) ; (lethality <- p_lethality)
(avoidance <- avoidance[c(2,4),])
(avoidance$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))
(lethality <- lethality[c(2,4),])
(lethality$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))

outcome_predict(traffic = lng_canada_mod,
                scale_factors=NULL,
                whale = bootstraps,
                seasonal = NULL,
                p_encounter_dir = 'tests/hw/lng_canada/',
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/hw/mitigation/2a/lng_canada/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

# Cedar LNG ====================================================================

bootstraps <- readRDS('tests/hw/pwhale_seasonal_boots.RData')
data(p_surface)
data(cedar_lng)
(vessels <- unique(cedar_lng$type))

cedar_lng_mod <- cedar_lng # Modify diel here
cedar_lng_mod$diel <- 'day'

data(p_collision) ; (avoidance <- p_collision)
data(p_lethality) ; (lethality <- p_lethality)
(avoidance <- avoidance[c(2,4),])
(avoidance$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))
(lethality <- lethality[c(2,4),])
(lethality$type <- c(paste(vessels[3:4], collapse=' | '),
                     paste(vessels[1:2], collapse=' | ')))

outcome_predict(traffic = cedar_lng_mod,
                scale_factors=NULL,
                whale = bootstraps,
                seasonal = NULL,
                p_encounter_dir = 'tests/hw/cedar_lng/',
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/hw/mitigation/2a/cedar_lng/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

################################################################################
# 2b	All traffic > 180m)

bootstraps <- readRDS('tests/hw/pwhale_seasonal_boots.RData')
seasonal_boot <- NULL
scale_factors <- readRDS('tests/ais/vessel_trends.RData')
data(p_collision)
data(p_lethality)
data(p_surface)
data(ais_2019)

ais_mod <- ais_2019 # Modify diel here
bads <- which(ais_mod$length > 180)
ais_mod$diel[bads] <- 'day'

outcome_predict(traffic = ais_mod,
                scale_factors=scale_factors,
                whale = bootstraps,
                seasonal = seasonal_boot,
                p_encounter_dir = 'tests/hw/ais_2019/',
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/hw/mitigation/2b/ais_2030/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                iterations = 1000)

################################################################################
################################################################################
################################################################################
################################################################################
# (3) Seasonal rescheduling of LNG traffic

if(TRUE){

  # Baseline scenario (no mitigation)   ==========================================

  results <- readRDS('tests/hw/results.RData')
  outcomes <- rbind(results$outcomes$lng_canda,
                    results$outcomes$cedar_lng)

  # 3a	One month  ===============================================================
  mr <- mitigate_loop(outcomes,
                      mitigation_duration = 1,
                      reschedule = TRUE)
  saveRDS(mr, file='tests/hw/mitigation/3a/mitigation.RData')

  # 3b	Two months  ==============================================================
  mr <- mitigate_loop(outcomes,
                      mitigation_duration = 2,
                      reschedule = TRUE)
  saveRDS(mr, file='tests/hw/mitigation/3b/mitigation.RData')

  # 3c	Three months  ============================================================
  mr <- mitigate_loop(outcomes,
                      mitigation_duration = 3,
                      reschedule = TRUE)
  saveRDS(mr, file='tests/hw/mitigation/3c/mitigation.RData')

  ################################################################################
  ################################################################################
  ################################################################################
  ################################################################################
  # (4) Seasonal moratorium on LNG traffic

  # 4a	One month ================================================================
  mr <- mitigate_loop(outcomes,
                      mitigation_duration = 1,
                      reschedule = FALSE)
  saveRDS(mr, file='tests/hw/mitigation/4a/mitigation.RData')

  # 4b	Two months  ==============================================================
  mr <- mitigate_loop(outcomes,
                      mitigation_duration = 2,
                      reschedule = FALSE)
  saveRDS(mr, file='tests/hw/mitigation/4b/mitigation.RData')

  # 4c	Three months  ============================================================
  mr <- mitigate_loop(outcomes,
                      mitigation_duration = 3,
                      reschedule = FALSE)
  saveRDS(mr, file='tests/hw/mitigation/4c/mitigation.RData')

}

################################################################################
################################################################################
# Put all outcomes in a list

mitigations <- list(m1 = list(lng_canada = list(p_encounter = readRDS('tests/hw/mitigation/1a/lng_canada/p_encounter.RData'),
                                                outcomes = readRDS('tests/hw/mitigation/1a/lng_canada/outcomes.RData'),
                                                grid = readRDS('tests/hw/mitigation/1a/lng_canada/outcomes_grid.RData')),
                              cedar_lng = list(p_encounter = readRDS('tests/hw/mitigation/1a/cedar_lng/p_encounter.RData'),
                                               outcomes = readRDS('tests/hw/mitigation/1a/cedar_lng/outcomes.RData'),
                                               grid = readRDS('tests/hw/mitigation/1a/cedar_lng/outcomes_grid.RData')),
                              ais_2030 = list(p_encounter = readRDS('tests/hw/mitigation/1b/ais_2030/p_encounter.RData'),
                                              outcomes = readRDS('tests/hw/mitigation/1b/ais_2030/outcomes.RData'),
                                              grid = readRDS('tests/hw/mitigation/1b/ais_2030/outcomes_grid.RData'))),
                    m2 = list(lng_canada = list(outcomes = readRDS('tests/hw/mitigation/2a/lng_canada/outcomes.RData'),
                                                grid = readRDS('tests/hw/mitigation/2a/lng_canada/outcomes_grid.RData')),
                              cedar_lng = list(outcomes = readRDS('tests/hw/mitigation/2a/cedar_lng/outcomes.RData'),
                                               grid = readRDS('tests/hw/mitigation/2a/cedar_lng/outcomes_grid.RData')),
                              ais_2030 = list(outcomes = readRDS('tests/hw/mitigation/2b/ais_2030/outcomes.RData'),
                                              grid = readRDS('tests/hw/mitigation/2b/ais_2030/outcomes_grid.RData'))),
                    m3 = list(one_month = readRDS('tests/hw/mitigation/3a/mitigation.RData'),
                              two_months = readRDS('tests/hw/mitigation/3b/mitigation.RData'),
                              three_months = readRDS('tests/hw/mitigation/3c/mitigation.RData')),
                    m4 = list(one_month = readRDS('tests/hw/mitigation/4a/mitigation.RData'),
                              two_months = readRDS('tests/hw/mitigation/4b/mitigation.RData'),
                              three_months = readRDS('tests/hw/mitigation/4c/mitigation.RData')))

saveRDS(mitigations, file='tests/hw/mitigation/mitigation.RData')