R Notebook

noaa <- readLines('http://www.ssd.noaa.gov/PS/TROP/2018/adt/archive.html')
library(stringr)
library(readr)
pattern = 'http[^\"]+text[^\"]+\\.txt'
urls <- str_extract(noaa,pattern)
urls <- urls[!is.na(urls)]
header1 <- c("date", "time", "intensity_ci", "intensity_mslp", "intensity_vmax","tno_fnl","tno_adj","tno_ini","constrnt","wkng","rpd","temp_ctr","temp_mean","scene","est_rmw","mw","lat","lon","mthd","sat","vza","comments")
widths1 <- c(9,7,5,7,6,5,4,4,10,5,5,8,7,8,7,6,8,8,6,9,5,NA)
header2 = header <- c("date", "time", "intensity_ci", "intensity_mslp", "intensity_mslpplat", "intensity_vmax","tno_fnl","tno_adj","tno_ini","constrnt","wkng","rpd","temp_ctr","temp_mean","scene","est_rmw","mw","lat","lon","mthd","sat","vza","comments")
widths2 <- c(9,7,5,7,8,6,5,4,4,10,5,5,8,7,8,7,6,8,8,6,9,5,NA)
first <- T
q <- 1
#REMOVE
for(url in urls) {
  print(q)
  q <- q + 1
  lines <- readLines(url)
  name <- strsplit(substr(url,52,10000), "-")[[1]][1]
  # check if "BiasAdj" in line 4 -- if so, need an extra column!
  if(length(grep("BiasAdj",lines))) {
    widths <- widths2
    h <- header2
  } else {
    widths <- widths1
    h <- header1
  }
  
  if(url=="http://www.ssd.noaa.gov/PS/TROP/DATA/2018/adt/text/17P-list.txt") {
    # Special rule for this file due to incorrect negative MSLP values offsetting columns
    for(i in 5:55) {
      lines[i] = paste0(substr(lines[i],1,22), substr(lines[i],24,nchar(lines[i])), sep='')
    }
  }
  lines <- lines[5:(length(lines)-3)]
  
  data <- read_fwf(paste0(lines,collapse="\n"), fwf_widths(widths, col_names=h), na="N/A")
  data$est_rmw <- as.character(data$est_rmw)
  data$typhoon_name <- name
  if(first) {
    first = F
    typhoon <- data
  } else {
    typhoon <- bind_rows(typhoon, data)
  }
}
write.csv(typhoon,"../../Data/shipping/typhoon.csv", row.names=F)

library(plotly)
library(dplyr)
# load data
df <- read.csv("../../Data/Shipping/shipping_dataset.csv", stringsAsFactors = F, na="NA")
df_ships <- read.csv("../../Data/Shipping/index.csv", stringsAsFactors = F, na="NA")
df_ports <- read.csv("../../Data/Shipping/port_dataset.csv", stringsAsFactors = F, na="NA")
typhoon <- read.csv("../../Data/Shipping/typhoon.csv", stringsAsFactors = F)
# rename lat and lon to avoid conflicts
df_ports <- df_ports %>% rename(port_lat=lat, port_lon=lon)
# Join in ships and ports
df <- left_join(df, df_ships)

Joining, by = "imo"

df <- df %>% filter(!imo == 9528029 | reported_destination == "AU ADL")
df <- left_join(df,df_ports, by=c("left_port"="port"))
df <- df %>% rename(left_lat=port_lat, left_lon=port_lon)
df <- left_join(df,df_ports, by=c("arrive_port"="port"))
df <- df %>% rename(arrive_lat=port_lat, arrive_lon=port_lon)
# fix dates
df$last_update = as.POSIXct(df$last_update, tz="UTC",origin="1970-01-01")
df$left_time = as.POSIXct(df$left_time, tz="UTC",origin="1970-01-01")
df$arrive_time = as.POSIXct(df$arrive_time, tz="UTC",origin="1970-01-01")
# Fix typhoon dates
typhoon$date <- as.POSIXct(paste(typhoon$date,typhoon$time), format="%Y%b%d %H%M%S", tz="UTC", origin="1970-01-01")
# Fix typhoon longitude -- they have lon * -1 in the data
typhoon$lon <- typhoon$lon * -1
typhoon_all <- typhoon
# filter to dates in sample
typhoon <- typhoon %>% filter(date > as.POSIXct("2018-08-30", tz="UTC"))
df <- df %>% filter(last_update > as.POSIXct("2018-08-30", tz="UTC"))
geo <- list(
  showland = TRUE,
  showlakes = TRUE,
  showcountries = TRUE,
  showocean = TRUE,
  countrywidth = 0.5,
  landcolor = toRGB("grey90"),
  lakecolor = toRGB("aliceblue"),
  oceancolor = toRGB("aliceblue"),
  projection = list(
    type = 'orthographic',  # detailed at https://plot.ly/r/reference/#layout-geo-projection
    rotation = list(
      lon = 100,
      lat = 1,
      roll = 0
    )
  ),
  lonaxis = list(
    showgrid = TRUE,
    gridcolor = toRGB("gray40"),
    gridwidth = 0.5
  ),
  lataxis = list(
    showgrid = TRUE,
    gridcolor = toRGB("gray40"),
    gridwidth = 0.5
  )
)
library(RColorBrewer)
# plot with boats, ports, and typhoons
palette = brewer.pal(8, "Dark2")[c(1,8,3,2)]
p <- plot_geo(colors=palette) %>%
  add_markers(data=df_ports, x = ~port_lon, y = ~port_lat, color = "Port") %>%
  add_markers(data=df[df$run == 1,], x = ~lon, y = ~lat, color = ~ship_type,
              text=~paste('Ship name',shipname)) %>%
  add_markers(data=typhoon[typhoon$date > as.POSIXct("2018-08-31 00:00:00", tz="UTC") & typhoon$date < as.POSIXct("2018-09-01 00:00:00", tz="UTC"),], x = ~lon, y = ~lat, color="TYPHOON", colors=toRGB("red"), text=~paste("Name", typhoon_name)) %>%
   layout(
    showlegend = TRUE, geo = geo,
    title = 'Singaporean owned container and tanker ships, August 31, 2018'
  )
p

# plot with boats, ports, and typhoons
palette = brewer.pal(8, "Dark2")[c(1,3,2)]
p <- plot_geo(colors=palette) %>%
  add_markers(data=df[df$run == 14,], x = ~lon, y = ~lat, color = ~ship_type,
              text=~paste('Ship name',shipname)) %>%
  add_markers(data=typhoon[typhoon$typhoon_name == "JEBI_Y",], x = ~lon, y = ~lat, color="Typhoon Jebi", text=~paste("Name", typhoon_name, "</br>Time: ", date)) %>%
   layout(
    showlegend = TRUE, geo = geo,
    title = 'Singaporean container/tanker ships, September 4, 2018, evening'
  )
p

world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))

p <- ggplot(data = world1) + 
  geom_sf() + 
  geom_point(data = df, aes(x = lon, y = lat, frame=run))
ggplotly(p) %>%
   animation_opts(
    1500, easing = "elastic", redraw = FALSE
  )

# distance between two geo coordinates
library(geosphere)

# Linear interpolation for matching hurricane locations
linear_interp <- function(from, to, mid, d_from, d_to) {
  d_from * (mid - from) / (to - from) + d_to * (to - mid) / (to - from)
}

# for line in df
# line_last_update, line$lat, line_lon
# typhoon %>% group_by(typhoon_name) %>%

df$date_only <- as.Date(df$last_update, "UTC")
typhoon$ty_date_only <- as.Date(typhoon$date, "UTC")
typhoon <- typhoon %>% rename(ty_lat=lat, ty_lon=lon, ty_date=date)

first <- TRUE
for(d in unique(df$date_only)) {
  df2 <- df[df$date_only == d,]
  ty2 <- typhoon[typhoon$ty_date_only == d,]
  dfm <- as.matrix(df2[,c("lon","lat")])
  tym <- as.matrix(ty2[,c("ty_lon","ty_lat")])
  mat <- distm(dfm, tym, fun = distVincentyEllipsoid)
  df2 <- cbind(df2, ty2[max.col(-mat),])
  if(first) {
    df3 <- df2
    first <- FALSE
  } else {
    df3 <- bind_rows(df3, df2)
  }
}
write.csv(df3,"../../Data/shipping/combined.csv", row.names=F)

summary(fit)

Call:
glm(formula = delayed ~ (typhoon_100 + typhoon_500 + typhoon_1000):typhoon_intensity, 
    family = binomial, data = df3)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5770  -0.4648  -0.4648  -0.4648   2.2974  

Coefficients: (1 not defined because of singularities)
                                         Estimate Std. Error  z value Pr(>|z|)    
(Intercept)                               -2.1709     0.0144 -150.788   <2e-16 ***
typhoon_100:typhoon_intensity(0,34]       -0.2838     0.6195   -0.458   0.6469    
typhoon_100:typhoon_intensity(34,53]      -0.2609     0.4545   -0.574   0.5660    
typhoon_100:typhoon_intensity(53,92.4]    -7.7742    72.4631   -0.107   0.9146    
typhoon_100:typhoon_intensity(92.4,150]        NA         NA       NA       NA    
typhoon_500:typhoon_intensity(0,34]       -0.1955     0.1779   -1.099   0.2717    
typhoon_500:typhoon_intensity(34,53]       0.3250     0.1949    1.668   0.0953 .  
typhoon_500:typhoon_intensity(53,92.4]     0.2382     0.2517    0.946   0.3440    
typhoon_500:typhoon_intensity(92.4,150]   -0.5626     0.7576   -0.743   0.4577    
typhoon_1000:typhoon_intensity(0,34]       0.1935     0.1001    1.932   0.0533 .  
typhoon_1000:typhoon_intensity(34,53]      0.1373     0.1468    0.935   0.3496    
typhoon_1000:typhoon_intensity(53,92.4]    0.1409     0.1610    0.875   0.3814    
typhoon_1000:typhoon_intensity(92.4,150]   0.1686     0.1889    0.892   0.3722    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 37178  on 55846  degrees of freedom
Residual deviance: 37157  on 55835  degrees of freedom
  (7204 observations deleted due to missingness)
AIC: 37181

Number of Fisher Scoring iterations: 8

OLD BELOW

library(readr)
library(dplyr)
library(ggplot2)

url_csv <- 'https://raw.githubusercontent.com/d4tagirl/R-Ladies-growth-maps/master/rladies.csv'
rladies <- read_csv(url(url_csv)) %>% 
  select(-1)
rladies2 <- rladies
rladies$frame <- 1
rladies2$frame <- 2
rladies2$lat <- rladies2$lat + 5
rladies3 <- rbind(rladies,rladies2)
rladies2$frame <- 3
rladies2$lat <- rladies2$lat + 5
rladies3 <- rbind(rladies3,rladies2)
rladies2$frame <- 4
rladies2$lat <- rladies2$lat + 5
rladies3 <- rbind(rladies3,rladies2)
rladies2$frame <- 5
rladies2$lat <- rladies2$lat + 5
rladies3 <- rbind(rladies3,rladies2)

library(plotly)

g <- list(resolution = 50,
          showland = TRUE,
          showlakes = TRUE,
          landcolor = toRGB("grey80"),
          countrycolor = toRGB("grey80"),
          lakecolor = toRGB("white"),
          projection = list(type = "equirectangular"),
          coastlinewidth = 2,
          lataxis = list(
            range = c(20, 60),
            showgrid = TRUE,
            tickmode = "linear",
            dtick = 10
          ),
          lonaxis = list(
            range = c(-100, 20),
            showgrid = TRUE,
            tickmode = "linear",
            dtick = 20
          ))

plot_geo() %>%
  add_markers(data = rladies3,
              x = ~lon,
              y = ~lat,
              frame = ~frame,
              color=I("purple"),
              hoverinfo = "text",
              text = ~paste('city: ', location,
                            '<br /> created: ', created_at,
                            '<br /> followers: ', followers)) %>%
  layout(geo = g) %>%
  animation_opts(
    1000, easing = "elastic", redraw = FALSE
  )

maps, maptools, sf, ggplot @v3, rgeos

world1 <- sf::st_as_sf(map('world', plot = FALSE, fill = TRUE))

p <- ggplot(data = world1) + 
  geom_sf() + 
  geom_point(data = rladies3, aes(x = lon, y = lat, frame=frame))
ggplotly(p) %>%
   animation_opts(
    1500, easing = "elastic", redraw = FALSE
  )
# ANIMATION: https://github.com/plotly/plotly.js/blob/master/src/plots/animation_attributes.js