In BenioffOceanInitiative/shipr: ships for whales functions

knitr::opts_chunk$set(echo = TRUE, comment = NA)

library(knitr)
library(kableExtra)
library(tidyverse)
library(leaflet)
library(RColorBrewer)
library(maptools)
library(sf)
library(rgeos)
library(dplyr)
library(RPostgreSQL)
library(RPostgres)
library(dbplyr)
library(lubridate)
library(units)
library(data.table)
library(parallel)
library(ggplot2)
library(here)

{.tabset}

Source R files or Github Installation

Source or Install

# Database connection function
source('~/github/whalesafe4r/R/db.R')

# Crawl https://ais.sbarc.org/logs_delimited/ for URLs
source('~/github/whalesafe4r/R/crawlers.R')

# Read ais.txt URL into ais_data df and read ais_data df into ais_segs df
source('~/github/whalesafe4r/R/readers.R')

# 1) Loop through ais.txt URLs to return ais_data
# 2) Create ais_segments from ais_data
# 3) Intersect ais_segments and vsr_zones table in the database to create vsr_segments table in the database
source('~/github/whalesafe4r/R/update_ais.R')

# Utility functions for crawler, date and spatial functions, etc...
source('~/github/whalesafe4r/R/utils.R')

# Create summary statistics for merged vsr_sements/ihs_data df
source('~/github/whalesafe4r/R/seg_stats.R')

#devtools::install_github("BenioffOceanInitiative/whalesafe4r")

Initiate Connection & Get AIS URLs

Initiate Database Connection

# Must be on UCSB IP address & database credentials for connection to work. 
# (add to repo s4w_amazon_rds.yml file)
# (add yml file to your gitignore)

con <- db_connect()

Get last read URL from ais_data table

#Query database to find the 'last read' url from the ais_data table
last_read_url <- dbGetQuery(con, "SELECT url AS last_read_url FROM ais_data WHERE datetime = (SELECT MAX(datetime) FROM ais_data) LIMIT 1;") %>% .$last_read_url

last_read_url

Get List of New Links

# Use last_read_url in the get_ais_urls() function to get all new ais urls

new_links <- get_ais_urls(last_read_url)

new_links

Example Links

# Example get_ais_urls() from November 10, 2019, 00 hour

example_path <- "https://ais.sbarc.org/logs_delimited/2019/191110/AIS_SBARC_191110-00.txt"

example_links <- get_ais_urls(example_path)

example_links %>% head() 
example_links%>% tail()

Get New AIS and Segments Data

New AIS Data

# Use update_ais_data() to loop through 'new_links' and return new_ais_data data.frame

new_ais_data <- get_ais_data(links = new_links)

example_ais_data <- get_ais_data(links = example_links)

example_ais_data_5000 <- example_ais_data %>% 
  head(5000)

Append ais_data Table in Database with new_ais_data

# Append new_ais_data to ais_data table in the database
# Appending to 'test_ais_data' for testing
update_ais_data(ais_data = example_ais_data_5000)

Get new_segs_data for global environment

# Use new_ais_data in the update_segments_data() function to return new_segs_data 
new_segs_data <- get_segment_data(ais_data = new_ais_data)

# Use example_ais_data in the update_segments_data() function to return example_segs_data 
example_segs_data <- get_segment_data(ais_data = example_ais_data_5000)

Write ais_segs_data to ais_segments Table in Database

# write new_segs_data to the database if new_segs_data exists
# Writing to 'test_ais_segments' for rmd
update_segments_data(segs_data = new_segs_data)

update_segments_data(segs_data = example_segs_data)

Update vsr_segments Table in the Database

Update vsr_segments

# Update vsr_segments on the database side
update_vsr_segments(segs_data = new_segs_data)

temp_vsr_segments <- update_vsr_segments(segs_data = example_segs_data)

Map VSR segments sample data

# con <- db_connect()
vsr_segments_ex <- sf::st_read(dsn = con, EWKB = TRUE,
                               query = "SELECT * FROM vsr_segments
                               limit 1000;")
# Get 2019 VSR Zone
vsr_zones <- sf::st_read(dsn = con, EWKB = TRUE, 
                               query = "SELECT * FROM vsr_zones;")

vsr_2019 <- head(vsr_zones,n = 1L)

# VSR Segments map example
ex_vsr_seg_map <- leaflet(vsr_segments_ex) %>% 
  setView(lng = '-119', lat = '34',zoom = 7) %>% 
  addTiles() %>% 
  addPolylines(label = ~sprintf("%0.03f km/hr on %s", seg_kmhr, beg_dt, name), group = "segments") %>% 
  addPolygons(data = vsr_2019, fillColor = 'blue', group = 'vsr') %>%   
  addLayersControl(overlayGroups = c("segments","vsr"))

ex_vsr_seg_map

Calculate Statistics

Merge vsr_segments data with ihs_data

# Merge vsr_segments with ihs_data from database to return vsr_ihs_data in global environment. (Takes a little time ~60 seconds)
system.time({
vsr_ihs_data <- get_vsr_ihs_data()
})

Calculate Ship Statistics

# Calculate ship stats for each mmsi using merged df, vsr_ihs_data
ship_stats <- ship_statistics(data = vsr_ihs_data)

ship_stats_2019 <- ship_statistics(data = vsr_ihs_data,
                                   date_start = '2019-01-01',
                                   date_end = '2019-12-31',
                                   tonnage = 300)

Overwrite the ship_stats Table in the Database

# And overwrite the ship_stats table in the database 
dbWriteTable(con = con, 
             name = "ship_stats", 
             value = ship_stats, 
             overwrite = TRUE)

Calculate Operator Statistics

# Calculate operator stats for each mmsi 
operator_stats = operator_statistics(data = vsr_ihs_data)

operator_stats_2019 = operator_statistics(data = vsr_ihs_data,
                                     date_start = '2019-01-01',
                                     date_end = '2019-12-31',
                                     tonnage = 300) %>% filter("total distance (km)" >= 100)

Overwrite the operator_stats Table in the Database

# And overwrite the operator_stats table in the database 
dbWriteTable(con = con, 
             name = "operator_stats", 
             value = operator_stats, 
             overwrite = TRUE)

Disconnect from database

Disconnect from Database

dbDisconnect(con)

Maps & Tables

example_ais_data_table <- kable(head(example_ais_data, n = 100L)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% 
  scroll_box(width = "100%", height = "300px")

Example AIS Points data table

example_ais_data_table

ex_ais_data_map <- leaflet(head(example_ais_data, n=1000L)) %>% 
        addTiles() %>% 
        addCircles(lat = example_ais_data$lat, lng = example_ais_data$lon) %>% 
        setView(lng = '-119', lat = '34',zoom = 7)

Example AIS Points Map

ex_ais_data_map

ex_seg_map <- leaflet(example_segs_data) %>% 
                    setView(lng = '-119', lat = '34',zoom = 7) %>% 
                    addTiles() %>% 
                    addPolylines()

Example Segments Map

#different df that what's shown in ex_ais_data_map
ex_seg_map

Ship & Operator Statistics Tables

ship_stats_2019_table <- ship_stats_2019 %>% 
  kable(escape = F) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "responsive")) %>% 
  scroll_box(width = "100%", height = "300px")

Ship Statistics Table

ship_stats_2019_table

operator_stats_2019_table <- operator_stats_2019 %>% 
  kable(escape = F) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "responsive")) %>% 
  scroll_box(width = "100%", height = "300px")

Operator Statistics Table

operator_stats_2019_table

mmsi & Operator Counts

mmsi Count (IHS Data)

con <- db_connect()
ihs_data <- dbGetQuery(con, "SELECT * FROM ihs_data;")

length(unique(ihs_data$mmsi))

mmsi Count (IHS Data) >= 300 gross tonnage

length(unique(ihs_data$mmsi[ihs_data$gt>=300]))

mmsi Count (VSR Segments Data)

# unique mmsi's in vsr_segments database table
dbGetQuery(con, "select count(distinct mmsi) from vsr_segments;")

2019 mmsi Count (VSR Segments Data)

# 897 unique mmsi's in vsr_segments database table in 2019
dbGetQuery(con, "select count(distinct mmsi) from vsr_segments where datetime >= '2019-01-01';")

dbDisconnect(con)

mmsi Count (VSR Segments & IHS Data)

# 1441 unique mmsi's (merged df) in vsr in 2018 & 2019
length(unique(vsr_ihs_data$mmsi))

2019 mmsi Count (VSR Segments & IHS Data)

# unique mmsi's in vsr in 2019
length(unique(vsr_ihs_data$mmsi[vsr_ihs_data$datetime>='2019-01-01']))

2019 mmsi Count (VSR Segments & IHS Data) >= 300 gross tonnage

# unique mmsi's in vsr in 2019 with gross tonnage >= 300
length(unique(vsr_ihs_data$mmsi[vsr_ihs_data$datetime>='2019-01-01' & vsr_ihs_data$gt>=300]))

DOUBLE CHECKIN'

# unique mmsi's in vsr in 2019
length(unique(ship_stats_2019$mmsi))

# unique mmsi's in vsr in 2019 with gross tonnage >= 300
length(unique(ship_stats_2019$mmsi[ship_stats_2019$gt>=300]))

2019 mmsi Count (VSR Segments & IHS Data) >= 300 gross tonnage & >= 100 km travelled

#  unique mmsi's in vsr in 2019 with gross tonnage >= 300 and travelled >= 100 km
length(unique(ship_stats_2019$mmsi[ship_stats_2019$gt>=300 & ship_stats_2019$`total distance (km)`>=100]))

Operator Counts

# 1095 operators total in IHS data
length(unique(ihs_data$operator))

Operators with vessels >= 300 gross tonnage

# 1062 operators over 300 gt
length(unique(ihs_data$operator[ihs_data$gt>=300]))

Operators in Merged VSR segments and IHS data

# 467 unique operators in vsr in 2018 & 2019
length(unique(vsr_ihs_data$operator))

Operators in Merged VSR segments and IHS data (2019)

# 267 unique operators in vsr in 2019
length(unique(vsr_ihs_data$operator[vsr_ihs_data$year==2019]))

Operators in Merged VSR segments and IHS data (2019 & gt >= 300)

# 257 unique operators in vsr in 2019 with gross tonnage >= 300
length(unique(vsr_ihs_data$operator[vsr_ihs_data$year==2019 & vsr_ihs_data$gt>=300]))

Double Checking

# DOUBLE CHECKIN'
# 257 unique operators in vsr in 2019 with gt>=300
length(unique(operator_stats_2019$operator))

# 257 unique operators in vsr in 2019 with vessels gt >= 300
length(unique(ship_stats_2019$operator[ship_stats_2019$gt>=300]))

Operators in Merged VSR segments and IHS data (2019, gt >= 300, distance >= 100 km)

# 208 operators with over 100 km travelled total
length(unique(operator_stats_2019$operator[operator_stats_2019$`total distance (km)`>=100]))

Grade Graphs

Operator Grades (basic rubric)

operator_stats_2019_100k=operator_stats_2019 %>% 
  filter(operator_stats_2019$`total distance (km)`>=100)

theme_set(theme_bw()) 

op_grade <- ggplot(data.frame(operator_stats_2019_100k), aes(x=grade)) +
  geom_bar() +  
  theme(plot.title = element_text(hjust = 0.5)) + 
  ggtitle("Operator VSR Cooperation Grade Distribution: 2019") + 
  labs(caption="Source: Benioff Ocean Initiative",
       x="VSR Cooperation Grade",
       y="Count") +
  geom_text(stat='count', aes(label=..count..), vjust=1.5, color="white")+
      # geom_text(stat='count', aes(label=..count..), position = position_stack(vjust = 0.5),size=4, color = "white") + 
    theme(plot.title = element_text(size = 14))

op_grade

Vessel Grades (modified NOAA Grades)

ship_stats_2019_100k=ship_stats_2019 %>% 
  filter(ship_stats_2019$`total distance (km)`>=100)

ship_grades <- ggplot(data.frame(ship_stats_2019_100k), aes(x=noaa_grade)) +
  geom_bar() +  
  theme(plot.title = element_text(hjust = 0.5)) + 
    theme(plot.subtitle = element_text(hjust = 0.5)) +
  ggtitle("Ship VSR Cooperation Grade Distribution: 2019", subtitle = "Modified NOAA Rubric") + 
  labs(caption="Source: Benioff Ocean Initiative",
       x="VSR Cooperation Grade",
       y="Count") +
    geom_text(stat='count', aes(label=..count..), vjust=1.5, color="white")

ship_grades

summary_stats = data.table(  
  nm_0_10=sum(ship_stats_2019$`distance (nautcal miles) 0-10 knots`),
    nm_10_12=sum(ship_stats_2019$`distance (nautcal miles) 10-12 knots`),
    nm_12_15=sum(ship_stats_2019$`distance (nautcal miles) 12-15 knots`),
    nm_over_15=sum(ship_stats_2019$`distance (nautcal miles) over 15 knots`))

summary_stats_percents = data.table(  
  `% Travelled 0-10 kn` = (sum(ship_stats_2019_100k$`distance (nautcal miles) 0-10 knots`)/sum(ship_stats_2019_100k$`total distance (nautcal miles)`)*100),
  `% Travelled 10-12 kn`= sum(ship_stats_2019_100k$`distance (nautcal miles) 10-12 knots`)/sum(ship_stats_2019_100k$`total distance (nautcal miles)`)*100,
  `% Travelled 12-15 kn` = sum(ship_stats_2019_100k$`distance (nautcal miles) 12-15 knots`)/sum(ship_stats_2019_100k$`total distance (nautcal miles)`)*100,
  `% Travelled over 15 kn` = sum(ship_stats_2019_100k$`distance (nautcal miles) over 15 knots`)/sum(ship_stats_2019_100k$`total distance (nautcal miles)`)*100)

sum_stats=melt(summary_stats_percents)

Speed Distribution Pie Chart

pie <- ggplot(sum_stats, aes(x="", y=value, fill=variable)) + geom_bar(stat="identity", width=1) + 
  coord_polar("y", start=0) +
   geom_text(aes(label = paste0(round(value), "%")), position = position_stack(vjust = 0.5)) +
  scale_fill_manual(values=c("palegreen1", "skyblue1", "gold1" ,"coral2")) + 
  labs(x = NULL, y = NULL, fill = NULL, title = "Summary of VSR cooperation for 2019") + 
  theme_classic() + theme(axis.line = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          plot.title = element_text(hjust = 0.5, color = "#666666")) +  theme(plot.title = element_text(hjust = -0.5)) + theme(plot.title = element_text(size = 16))

pie

BenioffOceanInitiative/shipr documentation built on Aug. 10, 2020, 1:39 a.m.