shipstrike

Eric Keen

Science co-director, North Coast Ceteacean Society

Sewanee: The University of the South

Introduction

The R package shipstrike offers utilities for spatially- and temporally-explicit predictions of whale-ship interactions and collisions. Its chief asset is its modular design, allowing researchers to parse the ship-strike analysis process - which can be complex – into discrete stages while still propagating uncertainty from every stage into the final result.

The package was developed into two stages: first, in association with this paper, which focused on the close-encounter rate simulation tools in shipstrike. The vignette originally published in tandem with that paper can be found here. We now publish the second version of the package in association with a forthcoming publication predicting ship-strike rates in the territorial waters of the Gitga’at First Nation. The link to that study – which we shall refer to as Keen et al. (2023) – will be provided once it is published.

The shipstrike framework involves the following analytical steps:

  1. Cooccurrences: To predict whale-vessel interaction rates, first prepare spatial grids of vessel traffic, density surface models of whale density, and seasonal models of whale abundance to predict the number of times in which a vessel and a whale occur in the same square km (hereafter, a ‘cooccurrence’).

  2. Close encounters: Second, cooccurrences are scaled by a ‘close-encounter’ rate, i.e., the rate at which a vessel and whale are expected to intersect in time and horizontal space assuming no avoidance measures are taken. The close-encounter rate is estimated using simulations that draw upon size and travel patterns for both actors.

  3. Strike-zone events: Third, close encounters are scaled according to vessel draft and whale depth distribution, which we inferred using whale-born time-depth-recording tags, to estimate the number of ‘strike-zone events’, i.e., the times in which the whale and vessel will collide unless avoidance measures are taken.

  4. Collisions: To predict collisions, strike-zone events are scaled according to an avoidance rate that was modeled as a function of vessel speed (Gende et al. 2011, Rockwood et al. 2017).

  5. Mortality: Finally, the number of collisions is scaled by the lethality rate, which we also model using equations from the ship-strike literature (Kelley et al. 2020), to predict the number of whale mortalities.

This analysis requires a variety of datasets, and most of the work ends up being in the preparation of those datasets. The general workflow for preparing and implementing a ship-strike analysis with shipstrike looks something like this:

  1. Prepare a spatial grid.

  2. Prepare marine traffic data.

  3. Prepare whale density surface.

  4. Prepare whale seasonal abundance curve.

  5. Estimate close-encounter rates.

  6. Predict whale-ship interaction rates & vessel strikes.

Here we demonstrate that workflow and the features of shipstrike using Keen et al. (2023) as a case study. The code below replicates the analysis therein.

Installation

The package can be installed directly from Github:

# Make sure you have "devtools" installed
if (!require('devtools')) install.packages('devtools')

# Install from GitHub
devtools::install_github('ericmkeen/shipstrike')

Now load the package:

library(shipstrike)
## Writing NAMESPACE
## Writing NAMESPACE
## Writing ais_table.Rd
## Writing vessel_grid.Rd

Other dependencies

The analysis from Keen et al. (2023) will rely on other packages as well:

library(readr)
library(dplyr)
library(devtools)
library(lubridate)
library(ggplot2)
library(knitr)

The analysis will also rely upon utilities specific to its study area, the Kitimat Fjord System (British Columbia, Canada), which are provided in the R package bangarang (also available on Github):

devtools::install_github('ericmkeen/bangarang')
library(bangarang)

Spatial grid

In shipstrike, the foundation of your analysis will be a spatial grid of your study area.

The 1-km2 grid used in Keen et al. (2023) is provided as a built-in package dataset:

data(grid_kfs, package='shipstrike')
head(grid_kfs)
##         y1      y2        y        x1        x2         x      km2  id
## 1 52.79656 52.8082 52.80238 -129.3193 -129.3077 -129.3135 1.016831 466
## 2 52.79656 52.8082 52.80238 -129.3077 -129.2961 -129.3019 1.016831 467
## 3 52.79656 52.8082 52.80238 -129.2146 -129.2029 -129.2088 1.016831 475
## 4 52.79656 52.8082 52.80238 -129.2029 -129.1913 -129.1971 1.016831 476
## 5 52.79656 52.8082 52.80238 -129.1913 -129.1797 -129.1855 1.016831 477
## 6 52.79656 52.8082 52.80238 -129.1797 -129.1680 -129.1738 1.016831 478
##           z       zmin     zmax       zsd    zrange block
## 1 163.81329  0.0000005 209.7133 59.861445 209.71330   CAA
## 2 223.28670 12.4876499 237.9200 44.017913 225.43235   CAA
## 3  44.95333  0.0000000  79.5200 21.548575  79.52000   CAA
## 4  68.68333 45.8133316  84.8300  6.757624  39.01667   CAA
## 5  77.22667 69.2300034  95.6100  6.129383  26.38000   CAA
## 6  97.62000 81.7600021 125.8667 11.170772  44.10670   CAA

In this data.frame, each row is a grid cell. The essential fields are y1 (lower bound of cell, in decimal degrees), y2 (right bound), y (latitudinal center), x1 (left bound), x2 (right bound), x (longitudinal center), km2 (the area of the cell, in square km), and id (a unique identifier). The block field is also necessary if your study area is split into distinct regions that you want to analyze separately. In our case, we wanted to investigate risk differences across 8 different waterways in the Kitimat Fjord System. The other z... fields are only necessary if needed for density surface modeling, which we will not cover here.

The grid_kfs object was produced using the shipstrike function make_grid():

grid_kfs <- make_grid(xlims = c(-129.75, -128.5),
                      ylims = c(52.75, 54.1),
                      grid_int = 1.287)

This function produces a data.frame for a rectangular grid of cells.

bangarang::gg_kfs() + 
  geom_point(data=grid_kfs, mapping=aes(x=x, y=y), size=.5)

Any cell that occurs on land (mainland or islands) needs to be removed. We did so using the bangarang functions inKFS() (determines whether a coordinate is inside the Kitimat Fjord System proper) and inwater() (determines whether a coordinate is on land or water) as well as some bangarang built-in datasets.

# Load datasets from `bangarang` pkg
data(kfs_boundary)
kfs <- kfs_boundary
data(kfs_land)
land <- kfs_land

# Perform tests
newgrids <- data.frame(x = grid_kfs$x, y = grid_kfs$y)
newgrids$inkfs <- inKFS(newgrids,kfs,toplot=TRUE)$inkfs
newgrids$inwater <- inwater(newgrids,land,toplot=TRUE)$valid

# Update grid_kfs
grid_kfs <-
  newgrids %>%
  dplyr::filter(inkfs == TRUE,
                inwater == TRUE)

Our spatial grid now looks like this:

bangarang::gg_kfs() + 
  geom_point(data=grid_kfs, mapping=aes(x=x, y=y), size=.5)

Marine traffic

In Keen et al. (2023), we use both AIS observations of marine traffic as well as simulated marine traffic.

AIS data

The chief AIS dataset in Keen et al. (2023) is from 2019. The final polished version of that dataset is provided as a built-in dataset within shipstrike:

data(ais_2019)

head(ais_2019)
## # A tibble: 6 × 17
## # Groups:   grid_id [2]
##   grid_id   vid type    speed length width draft datetime                x     y
##     <int> <dbl> <chr>   <dbl>  <dbl> <dbl> <dbl> <dttm>              <dbl> <dbl>
## 1    5713   643 Other …   4.5     50    11     5 2019-01-02 02:24:28 -129.  53.4
## 2    5713   643 Other …   8.3     50    11     5 2019-01-02 02:24:28 -129.  53.4
## 3    5713   643 Other …   8.6     50    11     5 2019-01-02 02:24:28 -129.  53.4
## 4    5606   643 Other …   8.6     50    11     5 2019-01-02 02:29:08 -129.  53.4
## 5    5606   643 Other …   9.3     50    11     5 2019-01-02 02:29:08 -129.  53.4
## 6    5606   643 Other …   6.6     50    11     5 2019-01-02 02:29:08 -129.  53.4
## # … with 7 more variables: km <dbl>, year <dbl>, month <dbl>, yday <dbl>,
## #   channel <chr>, sun <dbl>, diel <chr>

This is what your AIS data will need to look like in order to proceed with using shipstrike. We prepared our AIS data using the following code.

Read, filter, format, & re-group

Read in the raw AIS data:

fn <- '/Users/erickeen/Dropbox/Other WIPs/NCCS/Projects/ship-strike/data/ais/AIS_2019_w_month.csv'
df <- read_csv(fn)

head(df)
## # A tibble: 6 × 17
##    ...1     X UTC.Time    Local…¹ Type  Length  Beam Draught   SOG   COG Latit…²
##   <dbl> <dbl> <chr>       <chr>   <chr>  <dbl> <dbl>   <dbl> <dbl> <dbl>   <dbl>
## 1     1     1 2019-01-01… 2019-0… Tug       17     6       0     0  334.    54.0
## 2     2     2 2019-01-01… 2019-0… Tug       17     6       0     0  230.    54.0
## 3     3     3 2019-01-01… 2019-0… Tug       17     6       0     0  282     54.0
## 4     4     4 2019-01-01… 2019-0… Tug       17     6       0     0  312.    54.0
## 5     5     5 2019-01-01… 2019-0… Tug       17     6       0     0  210.    54.0
## 6     6     6 2019-01-01… 2019-0… Tug       17     6       0     0  301     54.0
## # … with 6 more variables: Longitude <dbl>, Destination <chr>, Status <chr>,
## #   Flag <chr>, id <dbl>, Month.lt <dbl>, and abbreviated variable names
## #   ¹​Local.Time, ²​Latitude

Now filter out obviously erroneous and irrelevant entries, and reformat column names:

ais <-
  df %>%
  dplyr::filter(SOG > 3,
                SOG < 40,
                Length > 5,
                Length < 500) %>%
  dplyr::filter(Longitude >= -129.68,
                Longitude < -128.66,
                Latitude >= 52.8,
                Latitude < 53.55) %>%
  dplyr::select(vid = id,
                type = Type,
                speed = SOG,
                length = Length,
                width = Beam,
                draft = Draught,
                datetime = UTC.Time,
                x = Longitude,
                y = Latitude) %>%
  dplyr::mutate(datetime = lubridate::ymd_hm(datetime, tz='UTC'))

# Fix obviously erroneous beam widths 
# using a length conversion factor of 0.125
bads <- which(ais$width < 2)
ais$width[bads] <- NA
bads <- which(is.na(ais$width))
if(length(bads)>0){ais$width[bads] <- ais$length[bads] * 0.125}

# Find invalid drafts
# any draft 50% of length or deeper is invalid
ld_ratio <- ais$draft / ais$length
bads <- which(ld_ratio > 0.5)
if(length(bads)>0){ais$draft[bads] <- NA}

# Any draft less than 0.2m is also invalid
bads <- which(ais$draft < .2)
if(length(bads)>0){ais$draft[bads] <- NA}

# Now fix invalid drafts
# Use a 1-to-0.05 ratio as the conversion factor
bads <- which(is.na(ais$draft))
if(length(bads)>0){ais$draft[bads] <- ais$length[bads] * 0.05}

Note that AIS data rarely follows the exact same formatting, so you will likely need to modify the code above in order to achieve the desired outcome with your own AIS data.

Check out the AIS data sample size:

# Number of records
nrow(ais)
## [1] 229452
# Number of unique vessel IDs
ais$vid %>% unique %>% length
## [1] 992

This dataset contains over 20 vessel types:

(types <- ais$type %>% unique)
##  [1] "Search/rescue"              "Passenger ship"            
##  [3] "Tug"                        "Towing"                    
##  [5] "Cargo ship"                 "Fishing"                   
##  [7] "Towing(200/25)"             "Dredging or underwater op."
##  [9] "Pleasure Craft"             "Other"                     
## [11] "Pilot"                      "Law enforcement"           
## [13] "Sailing"                    "Tanker"                    
## [15] "HSC"                        "Diving op."                
## [17] "Passenger ship:DG,HS,MP(Y)" "Local ship"                
## [19] "Port tender"                "Military op."              
## [21] "Cargo ship:DG,HS,MP(X)"

In Keen et al. (2023) we pool them into just 10 types:

newtype <- rep('Other < 40m',times=nrow(ais))
newtype[ais$length > 40] <- 'Other > 40m'
newtype[ais$length > 100] <- 'Other > 100m'
newtype[ais$type %in% types[c(9)] & ais$length < 40] <- 'Pleasurecraft < 40m'
newtype[ais$type %in% types[c(4, 7)] & ais$length < 50] <- 'Towing < 50m'
newtype[ais$type %in% types[c(3)] & ais$length < 50] <- 'Tug < 50m'
newtype[ais$type %in% types[c(2, 17)] & ais$length >= 180] <- 'Passenger > 180m'
newtype[ais$type %in% types[c(5, 21)] & ais$length >= 180] <- 'Cargo > 180m'
newtype[ais$type %in% types[c(13)]] <- 'Sailing'
newtype[ais$type %in% types[c(6)] & ais$length < 60] <- 'Fishing < 60m'

# Update with new type categories
ais$type <- newtype

New vessel type names:

ais$type %>% unique
##  [1] "Other > 40m"         "Other > 100m"        "Tug < 50m"          
##  [4] "Towing < 50m"        "Cargo > 180m"        "Fishing < 60m"      
##  [7] "Other < 40m"         "Pleasurecraft < 40m" "Sailing"            
## [10] "Passenger > 180m"

Now summarize the AIS dataset by vessel type, using the shipstrike function ais_table():

ais_table(ais) 
## # A tibble: 10 × 22
##    type       vids fixes trans…¹ trans…² dates…³ prese…⁴ speed…⁵ speed…⁶ speed…⁷
##    <chr>     <int> <int>   <dbl>   <dbl>   <int>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 Cargo > …    38  4480      94    0.26      88    0.24    13.1     1.4    17.1
##  2 Fishing …   305 41720     822    2.25     282    0.77     8.4     3.2    32.2
##  3 Other < …    70 15106     565    1.55     296    0.81    11.5     7.2    39.2
##  4 Other > …    23 17155     378    1.04     268    0.73    16       3.4    23  
##  5 Other > …    46 15551     340    0.93     229    0.63     9.5     2.6    16.6
##  6 Passenge…     6  2845      73    0.2       65    0.18    17.4     3.4    23.3
##  7 Pleasure…   263 42416    1123    3.08     193    0.53     8.1     3.8    37.5
##  8 Sailing     117 19520     426    1.17     150    0.41     6       1.3    12.6
##  9 Towing <…    74 33783     738    2.02     322    0.88     8       1.9    21.7
## 10 Tug < 50m    61 36876     835    2.29     328    0.9      7.3     1.8    13.9
## # … with 12 more variables: length_mn <dbl>, length_sd <dbl>, length_min <dbl>,
## #   length_max <dbl>, beam_mn <dbl>, beam_sd <dbl>, beam_min <dbl>,
## #   beam_max <dbl>, draft_mn <dbl>, draft_sd <dbl>, draft_min <dbl>,
## #   draft_max <dbl>, and abbreviated variable names ¹​transits, ²​transit_rate,
## #   ³​dates_present, ⁴​present_rate, ⁵​speed_mn, ⁶​speed_sd, ⁷​speed_max

Interpolate the vessel grid

Now that the data are clean and re-grouped, we can proceed with data processing. First we interpolate vessel position fixes and associate each interpolated fix to a grid cell, using a function from shipstrike. In the result, each row is an interpolated vessel position fix. This process can take several minutes.

vgrid <- vessel_grid(grid_kfs, ais)

We then used the grid cell to assign each vessel position fix to a waterway in our study area, by joining to our spatial grid:

grid_join <- grid_kfs %>% select(grid_id = id, channel = block)
vgrid2 <- left_join(vgrid, grid_join, by='grid_id')

Determine sun angle

Now add columns for the elevation of the sun during each record, using a shipstrike function:

vgrid2$sun <- vessel_sun_angle(vgrid2, verbose=TRUE)

We can then use solar elevation to determine which events occur during day (sun >= -12 degrees, which is nautical dawn/dusk) and which during night (sun < -12 degrees).

vgrid2$diel <- 'day'
vgrid2$diel[vgrid2$sun < -12] <- 'night' # nautical dawn/dusk

Your AIS dataset is now ready for use in subsequent shipstrike functions. Again, note that the package provides a built-in version of this dataset:

data(ais_2019)

# Map it
gg_kfs() + 
  geom_point(data= ais_2019 %>% group_by(grid_id) %>% 
               summarize(n=n(), x=x[1], y=y[1]),
             mapping=aes(x=x,y=y,color=sqrt(n))) + 
  scale_colour_gradient(low='white', high='firebrick')

We use this exact same code to produce datasets for AIS data in 2014, 2015, and 2018.

data(ais_2014)
data(ais_2015)
data(ais_2018)

Simulated traffic

In Keen et al. (2023), we simulate traffic from two LNG-trafficking projects: (1) LNG Canada and (2) Cedar LNG. Both of these shipping projects use the exact same route:

data(tanker_route) # from `bangarang` package
tanker_route %>% head
##           x        y
## 1 -128.6883 53.99415
## 2 -128.6878 53.99239
## 3 -128.6873 53.99062
## 4 -128.6868 53.98885
## 5 -128.6862 53.98709
## 6 -128.6857 53.98532
bangarang::gg_kfs() +
  geom_point(data=tanker_route, mapping=aes(x=x, y=y))
## Warning: Removed 512 rows containing missing values (geom_point).

See the Keen et al. study for all references.

LNG Canada

The LNG Canada project will involve 350 ship calls to the Port of Kitimat:

new_calls <- 350

The dimensions of these ships will vary slightly throughout their fleet. The possibilities are provided in this data.frame:

(vessel_options <- data.frame(length = c(290, 298, 286, 288, 286),
                              width = c(45, round((298/315)*50), 43, 44, 41),
                              draft = c(12, 12, 11.9, 11.5, 11.8)))

These ship calls will be evenly distributed throughout the year, and the exit transit (‘in product’) will always occur one day after the entry transit (‘in heel’):

(julians <- sample(1:365, size=new_calls, replace=FALSE) %>% sort)
julian_product <- julians + 1

First we will simulate transits for the in-heel entry transits:

# Arrange tanker route coordinates in reverse order
traffic <-
    tanker_route %>%
    select(x, y) %>%
    mutate(id = 1:n()) %>%
    arrange(desc(id))

# Loop through each transit
heels <- data.frame()
t0 <- lubridate::as_datetime('2030-01-01 00:00:01', tz='UTC')
for(i in 1:length(julians)){
  traffi <- traffic
  (start_time <- t0 + lubridate::days(julians[i]) + lubridate::seconds(sample(1:86400, 1)))

  # Specify speed along route, randomly drawn from a uniform distribution
  (speedi <- runif(1,8,14))
  traffi$speed <- speedi

  # Add slight course variability
  (delta_y <- rnorm(1, 0, .002))
  traffi$y <- traffi$y + delta_y
  (delta_x <- rnorm(1, 0, .002))
  traffi$x <- traffi$x + delta_x

  # Compile a dataframe of position fixes
  (xy <- data.frame(lat1 = traffi$y[1:(nrow(traffi)-1)],
                    lon1 = traffi$x[1:(nrow(traffi)-1)],
                    lat2 = traffi$y[2:(nrow(traffi))],
                    lon2 = traffi$x[2:(nrow(traffi))]))
  
  # Determine the distance between each position fix
  (nmi <- apply(xy, 1, function(x){swfscMisc::distance(x[1],x[2],x[3],x[4], units='nm')}))
  (kms <- nmi * 1.852)
  kms <- c(kms,0.1)
  traffi$km <- kms
  
  # Simulate timestamps for each fix
  (traffi$hours <- traffi$km / (traffi$speed * 1.852))
  (traffi$tot_hours <- cumsum(traffi$hours))
  traffi$datetime <- start_time + lubridate::seconds(round(3600*traffi$tot_hours))
  traffi$datetime
  head(traffi)

  # Select LNG ship dimension
  (vessi <- vessel_options[sample(1:nrow(vessel_options), size = 1),])

  # Format to match AIS data expectation
  ais <-
    traffi %>%
    mutate(type = 'LNG Canada tanker in-heel',
           length = vessi$length,
           width = vessi$width,
           draft = vessi$draft,
           vid = i) %>%
    dplyr::filter(x >= -129.68,
                  x < -128.66,
                  y >= 52.8,
                  y < 53.55) %>%
    dplyr::select(vid,
                  type,
                  speed,
                  length,
                  width,
                  draft,
                  datetime,
                  x,
                  y)

  heels <- rbind(heels, ais)
  message(i)
}

Now do the same for in-product transits:

traffic <-
   tanker_route %>%
   select(x, y, speed) %>%
   mutate(id = 1:n())) 

products <- data.frame()
t0 <- lubridate::as_datetime('2030-01-01 00:00:01', tz='UTC')
for(i in 1:length(julian_product)){
  traffi <- traffic
  (start_time <- t0 + lubridate::days(julian_product[i]) + lubridate::seconds(sample(1:86400, 1)))

  (speedi <- runif(1,8,14))
  traffi$speed <- speedi

  (delta_y <- rnorm(1, 0, .002))
  traffi$y <- traffi$y + delta_y
  (delta_x <- rnorm(1, 0, .002))
  traffi$x <- traffi$x + delta_x

  (xy <- data.frame(lat1 = traffi$y[1:(nrow(traffi)-1)],
                    lon1 = traffi$x[1:(nrow(traffi)-1)],
                    lat2 = traffi$y[2:(nrow(traffi))],
                    lon2 = traffi$x[2:(nrow(traffi))]))
  
  (nmi <- apply(xy, 1, function(x){swfscMisc::distance(x[1],x[2],x[3],x[4], units='nm')}))
  (kms <- nmi * 1.852)
  kms <- c(kms,0.1)
  traffi$km <- kms
  
  (traffi$hours <- traffi$km / (traffi$speed * 1.852))
  (traffi$tot_hours <- cumsum(traffi$hours))
  traffi$datetime <- start_time + lubridate::seconds(round(3600*traffi$tot_hours))
  traffi$datetime
  
  # Get dimensions based on corresponding in-heel transit
  (heeli <- heels[heels$vid == i, 4:6])
  vessi <- heeli[!duplicated(heeli),]

  ais <-
    traffi %>%
    mutate(type = 'LNG Canada tanker in-product',
           length = vessi$length,
           width = vessi$width,
           draft = vessi$draft,
           vid = i) %>%
    dplyr::filter(x >= -129.68,
                  x < -128.66,
                  y >= 52.8,
                  y < 53.55) %>%
    dplyr::select(vid,
                  type,
                  speed,
                  length,
                  width,
                  draft,
                  datetime,
                  x,
                  y)
  ais %>% head
  products <- rbind(products, ais)
  message(i)
}

Now we account for the escort tugs that will accompany these ships:

tug_length <- 35
4/27 # width:length ratio from AIS 2019 tugs
tug_width <- (4 / 27)* 35
tug_draft <- 3.1 # mead draft from AIS 2019 tugs

tug_heels <- heels
tug_heels$type <- 'LNG Canada tug in-heel'
tug_heels$length <- tug_length
tug_heels$width <- tug_width
tug_heels$draft <- tug_draft

tug_product <- products
tug_product$type <- 'LNG Canada tug in-product'
tug_product$length <- tug_length
tug_product$width <- tug_width
tug_product$draft <- tug_draft

We combine all transits into the same ais object:

ais <- rbind(heels, products, tug_heels, tug_product)

Now we proceed with the exact same processing steps as used with actual AIS data:

# Interpolate vessel grid
vgrid <- vessel_grid(grid_kfs, ais, verbose=TRUE)

# Join channel to each grid number
grid_join <- grid_kfs %>% select(grid_id = id, channel = block)
vgrid2 <- left_join(vgrid, grid_join, by='grid_id')

# Add sun angles
vgrid2$sun <- vessel_sun_angle(vgrid2, verbose=TRUE)
vgrid2$diel <- 'day'
vgrid2$diel[vgrid2$sun < -12] <- 'night' # nautical dawn/dusk

We provide this simulated LNG Canada dataset as a built-in dataset in shipstrike:

data(lng_canada)

# Map it
gg_kfs() + 
  geom_point(data= lng_canada %>% group_by(grid_id) %>% 
               summarize(n=n(), x=x[1], y=y[1]),
             mapping=aes(x=x,y=y,color=sqrt(n))) + 
  scale_colour_gradient(low='white', high='firebrick')

Cedar LNG

We perform the exact same procedure for the Cedar LNG project, the only differences being that we expect only 50 ship calls regularly spaced throughout the year for this project.

# Total tankers expected
new_calls <- 50

# Get calendar dates of transits
set.seed(126)
(spacing <- sample(5:9, size=new_calls, replace=TRUE))
(julians <- cumsum(spacing))
julian_product <- julians + 1

We also provide this simulated Cedar LNG traffic as a built-in dataset:

data(cedar_lng)

Fin whales

Keen et al. (2023) predicts whale-ship interactions for fin whales first, then humpback whales second.

Whale density surface

Keen et al. (2023) use a density surface model to prepare a density surface for fin whales. That dataset looks like this:

load('../tests/fw/dsm-estimate.RData')
grids %>% head
##         y1      y2        y        x1        x2         x      km2  id
## 1 52.79656 52.8082 52.80238 -129.3193 -129.3077 -129.3135 1.016831 466
## 2 52.79656 52.8082 52.80238 -129.3077 -129.2961 -129.3019 1.016831 467
## 3 52.79656 52.8082 52.80238 -129.2146 -129.2029 -129.2088 1.016831 475
## 4 52.79656 52.8082 52.80238 -129.2029 -129.1913 -129.1971 1.016831 476
## 5 52.79656 52.8082 52.80238 -129.1913 -129.1797 -129.1855 1.016831 477
## 6 52.79656 52.8082 52.80238 -129.1797 -129.1680 -129.1738 1.016831 478
##           z       zmin     zmax       zsd    zrange block     d_mean
## 1 163.81329  0.0000005 209.7133 59.861445 209.71330   CAA 0.02475727
## 2 223.28670 12.4876499 237.9200 44.017913 225.43235   CAA 0.02475727
## 3  44.95333  0.0000000  79.5200 21.548575  79.52000   CAA 0.02072076
## 4  68.68333 45.8133316  84.8300  6.757624  39.01667   CAA 0.02849723
## 5  77.22667 69.2300034  95.6100  6.129383  26.38000   CAA 0.03504774
## 6  97.62000 81.7600021 125.8667 11.170772  44.10670   CAA 0.03268333

Essentially, it is the data(grid_kfs) dataset with a single new column: d_mean, providing the season-wide mean density expected for each respective grid cell.

gg_kfs() + 
  geom_point(data=grids, mapping=aes(x=x, y=y, color=d_mean)) + 
  scale_colour_gradient(low='white', high='firebrick')

Keen et al. (2023) uses a bootstrapping routine to develop uncertainty distributions for each grid cell estimate:

load('../tests/fw/dsm-bootstraps.RData')
bootstraps %>% head
##   grid_id iteration month            D
## 1     466       711     0 0.0025555857
## 2     467       711     0 0.0032540102
## 3     475       711     0 0.0017659839
## 4     476       711     0 0.0017545198
## 5     477       711     0 0.0017320005
## 6     478       711     0 0.0001134392

This dataset has 1,000 estimates of D for each grid_id. Note that month is set to 0 because this density estimate for the entire summer season.

Whale seasonal abundance curve

Keen et al. (2023) uses a seasonal abundance model to scale the fin whale density surface, which shows average density from June to September, for each month of the year.

load('../tests/fw/seasonal_posterior.RData')
seasonal_boot %>% head(12)
##    iteration month doy      scaled
## 1          2     1  16 0.004268919
## 2          2     2  46 0.012001038
## 3          2     3  76 0.019733158
## 4          2     4 106 0.098665789
## 5          2     5 136 0.493328947
## 6          2     6 166 0.843473288
## 7          2     7 196 1.054233576
## 8          2     8 226 1.088618045
## 9          2     9 256 0.910483165
## 10         2    10 286 0.533614878
## 11         2    11 316 0.106722976
## 12         2    12 345 0.021344595
# Plot it
ggplot(seasonal_boot, 
       aes(x=factor(month), y=scaled)) + 
  geom_jitter(alpha=.1) +
  geom_hline(yintercept=1, lty=3, col='blue')

This dataset provides a scaling factor (column scaled) for each month, which is used to scale the season-wide density estimate to predict whale density in a given month. This dataset contains 1,000 bootstrap iterations (column iteration) of these estimates.

Close-encounter rates

To demonstrate the close-encounter rate estimation process, we will use the LNG Canada traffic:

data(lng_canada)

Whale dimensions and movement parameters are specified using the shipstrike functions fin_params() – which uses as defaults the values from Keen et al. for fin whales – or humpback_params – which provides humpback whale values as defaults. Either function can be used to define new parameters of your own.

fin_params()
## $length_min
## [1] 10
## 
## $length_max
## [1] 26
## 
## $length_mean
## [1] 20
## 
## $length_sd
## [1] 1.65
## 
## $width_factor
## [1] 0.2074
## 
## $speed_day_min
## [1] 0.27
## 
## $speed_day_max
## [1] 2.22
## 
## $speed_day_mean
## [1] 1.140301
## 
## $speed_day_sd
## [1] 0.4375298
## 
## $speed_night_min
## [1] 0.55
## 
## $speed_night_max
## [1] 3.04
## 
## $speed_night_mean
## [1] 1.809226
## 
## $speed_night_sd
## [1] 0.5718535
## 
## $delta_day_mean
## [1] 0.8087532
## 
## $delta_day_sd
## [1] 1.332026
## 
## $delta_night_mean
## [1] 1.494988
## 
## $delta_night_sd
## [1] 1.386166

The shipstrike function for creating simulations that estimates the close-encounter rate is encounter_rate(). To demonstrate this function, we conduct only 1 run of the simulator (100 iterations per run). The simulator will return an encounter rate estimate for each vessel-month-diel scenario. The diagnostic plot shows each iteration wtihin the run (black lines = whale tracks; blue line up center = ship track; large red dots = close-encounters).

mr <- encounter_rate(vessels = lng_canada,
                     whales = fin_params(),
                     outcome_dir = 'tests/',
                     month_batches = list(all = 1:12),
                     runs  = 1, 
                     iterations = 100, 
                     toplot = TRUE)
## 
## Type LNG Canada tanker in-heel :: Months of all :: diel period day :: simulating encounters ....
## 
## Encounter rates: 0.08
## 
## Type LNG Canada tanker in-heel :: Months of all :: diel period night :: simulating encounters ....

## 
## Encounter rates: 0.13
## 
## Type LNG Canada tanker in-product :: Months of all :: diel period day :: simulating encounters ....

## 
## Encounter rates: 0.1
## 
## Type LNG Canada tanker in-product :: Months of all :: diel period night :: simulating encounters ....

## 
## Encounter rates: 0.03
## 
## Type LNG Canada tug in-heel :: Months of all :: diel period day :: simulating encounters ....

## 
## Encounter rates: 0.02
## 
## Type LNG Canada tug in-heel :: Months of all :: diel period night :: simulating encounters ....

## 
## Encounter rates: 0.02
## 
## Type LNG Canada tug in-product :: Months of all :: diel period day :: simulating encounters ....

## 
## Encounter rates: 0.01
## 
## Type LNG Canada tug in-product :: Months of all :: diel period night :: simulating encounters ....

## 
## Encounter rates: 0.02
## 
## Finished!

To explore these results more closely, shipstrike will return a highly detailed data.frame when encounter_rate() is called with a single run (runs argument is 1). You can then pass that dataframe to the shipstrike function named encounter_diagnostics(), which will step through each iteration and display a detailed series of plots. Here is one example:

# Lots of details:
mr %>% nrow
## [1] 187913
# Filter to single vessel type and diel period
mri <- mr %>% filter(type=='LNG Canada tanker in-heel',
                         diel=='day')

# Try out diagnostics for last iteration in run
ggenc <- encounter_diagnostics(mri, id=100)
## finding 8 transit(s) with encounters: 17, 42, 56, 69, 71, 76, 81, 84
## 
## Transit 100 : no encounter.
## 
## Analyzing ellipses ...
## Creating timeline plot ...
## Creating map ...

Here is the code for running the full fin whale close-encounter rate analysis as done in Keen et al. (2023):

# Vessels 2019
data(ais_2019) 
encounter_rate(vessels = ais_2019, 
               whales = fin_params(),
               outcome_dir = 'tests/fw/ais_2019/',
               month_batches = list(winter = c(0:4, 11:12),
                                    summer = 5:10),
               runs  = 100, iterations = 100, toplot = FALSE)

# Vessels 2015
data(ais_2015)
encounter_rate(vessels = ais_2015,
               whales = fin_params(),
               outcome_dir = 'tests/fw/ais_2015/',
               month_batches = list(winter = c(0:4, 11:12),
                                    summer = 5:10),
               runs  = 100, iterations = 100, toplot = FALSE)

# LNG Canada
data(lng_canada) 
encounter_rate(vessels = lng_canada, 
               whales = fin_params(),
               outcome_dir = 'tests/fw/lng_canada/',
               month_batches = list(all = 1:12),
               runs  = 100, iterations = 100, toplot = FALSE)

# Cedar LNG
data(cedar_lng) 
encounter_rate(vessels = cedar_lng, 
               whales = fin_params(),
               outcome_dir = 'tests/fw/cedar_lng/',
               month_batches = list(all = 1:12),
               runs  = 100, iterations = 100, toplot = FALSE)

Note that the function is saving results to a file named p_encounter.RData within the directory specified by the argument outcome_dir. The result looks like this:

penc <- readRDS('../tests/fw/lng_canada/p_encounter.RData')
penc %>% head(10)
##                         type month diel  i p_encounter
## 1  LNG Canada tanker in-heel   all  day  1        0.11
## 2  LNG Canada tanker in-heel   all  day  2        0.09
## 3  LNG Canada tanker in-heel   all  day  3        0.08
## 4  LNG Canada tanker in-heel   all  day  4        0.09
## 5  LNG Canada tanker in-heel   all  day  5        0.06
## 6  LNG Canada tanker in-heel   all  day  6        0.08
## 7  LNG Canada tanker in-heel   all  day  7        0.07
## 8  LNG Canada tanker in-heel   all  day  8        0.10
## 9  LNG Canada tanker in-heel   all  day  9        0.11
## 10 LNG Canada tanker in-heel   all  day 10        0.05
# Plot it
ggplot(penc, 
       aes(x=p_encounter, fill=type, color=type)) + 
  geom_density(alpha=.5)

Ship-strike analysis

To predict whale-ship interaction outcomes, use the shipstrike function outcome_predict(). This function requires the following datasets:

  1. Marine traffic grid: prepared above.

  2. Whale density surface (either a single estimate or iterated bootstraps): shown above.

  3. Close-encounter rate distribution: prepared above.

  4. Whale depth distribution: we provide the fin whale depth distribution used in Keen et al. (2023) as a built-in dataset for shipstrike:

data(p_surface)
p_surface %>% head(10)
## # A tibble: 10 × 4
## # Groups:   diel [1]
##    diel      z  p_mean   p_sd
##    <chr> <dbl>   <dbl>  <dbl>
##  1 day     0   0       0     
##  2 day     0.5 0.00573 0.0152
##  3 day     1   0.0873  0.0480
##  4 day     1.5 0.109   0.0765
##  5 day     2   0.144   0.0717
##  6 day     2.5 0.147   0.0770
##  7 day     3   0.174   0.0697
##  8 day     3.5 0.192   0.0656
##  9 day     4   0.228   0.0500
## 10 day     4.5 0.230   0.0566
# Plot it
ggplot(p_surface, aes(y=z, x=p_mean, color=diel)) + 
  geom_point() + ylim(250, 0)

  1. Model parameters for the logistic relationship between P(Collision) and vessel speed: we provide the model parameters used in Keen et al. (2023) as a shipstrike dataset:
data(p_collision)
p_collision[4,]
##                                                             type         c1
## 4 Cargo > 180m | Passenger > 180m | Tanker > 180m | Other > 100m -0.2002002
##         c2 asymptote
## 4 11.80241       0.9

You can explore this model and/or develop your own version of it using the helper function, collision_curve(), which provides P(Collision) at every speed from 0 to 30 knots (this can be changed using the function arguments) according to the model parameters c1, c2, and asymptote.

collision_curve() %>% head
##        speed Pcollision
## 1 0.00000000 0.07744518
## 2 0.03003003 0.07787177
## 3 0.06006006 0.07830050
## 4 0.09009009 0.07873135
## 5 0.12012012 0.07916435
## 6 0.15015015 0.07959951
# Plot it
ggplot(collision_curve(), aes(x=speed, y=Pcollision)) + 
  geom_line(lwd=2, color='firebrick') + 
  ylim(0,1)

  1. Model parameters for the logistic relationship between P(Lethality) and vessel speed: we provide the model used in Keen et al. (2023) as a shipstrike dataset:
data(p_lethality)
p_lethality[4,]
##                                                             type        c1
## 4 Cargo > 180m | Passenger > 180m | Tanker > 180m | Other > 100m -1.241241
##          c2 asymptote
## 4 0.2712432         1

Similar to above, you can explore this model and/or develop your own version of it using the helper function, lethality_curve(), which provides P(Lethality) at every speed from 0 to 30 knots (this can be changed using the function arguments) according to the model parameters c1, c2, and asymptote.

lethality_curve() %>% head
##        speed Plethality
## 1 0.00000000  0.2242200
## 2 0.03003003  0.2256401
## 3 0.06006006  0.2270665
## 4 0.09009009  0.2284993
## 5 0.12012012  0.2299384
## 6 0.15015015  0.2313838
# Plot it
ggplot(lethality_curve(), aes(x=speed, y=Plethality)) + 
  geom_line(lwd=2, color='firebrick') + 
  ylim(0,1)

 

Below is the code that replicates the fin whale analysis in Keen et al. (2023):

# Common datasets to all analyses ==============================================
species <- 'fw'

load('tests/fw/dsm-bootstraps.RData')
whale <- bootstraps

load('tests/fw/seasonal_posterior.RData')
seasonal <- seasonal_boot


# 2019 traffic  ================================================================

data(ais_2019)
outcome_predict(traffic = ais_2019,
                scale_factors=NULL,
                whale = whale,
                seasonal = seasonal,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/fw/ais_2019/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species = 'fw',
                year = 2019,
                iterations = 1000)


# 2015 traffic  ================================================================

data(ais_2015)
outcome_predict(traffic = ais_2015,
                scale_factors=NULL,
                whale = whale,
                seasonal = seasonal,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/fw/ais_2015/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species = 'fw',
                year = 2015,
                iterations = 1000)

# 2030 traffic =================================================================

data(ais_2019) 
scale_factors <- readRDS('tests/ais/vessel_trends.RData')
outcome_predict(traffic = ais_2019,
                scale_factors=scale_factors,
                whale = whale,
                seasonal = seasonal,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/fw/ais_2030/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species = 'fw',
                year = 2030,
                iterations = 1000)


# LNG Canada (8 - 14 knots)  ===================================================

data(lng_canada) 
traffic <- lng_canada
vessels <- unique(traffic$type)

# Prepare avoidance / lethality models to recognize LNG traffic
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=' | '))

# Run outcome predictions
outcome_predict(traffic = lng_canada,
                scale_factors=NULL,
                whale = whale,
                seasonal = seasonal,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/fw/lng_canada/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species = 'fw',
                year = 2030,
                iterations = 1000)

# Cedar LNG (8 - 14 knots)  ====================================================

data(cedar_lng) 
traffic <- cedar_lng
vessels <- unique(traffic$type)

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=' | '))

# Run outcome predictions
outcome_predict(traffic = cedar_lng,
                scale_factors=NULL,
                whale = whale,
                seasonal = seasonal,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/fw/cedar_lng/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species = 'fw',
                year = 2030,
                iterations = 1000)

A few details of note:

  • Predictions will be saved in two files: outcomes.RData (all predictions for each iteration, summed across the spatial grid, and outcomes_grid.RData (the spatial grid of predictions) within the directory specified by the output_dir argument.

  • The argument asymptote_scaling can be used to manually scale the asymptotes of the collision and lethality models; for example, if the P(Collision) asymptote is 0.9 and asymptote_scaling is set to 0.5, the P(Collision) asymptote will be modified to 0.45.

  • The argument scale_factors can be used to increase/decrease the number of transits for each vessel type within traffic. This is how we use 2019 traffic data to predict 2030 AIS traffic. That scale_factors object looks like this:

scale_factors <- readRDS('../tests/ais/vessel_trends.RData')
scale_factors
## # A tibble: 10 × 10
##    type   km_2014 km_2015 km_2018 km_2019 km_2030 scale…¹ rate_…² rate_…³ pvalue
##    <chr>    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>
##  1 Cargo…      0       0      0     9444.  16943.    1.79  0.147   0.0820 0.300 
##  2 Fishi…  17165.  16377. 31949.   45784.  86437.    1.89  0.122   0.0646 0.0477
##  3 Other…  62201.  87269. 57399.   18723.      0     0    -0.482  -0.235  0.241 
##  4 Other…  25213.  19011. 27281.   26166.  33552.    1.28  0.0332  0.0259 0.441 
##  5 Other…  27391.   9943. 23408.   29184.  37725.    1.29  0.0497  0.0385 0.603 
##  6 Passe…   6212.   5489. 10403.    6168.  11554.    1.87  0.0693  0.0370 0.548 
##  7 Pleas…      0       0     55.8  51215.  91952.    1.80  0.147   0.0820 0.299 
##  8 Saili…    357.   1553. 12538.   20934.  50795.    2.43  0.191   0.0787 0.0242
##  9 Towin…  35311.  15299. 39409.   42838.  67175.    1.57  0.0755  0.0481 0.376 
## 10 Tug <…  30971.  44685. 46344.   43301.  61901.    1.43  0.0453  0.0317 0.335 
## # … with abbreviated variable names ¹​scale_factor, ²​rate_2019, ³​rate_2030

The only required fields are type and scale_factor.

Finally, we combine all of our fin whale results into a single list, to make them easier to find and use later:

fw <- list(params = fin_params(),
           p_encounter = list(ais_2015 = readRDS('tests/fw/ais_2015/p_encounter.RData'),
                              ais_2019 = readRDS('tests/fw/ais_2019/p_encounter.RData'),
                              lng_canada = readRDS('tests/fw/lng_canada/p_encounter.RData'),
                              cedar_lng = readRDS('tests/fw/cedar_lng/p_encounter.RData')),
           outcomes = list(ais_2015 = readRDS('tests/fw/ais_2015/outcomes.RData'),
                           ais_2019 = readRDS('tests/fw/ais_2019/outcomes.RData'),
                           ais_2030 = readRDS('tests/fw/ais_2030/outcomes.RData'),
                           lng_canada = readRDS('tests/fw/lng_canada/outcomes.RData'),
                           cedar_lng = readRDS('tests/fw/cedar_lng/outcomes.RData')),
           grid = list(ais_2015 = readRDS('tests/fw/ais_2015/outcomes_grid.RData'),
                       ais_2019 = readRDS('tests/fw/ais_2019/outcomes_grid.RData'),
                       ais_2030 = readRDS('tests/fw/ais_2030/outcomes_grid.RData'),
                       lng_canada = readRDS('tests/fw/lng_canada/outcomes_grid.RData'),
                       cedar_lng = readRDS('tests/fw/cedar_lng/outcomes_grid.RData')))

Results

The shipstrike package provides a number of functions to aid in summarizing, visualizing, and interpreting the results of outcome_predict(). To demonstrate these functions, we’ll focus on the predictions for all large ships in 2030 (AIS + LNG Canada + Cedar LNG):

large_ships <- c('Passenger > 180m', 'Cargo > 180m', 
                 'LNG Canada tanker in-heel',
                 'LNG Canada tanker in-product',
                 'Cedar LNG tanker in-heel',
                 'Cedar LNG tanker in-product')
# Spatially summarized outcomes
results <- rbind(fw$outcomes$ais_2030,
                 fw$outcomes$lng_canada,
                 fw$outcomes$cedar_lng) 
results %>% head
##   species year channel      vessel month diel iteration cooccurrence encounter
## 1      fw 2030     WRI Other > 40m     1  day         1            0         0
## 2      fw 2030     WRI Other > 40m     1  day         2            0         0
## 3      fw 2030     WRI Other > 40m     1  day         3            0         0
## 4      fw 2030     WRI Other > 40m     1  day         4            0         0
## 5      fw 2030     WRI Other > 40m     1  day         5            0         0
## 6      fw 2030     WRI Other > 40m     1  day         6            0         0
##   surface surface2 collision1.1 collision1.2 collision1.3 collision1.4
## 1       0        0            0            0            0            0
## 2       0        0            0            0            0            0
## 3       0        0            0            0            0            0
## 4       0        0            0            0            0            0
## 5       0        0            0            0            0            0
## 6       0        0            0            0            0            0
##   collision2.1 collision2.2 collision2.3 collision2.4 mortality1.1 mortality1.2
## 1            0            0            0            0            0            0
## 2            0            0            0            0            0            0
## 3            0            0            0            0            0            0
## 4            0            0            0            0            0            0
## 5            0            0            0            0            0            0
## 6            0            0            0            0            0            0
##   mortality1.3 mortality1.4 mortality2.1 mortality2.2 mortality2.3 mortality2.4
## 1            0            0            0            0            0            0
## 2            0            0            0            0            0            0
## 3            0            0            0            0            0            0
## 4            0            0            0            0            0            0
## 5            0            0            0            0            0            0
## 6            0            0            0            0            0            0
# Filter to only large ships > 180m
results <- results %>% filter(vessel %in% large_ships)
# Spatially explicit outcomes (summed across iterations)
results_grid <- rbind(fw$grid$ais_2030,
                      fw$grid$lng_canada,
                      fw$grid$cedar_lng)
results_grid %>% head
##   species year channel      vessel month diel grid_id cooccurrence encounter
## 1      fw 2030     WRI Other > 40m     1  day    4422            0         0
## 2      fw 2030     WRI Other > 40m     1  day    4423            0         0
## 3      fw 2030     WRI Other > 40m     1  day    4531            0         0
## 4      fw 2030     WRI Other > 40m     1  day    4639            0         0
## 5      fw 2030     WRI Other > 40m     1  day    4640            0         0
## 6      fw 2030     WRI Other > 40m     1  day    4748            0         0
##   surface surface2 collision1.1 collision1.2 collision1.3 collision1.4
## 1       0        0            0            0            0            0
## 2       0        0            0            0            0            0
## 3       0        0            0            0            0            0
## 4       0        0            0            0            0            0
## 5       0        0            0            0            0            0
## 6       0        0            0            0            0            0
##   collision2.1 collision2.2 collision2.3 collision2.4 mortality1.1 mortality1.2
## 1            0            0            0            0            0            0
## 2            0            0            0            0            0            0
## 3            0            0            0            0            0            0
## 4            0            0            0            0            0            0
## 5            0            0            0            0            0            0
## 6            0            0            0            0            0            0
##   mortality1.3 mortality1.4 mortality2.1 mortality2.2 mortality2.3 mortality2.4
## 1            0            0            0            0            0            0
## 2            0            0            0            0            0            0
## 3            0            0            0            0            0            0
## 4            0            0            0            0            0            0
## 5            0            0            0            0            0            0
## 6            0            0            0            0            0            0
# Filter to only large ships > 180m
results_grid <- results_grid %>% filter(vessel %in% large_ships)

outcome_table()

This function summarizes the predictions for each stage of whale-vessel interaction, from cooccurrence to mortality.

outcome_table(results)
## # A tibble: 20 × 6
##    event          mean median    q5   q95   q20
##    <fct>         <dbl>  <dbl> <dbl> <dbl> <dbl>
##  1 cooccurrence 166.      166   143   189   153
##  2 encounter     11.7      12     6    17     9
##  3 surface        5.07      5     2     9     3
##  4 surface2       5.1       5     2     9     3
##  5 collision1.1   2.29      2     0     5     1
##  6 collision1.2   2.63      2     0     6     1
##  7 collision1.3   2.57      2     0     5     1
##  8 collision1.4   5.07      5     2     9     3
##  9 collision2.1   2.27      2     0     5     1
## 10 collision2.2   2.68      3     0     6     1
## 11 collision2.3   2.65      3     0     5     1
## 12 collision2.4   5.1       5     2     9     3
## 13 mortality1.1   2.06      2     0     5     1
## 14 mortality1.2   2.41      2     0     5     1
## 15 mortality1.3   2.36      2     0     5     1
## 16 mortality1.4   4.54      4     1     8     3
## 17 mortality2.1   2.03      2     0     4     1
## 18 mortality2.2   2.44      2     0     5     1
## 19 mortality2.3   2.42      2     0     5     1
## 20 mortality2.4   4.54      4     1     8     3

The events with numbers can be read as follows: the first digit is the ‘strike zone’ scenario (1= 1x vessel draft, 2= 1.5x vessel draft); the second digit is the avoidance scenario (1 = P(Avoidance) is constant at 0.55; 2 = P(Avoidance) is a function of vessel speed; 3 = a deprecated duplicated of scenario 2; and 4 = No avoidance).

outcome_histograms()

This function produces pretty histograms of the posterior distributions of outcome predictions (the default is to only show the outcomes for scenarios in which strike zone is 1.5x draft and avoidance is a function of vessel speed).

outcome_histograms(results)$simple

outcome_map()

This function produces a map of a risk metric. The default is to show cooccurrences for all month-diel-vessels.

outcome_map(results_grid)

The default base map is of the Kitimat Fjord System (using the bangarang function gg_kfs()), but you can supply your own base map as an input.

outcome_shares()

This function uses a bootstrapping routine (described in Keen et al. (2023)) to estimate the share of risk attributable to each level for vessel type, channel, month, and diel period. It does this for each outcome type (cooccurrence, close encounter, …, mortality).

shares <- outcome_shares(results)
## Melting outcome data ...
## Calculating share of risks by vessel type ...
## Calculating share of risks by channel ...
## Calculating share of risks by month ...
## Calculating share of risks by diel period ...
shares$vessel 
## # A tibble: 120 × 3
## # Groups:   event [20]
##    event        vessel                        prop
##    <fct>        <chr>                        <dbl>
##  1 cooccurrence Cargo > 180m                    11
##  2 cooccurrence Cedar LNG tanker in-heel         3
##  3 cooccurrence Cedar LNG tanker in-product      3
##  4 cooccurrence LNG Canada tanker in-heel       21
##  5 cooccurrence LNG Canada tanker in-product    21
##  6 cooccurrence Passenger > 180m                42
##  7 encounter    Cargo > 180m                    10
##  8 encounter    Cedar LNG tanker in-heel         3
##  9 encounter    Cedar LNG tanker in-product      3
## 10 encounter    LNG Canada tanker in-heel       25
## # … with 110 more rows
shares$channel
## # A tibble: 160 × 3
## # Groups:   event [20]
##    event        channel  prop
##    <fct>        <chr>   <dbl>
##  1 cooccurrence CAA        30
##  2 cooccurrence CMP        13
##  3 cooccurrence EST         0
##  4 cooccurrence MCK         0
##  5 cooccurrence SQU        56
##  6 cooccurrence VER         0
##  7 cooccurrence WHA         0
##  8 cooccurrence WRI         1
##  9 encounter    CAA        26
## 10 encounter    CMP        10
## # … with 150 more rows
shares$month
## # A tibble: 240 × 3
## # Groups:   event [20]
##    event        month  prop
##    <fct>        <int> <dbl>
##  1 cooccurrence     1     0
##  2 cooccurrence     2     0
##  3 cooccurrence     3     0
##  4 cooccurrence     4     1
##  5 cooccurrence     5     6
##  6 cooccurrence     6    18
##  7 cooccurrence     7    22
##  8 cooccurrence     8    26
##  9 cooccurrence     9    19
## 10 cooccurrence    10     6
## # … with 230 more rows
shares$diel
## # A tibble: 40 × 3
## # Groups:   event [20]
##    event        diel   prop
##    <fct>        <chr> <dbl>
##  1 cooccurrence day      82
##  2 cooccurrence night    18
##  3 encounter    day      79
##  4 encounter    night    21
##  5 surface      day      70
##  6 surface      night    30
##  7 surface2     day      69
##  8 surface2     night    31
##  9 collision1.1 day      70
## 10 collision1.1 night    30
## # … with 30 more rows

outcome_chances()

This function predicts the chances of certain outcome severities.

outcome_chances(results)
## $at_least
##   Chances_of Collisions Mortalities
## 1       Zero        7.7         9.5
## 2 At least 1       92.3        90.5
## 3 At least 2       74.1        69.0
## 4 At least 3       51.0        44.3
## 5 At least 4       28.4        23.6
## 6 At least 5       14.2        11.0
## 
## $no_more_than
##   Chances_of Collisions Mortalities
## 1       Zero        7.7         9.5
## 2   Max of 1       25.9        31.0
## 3   Max of 2       49.0        55.7
## 4   Max of 3       71.6        76.4
## 5   Max of 4       85.8        89.0
## 6   Max of 5       94.2        96.2

outcome_melt()

This function may be handy if you want to conduct your own tailored analyses. It melts the results data into a tidyverse-friendly format in which every row is a single prediction. Every function above calls this function internally.

results_melted <- outcome_melt(results)

results_melted %>% head
##   species year       vessel channel month diel iteration        event outcome
## 1      fw 2030 Cargo > 180m     WRI     1  day         1 cooccurrence       0
## 2      fw 2030 Cargo > 180m     WRI     1  day         2 cooccurrence       0
## 3      fw 2030 Cargo > 180m     WRI     1  day         3 cooccurrence       0
## 4      fw 2030 Cargo > 180m     WRI     1  day         4 cooccurrence       0
## 5      fw 2030 Cargo > 180m     WRI     1  day         5 cooccurrence       0
## 6      fw 2030 Cargo > 180m     WRI     1  day         6 cooccurrence       0

outcome_validation()

This function replicates the validation exercise from Keen et al. (2023) to determine what our strike detection rate would need to be in order for our results to be plausible given that we have never observed a strike in our study area.

vali <- outcome_validation(results)
## Melting outcomes & prepping the posterior ...
## Determining the probability of your observations ...
## Likelihood of your observation, assuming perfect detection = 0
## Finding the strike detection rate (SDR) that would make your observations plausible ...
## preparing L ~ SDR plot ...
## --- SDR needed for P(Observation) of 0.05 = 0.14
## --- SDR needed for P(Observation) of 0.10 = 0.115
## --- SDR needed for P(Observation) of 0.20 = 0.08
## --- SDR needed for P(Observation) of 0.55 = 0.03
vali$L_plot

vali$SDR_plot

Mitigation

Keen et al. (2023) evaluates the efficacy of four main types of mitigation measures: speed reductions, daytime-only transits, seasonal displacement/rescheduling of ship transits, and seasonal moratoria on ship transits.

Speed reductions

To test for the efficacy of speed reductions, Keen et al. compares measures focused solely upon new LNG traffic to measures applied to all large ships > 180m. Changes in speed may effect both the close-encounter rate as well as the rate of collision/mortality.

LNG-only

First estimate the close-encounter rate with the new speed profile. This is possible using a convenience input, new_speeds, for the encounter_rate() function, which allows us to avoid re-simulating an entire new traffic scheme for LNG Canada.

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

We then pass these new encounter rates to the outcome_predict() function, making sure to also modify the speeds that will be used in predicting collision/lethality:

load('tests/fw/dsm-bootstraps.RData')
load('tests/fw/seasonal_posterior.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=' | ')))

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

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

That was for LNG Canada. Now repeat for Cedar LNG:

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

# Outcomes
load('tests/fw/dsm-bootstraps.RData')
load('tests/fw/seasonal_posterior.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 = seasonal_boot,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/fw/mitigation/1a/cedar_lng/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

All large-ship traffic

To assess the effect of speed restrictions on all large ships, we re-run the 2019 AIS results, starting with encounter_rate(), using the arguments speed_restriction and lengths_restricted to specify that this restriction only applies to ships over 180m.

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

# Predict outcomes on AIS 2030 with < 9 kn
load('tests/fw/dsm-bootstraps.RData')
load('tests/fw/seasonal_posterior.RData')
data(p_surface)
data(ais_2019)
scale_factors <- readRDS('tests/ais/vessel_trends.RData')
data(p_collision)
data(p_lethality)
#
# # Modify AIS data
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 = scale_factors,
                whale = bootstraps,
                seasonal = seasonal_boot,
                p_encounter_dir = 'tests/fw/mitigation/1b/ais_2030/',
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/fw/mitigation/1b/ais_2030/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species='fw',
                year=2030,
                iterations = 1000)

Daytime-only transits

For this mitigation measure, we do not need to re-run the encounter rate – just outcome predictions.

LNG only

First for LNG Canada …

load('tests/fw/dsm-bootstraps.RData')
load('tests/fw/seasonal_posterior.RData')
data(p_surface)
data(lng_canada)
(vessels <- unique(lng_canada$type))

# Coerce all traffic to daytime
lng_canada_mod <- lng_canada 
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 = seasonal_boot,
                p_encounter_dir = 'tests/fw/lng_canada/',
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/fw/mitigation/2a/lng_canada/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

… then for Cedar LNG:

load('tests/fw/dsm-bootstraps.RData')
load('tests/fw/seasonal_posterior.RData')
data(p_surface)
data(cedar_lng)
(vessels <- unique(cedar_lng$type))

# Coerce all traffic to daytime
cedar_lng_mod <- cedar_lng 
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 = seasonal_boot,
                p_encounter_dir = 'tests/fw/cedar_lng/',
                surface = p_surface,
                avoidance = avoidance,
                lethality = lethality,
                outcome_dir = 'tests/fw/mitigation/2a/cedar_lng/',
                asymptote_scaling = NULL,
                month_batches = list(all = 0:12),
                iterations = 1000)

All large-ship traffic

load('tests/fw/dsm-bootstraps.RData')
load('tests/fw/seasonal_posterior.RData')
scale_factors <- readRDS('tests/ais/vessel_trends.RData')
data(p_collision)
data(p_lethality)
data(p_surface)
data(ais_2019)

# Coerce all traffic to daytime
ais_mod <- ais_2019 
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/fw/ais_2019/',
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/fw/mitigation/2b/ais_2030/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species = 'fw',
                year = 2030,
                iterations = 1000)

The mitigated outcomes for speed reductions and daytime transits have the exact same structure as the unmitigated ship-strike predictions, and all of the same outcome_ results functions can be used to explore them.

Seasonal displacement of LNG traffic

To test the effects of seasonally displacing LNG traffic (i.e., rescheduling transits from a given month into other months), we use the shipstrike function mitigate_loop(), which applies the displacement window to each candidate month in a loop. This allows us to determine which month would be the most efficacious target for mitigation.

We specify that reschedule is TRUE to make it clear that we are displacing traffic during the mitigation window, not canceling it outright.

# Baseline scenario (no mitigation)
results <- readRDS('tests/fw/results.RData')
outcomes <- rbind(results$outcomes$lng_canda,
                  results$outcomes$cedar_lng)

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


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

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

The mitigate_loop function returns a list of results. For example:

# Result of mitigate loop
mr <- readRDS('../tests/fw/mitigation/3c/mitigation.RData')
mr %>% names

The $outcomes slot is the result of outcome_table for each mitigation window (field test):

mr$outcomes

The $chances slot is the result of outcome_chances()$at_least for each mitigation window:

mr$chances %>% head(10)

And the $hist slot provides all iterations of the mortality2.2 outcome for each mitigation window:

mr$hist

ggplot(mr$hist, aes(x=test, y=outcome)) + 
  geom_jitter(alpha = .3)

Seasonal moratorium on LNG traffic

To assess the effects of a sesonal moratorium, we use the same mitigate_loop() function, this time setting the reschedule argument to FALSE.

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

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

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

Now compile all of these results into a single list, for ease of access and use later on:

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

Humpback whales

The humpback whale analysis is a near-replication of the fin whale analysis, with only a few minor changes. The chief difference is that the humpback whale density surface model produced a different density surface for each month. This means that a seasonal scaling curve was not needed, in contrast to fin whales above. Instead, we have a bootstrapped density surface for each month.

Close-encounter rates

# Common inputs
(whales <- humpback_params())
runs <- 100
iterations <- 100

# Vessels 2019
data(ais_2019) ; vessels <- ais_2019
outcome_dir <- 'tests/hw/ais_2019/'
encounter_rate(vessels, whales, outcome_dir,
               month_batches = list(winter = c(0:4, 11:12),
                                    summer = 5:10),
               runs  = runs, iterations = iterations, toplot = FALSE)

# Vessels 2015
data(ais_2015)
encounter_rate(vessels = ais_2015,
               whales = whales,
               outcome_dir = 'tests/hw/ais_2015/',
               month_batches = list(winter = c(0:4, 11:12),
                                    summer = 5:10),
               runs  = 100, iterations = 100, toplot = FALSE)

# LNG Canada
data(lng_canada) ; vessels <- lng_canada
outcome_dir <- 'tests/hw/lng_canada/'
encounter_rate(vessels, whales, outcome_dir,
               month_batches = list(all = 1:12),
               runs  = runs, iterations = iterations, toplot = TRUE)

# Cedar LNG
data(cedar_lng) ; vessels <- cedar_lng
outcome_dir <- 'tests/hw/cedar_lng/'
encounter_rate(vessels, whales, outcome_dir,
               month_batches = list(all = 1:12),
               runs  = runs, iterations = iterations, toplot = FALSE)

Ship-strike analysis

# Common datasets to all analyses ==============================================

species <- 'hw'
whale <- readRDS('tests/hw/pwhale_seasonal_boots.RData')
seasonal <- NULL

data(p_surface)
surface <- p_surface

asymptote_scaling <- NULL

# 2019 traffic  ================================================================

data(ais_2019) ; traffic <- ais_2019
(vessels <- unique(traffic$type))

p_encounter_dir <- NULL
(outcome_dir <- 'tests/hw/ais_2019/') %>% dir

data(p_collision) ; (avoidance <- p_collision)
data(p_lethality) ; (lethality <- p_lethality)

month_batches <- list(winter = c(0:4, 11:12), summer = 5:10)

outcome_predict(traffic,
                scale_factors=NULL,
                 whale,
                 seasonal,
                 p_encounter_dir,
                 surface,
                 avoidance,
                 lethality,
                 outcome_dir,
                 asymptote_scaling = asymptote_scaling,
                 month_batches = month_batches,
                 species=species,
                 year=2019,
                 iterations = 1000)


# 2015 traffic  ================================================================

data(ais_2015)
data(p_collision)
data(p_lethality)
outcome_predict(traffic = ais_2015,
                scale_factors=NULL,
                whale = whale,
                seasonal = NULL,
                p_encounter_dir = NULL,
                surface = p_surface,
                avoidance = p_collision,
                lethality = p_lethality,
                outcome_dir = 'tests/hw/ais_2015/',
                asymptote_scaling = NULL,
                month_batches = list(winter = c(0:4, 11:12), summer = 5:10),
                species = species,
                year = 2015,
                iterations = 1000)

# 2030 traffic =================================================================

data(ais_2019) ; traffic <- ais_2019
(vessels <- unique(traffic$type))
(scale_factors <- readRDS('tests/ais/vessel_trends.RData'))

p_encounter_dir <- 'tests/hw/ais_2019/'
outcome_dir <- 'tests/hw/ais_2030/'

data(p_collision) ; (avoidance <- p_collision)
data(p_lethality) ; (lethality <- p_lethality)

month_batches <- list(winter = c(0:4, 11:12), summer = 5:10)

# Run outcome predictions
outcome_predict(traffic,
                scale_factors,
                 whale,
                 seasonal,
                 p_encounter_dir,
                 surface,
                 avoidance,
                 lethality,
                 outcome_dir,
                 asymptote_scaling = asymptote_scaling,
                 month_batches = month_batches,
                 species=species,
                 year=2030,
                 iterations = 1000)


# LNG Canada (8 - 14 knots)  ===================================================

data(lng_canada) ; traffic <- lng_canada
(vessels <- unique(traffic$type))

p_encounter_dir <- NULL
outcome_dir <- 'tests/hw/lng_canada/'

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=' | ')))
avoidance

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

month_batches <- list(all = 0:12)

# Run outcome predictions
outcome_predict(traffic,
                scale_factors=NULL,
                 whale,
                 seasonal,
                 p_encounter_dir,
                 surface,
                 avoidance,
                 lethality,
                 outcome_dir,
                 asymptote_scaling = asymptote_scaling,
                 month_batches = month_batches,
                 species=species,
                 year=2030,
                 iterations = 1000)


# Cedar LNG (8 - 14 knots)  ====================================================

data(cedar_lng) ; traffic <- cedar_lng
(vessels <- unique(traffic$type))

p_encounter_dir <- NULL
outcome_dir <- 'tests/hw/cedar_lng/'

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=' | ')))
avoidance

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

month_batches <- list(all = 0:12)

# Run outcome predictions
outcome_predict(traffic,
                scale_factors=NULL,
                whale,
                seasonal,
                p_encounter_dir,
                surface,
                avoidance,
                lethality,
                outcome_dir,
                asymptote_scaling = asymptote_scaling,
                month_batches = month_batches,
                species=species,
                year=2030,
                iterations = 1000)

Now compile all of those outcomes into a list:

hw <- list(params = humpback_params(),
           p_encounter = list(ais_2015 = readRDS('tests/hw/ais_2015/p_encounter.RData'),
                              ais_2019 = readRDS('tests/hw/ais_2019/p_encounter.RData'),
                              lng_canada = readRDS('tests/hw/lng_canada/p_encounter.RData'),
                              cedar_lng = readRDS('tests/hw/cedar_lng/p_encounter.RData')),
           outcomes = list(ais_2015 = readRDS('tests/hw/ais_2015/outcomes.RData'),
                           ais_2019 = readRDS('tests/hw/ais_2019/outcomes.RData'),
                           ais_2030 = readRDS('tests/hw/ais_2030/outcomes.RData'),
                           lng_canada = readRDS('tests/hw/lng_canada/outcomes.RData'),
                           cedar_lng = readRDS('tests/hw/cedar_lng/outcomes.RData')),
           grid = list(ais_2015 = readRDS('tests/hw/ais_2015/outcomes_grid.RData'),
                       ais_2019 = readRDS('tests/hw/ais_2019/outcomes_grid.RData'),
                       ais_2030 = readRDS('tests/hw/ais_2030/outcomes_grid.RData'),
                       lng_canada = readRDS('tests/hw/lng_canada/outcomes_grid.RData'),
                       cedar_lng = readRDS('tests/hw/cedar_lng/outcomes_grid.RData')))

Mitigation

The mitigation analysis is a mirror of that used for fin whales above:

Speed reductions

LNG-only

First for 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)

… then for 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)

All large ships

# 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)

Daytime-only transits

LNG-only

First for 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)

… then for 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)

All large ships

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)

Seasonal displacement of LNG traffic

# 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')

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')
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')))

Refer to the Supplemental Material for Keen et al. (2023) for all code that produces the tables and figures within that publication.