Project goal: The purpose of this assignment is to draw on your knowledge of the grid and ggplot2 package to build a new geom. You will then need to load and tidy the provided dataset in order to plot the geom on a map.

The submission of this assignment will include the following 2 files:

  1. An R script containing the code implementing the geom
  2. A plot file (in PNG or PDF format) containing a map that has the wind radii overlayed for Hurricane Ike

Build a custom geom for ggplot2 that can be used to add the hurricane wind radii chart for a single storm observation to a map (i.e., could be used to recreate the figure shown above). Use the geom to map the create a map showing the wind radii chart at one observation times for Hurricane Ike, which occurred in September 2008. Use an observation time when the storm was near or over the United States. Download the data on all storms in the Atlantic basin from 1988–2015 (extended best tracks), then tidy the dataset into “long” format.

Clean data to spec: (1) add a column for storm_id that combines storm name and year (the same storm name can be used in different years, so this will allow for the unique identification of a storm); (2) format the longitude to ensure it is numeric and has negative values for locations in the Western hemisphere (this will make it easier to use the longitude for mapping); (3) format and combine columns describing the date and time to create a single variable with the date-time of each observation; and (4) convert the data to a “long” format, with separate rows for each of the three wind speeds for wind radii (34 kts, 50 kts, and 64 kts), where kts is Knots, or nautical miles per hour

packages <- c('readr', 'stringr', 'tidyr', 'ggplot2', 'ggmap', 'geosphere', 'tidyverse', 'dplyr')

# Load packages
lapply(packages, require, character.only = TRUE)
## Loading required package: readr
## Loading required package: stringr
## Loading required package: tidyr
## Loading required package: ggplot2
## Loading required package: ggmap
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
## Loading required package: geosphere
## Loading required package: tidyverse
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ tibble  3.1.2     ✓ dplyr   1.0.6
## ✓ purrr   0.3.4     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## [[1]]
## [1] TRUE
## 
## [[2]]
## [1] TRUE
## 
## [[3]]
## [1] TRUE
## 
## [[4]]
## [1] TRUE
## 
## [[5]]
## [1] TRUE
## 
## [[6]]
## [1] TRUE
## 
## [[7]]
## [1] TRUE
## 
## [[8]]
## [1] TRUE
library(readr)
library(stringr)
library(tidyr)
library(ggplot2)
library(geosphere)
library(tidyverse)
library(dplyr)
citation("ggmap")
## 
## To cite ggmap in publications, please use:
## 
##   D. Kahle and H. Wickham. ggmap: Spatial Visualization with ggplot2.
##   The R Journal, 5(1), 144-161. URL
##   http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf
## 
## A BibTeX entry for LaTeX users is
## 
##   @Article{,
##     author = {David Kahle and Hadley Wickham},
##     title = {ggmap: Spatial Visualization with ggplot2},
##     journal = {The R Journal},
##     year = {2013},
##     volume = {5},
##     number = {1},
##     pages = {144--161},
##     url = {https://journal.r-project.org/archive/2013-1/kahle-wickham.pdf},
##   }
D. Kahle and H. Wickham. ggmap: Spatial Visualization with ggplot2.
The R Journal, 5(1), 144-161. URL
http://journal.r-project.org/archive/2013-1/kahle-wickham.pdf

`

————————————-functions———————————

# function name:      tidy_hurdata(data)
# this function does: subsets data to a single object
# input:       raw data
# output:      tidied data in long format
# usage:       tidy_hurdata <- function(data)

# parameter:   ext_tracks - table of fixed width formatted data in a data.frame

# R functions: stringr str_c str_to_title dplyr mutate_ select_ tidyr gather spread
# note:        none

tidy_hurdata <- function(data) {
  data %>% 
    # Config storm_id and date 
    dplyr::mutate_(storm_id = ~stringr::str_c(stringr::str_to_title(storm_name), year, sep = '-'),
                   date = ~stringr::str_c(year, '-', month, '-', day, ' ', hour, ':', '00', ':', '00'),
                   longitude = ~-longitude
    ) %>% 
    # Select relevant columns
    dplyr::select_(.dots = c('storm_id', 'date', 'longitude', 'latitude', 
                            'radius_34_ne', 'radius_34_se', 'radius_34_sw', 'radius_34_nw',
                            'radius_50_ne', 'radius_50_se', 'radius_50_sw', 'radius_50_nw',
                            'radius_64_ne', 'radius_64_se', 'radius_64_sw', 'radius_64_nw')
    ) %>%
    
    #convert wide to long format
    tidyr::gather(variable, value, -storm_id, -date,-latitude, -longitude, -storm_id, -date) %>% dplyr::mutate_(wind_speed = ~str_extract(variable, "(34|50|64)"),
              variable = ~str_extract(variable, "(ne|nw|se|sw)")) %>% tidyr::spread(variable, value) %>% select_(.dots = c('storm_id', 'date', 'latitude', 'longitude', 'wind_speed', 'ne', 'nw', 'se', 'sw'))
  
}
# function name:      filter_hurdata(data)
# this function does: filters data by hurricane name, datentime
# input:     tidied data
# output:    filtered by name and time
# usage:     filter_hurdata <- function(data, hurricane, observation) {
#             data <- filter_(data, ~storm_id == hurricane & date == observation)}

# parameters: data:         tidied data, 'cleandata' rows, cols to spec
#             hurricane:    Ike-2008 selected from data frame
#             observation:  hurricane & date like 'Ike-2008 '2008-09-13 12:00:00''
# R functions:     dplyr filter_
# note:            none

filter_hurdata <- function(data, hurricane, observation) {
  data <- filter_(data, ~storm_id == hurricane & date == observation)
  
}

———————- GEOM ———————————————-

# function name:      geom_hurricane()
# this function does: create new geom with radii using Google's get_map() & API key
# input:     data, lon, lat, ne, se, nw, sw, fill, color
# output:    geom plotted in ggplot2, get_map map rendered as background
# usage:     geom_hurricane(data = stormdf) +
#       scale_color_manual(name = "Wind speed (kts)",
#                          values = c("red",
#                                     "orange",
#                                     "yellow")) +
#       scale_fill_manual(name = "Wind speed (kts)",
#                         values = c("red",
#                               "orange",
#                               "yellow"))

# parameters:  data                storm data tidied to spec as dataframe
#              lon, lat            initial x,y point definition juxtaposed on map layer
#              ne, se, nw, sw      quadrant positions used in various calculations
#              fill, color:        color coded wind speeds
# R functions: cbind, destPoint(), rbind, bind_rows, ggplot, geom_polygon, as_tibble, mutate, ggmap, devtools, devtools::install_github

# note: The new geom uses the center of the hurricane as the center of the plot, with 4 separate quadrants. Each quadrant has its' own wind speed, given per quadrant with 3 categories to show the area covered per wind speed. From the geosphere package the function 'destPoint()' is used which calculates a "circle" around a coordinate (x,y), a start point. This point is the initial bearing (direction), and distance. The function computes the destination point traveling along the shortest path on an ellipsoid.

# note: I'm new at this, but have found that Google maps API's on Google Cloud will overload and crash.  Three keys are here, as alternatives, although using your own is advisable.  Daily cloud hit quotas exist. Best to be frugal to keep the API from the brink. It can be quite a wait. 

geom_hurricane <- function(data = data,
                      x = longitude,
                      y = latitude,
                      r_ne = ne,
                      r_se = se,
                      r_nw = nw,
                      r_sw = sw,
                      fill = wind_speed,
                      color = wind_speed)
       {
# make list
dframe <- as_tibble()

# create center o hurricane - hcenter
hcenter <- cbind(data$longitude,
                data$latitude)

# Store speed in data frame then make radius scale offset
rad_offset <- cbind(data$ne, # 0 deg at 12:00 -always in this sequence
                    data$se, 
                    data$sw, 
                    data$nw) 

# create the 4 quadrants
for (i in 1:4)
       {
       # we have 34, 50 and 64 kts areas
       for (j in 1:nrow(data))
              {
              # Generate the points
              destPoint(hcenter[j,],
                        b=((i-1)*90):(90*i),           # full circle -360 deg
                        d=rad_offset[j,i] * 1852) %>%  # radius-scale offset
                     
                     rbind(hcenter) %>% # Add hcenter & origins
                     
                            as_tibble() %>% # Convert regular dframe to tibble
                     
                                   mutate(i = i,      # Add columns i and j
                                          j = j) %>%  # use to filter
                                          
                                          bind_rows(dframe) -> dframe
              }
       }

# Load ggmap package
if(!requireNamespace("devtools")) install.packages("devtools")
devtools::install_github("dkahle/ggmap", ref = "tidyup")
library(ggmap)
library(devtools)

# register API Key recently generated on Google cloud -referenced to 3 API's
#register_google(key="AIzaSyAREeaI08wldmSYqZewuf3PtOjQGSdByKA")
register_google(key="AIzaSyC3SXSFPlS_Q0ZAx_z2x0ZXjui9Af7XoqY")
#register_google(key= "AIzaSyB0fKSElDN-a0LpvhvvWlFNP5CWCFf3jZM")

# Google Maps/Stratmen
get_map("Louisiana",
        zoom = 6,
        maptype = "toner-background") %>%
       
       # ggmap plot (citation above)
       ggmap(extent = "device") + 
              # all 4 quadrants at 64 knots
              geom_polygon(data = (dframe %>%        # NE quadrant, 64 kts
                                  filter(i %in% 1,   
                                         j %in% 1)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[1]),
                           alpha = 0.5) +
              
              geom_polygon(data = (dframe %>%        # SE quadrant, 64 kts
                                  filter(i %in% 2,   
                                         j %in% 1)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[1]),
                           alpha = 0.5) + 
              
              geom_polygon(data = (dframe %>%        # SW quadrant, 64 kts
                                  filter(i %in% 3,   
                                         j %in% 1)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[1]),
                           alpha = 0.5) +
              
              geom_polygon(data = (dframe %>%        # NW quadrant, 64 kts
                                  filter(i %in% 4,   
                                         j %in% 1)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[1]),
                           alpha = 0.5)+ 
              # all 4 quadrants at 50 knots
              geom_polygon(data = (dframe %>%        # NE quadrant, 50 kts
                                  filter(i %in% 1,    
                                         j %in% 2)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[2]),
                           alpha = 0.5) + 
              
              geom_polygon(data = (dframe %>%        # SE quadrant, 50 kts
                                  filter(i %in% 2,   
                                         j %in% 2)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[2]),
                           alpha = 0.5) + 
              
              geom_polygon(data = (dframe %>%        # SW quadrant, 50 kts
                                  filter(i %in% 3,   
                                         j %in% 2)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[2]),
                           alpha = 0.5) + 
              
              geom_polygon(data = (dframe %>%        # NW quadrant, 50 kts
                                  filter(i %in% 4,   
                                         j %in% 2)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[2]),
                           alpha = 0.5) + 
              # all 4 quadrants at 34 knots
              geom_polygon(data = (dframe %>%        # NE quadrant, 34 kts
                                  filter(i %in% 1,   
                                         j %in% 3)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[3]),
                           alpha = 0.5) + 
              
              geom_polygon(data = (dframe %>%        # SE quadrant, 34 kts
                                  filter(i %in% 2,   
                                         j %in% 3)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[3]),
                           alpha = 0.5) + 
              
              geom_polygon(data = (dframe %>%        # SW quadrant, 34 kts
                                  filter(i %in% 3,   
                                         j %in% 3)),
                           aes(x = lon, y = lat,
                           fill = data$wind_speed[3]),
                           alpha = 0.5) + 
              
              geom_polygon(data = (dframe %>%        # NW quadrant, 34 kts
                                  filter(i %in% 4,   
                                         j %in% 3)),
                           aes(x = lon, y = lat,
                               fill = data$wind_speed[3]),
                           alpha = 0.5)
}

———————–begin execution——————————————

# Load data 
# input:  unzipped to ./data/ebtrk_atlc_1988_2015.txt
# output: create ext_tracks as input to tidy section
# R functions: readr, read.fwf and fwf_widths
# note:        formatting of widths and columns designated by Roger Peng

ext_tracks_file <- paste0("http://rammb.cira.colostate.edu/research/",
                         "tropical_cyclones/tc_extended_best_track_dataset/",
                          "data/ebtrk_atlc_1988_2015.txt")

# vectors as specified
ext_tracks_widths <- c(7, 10, 2, 2, 3, 5, 5, 6, 4, 5, 4, 4, 5, 3, 4, 3, 3, 3,
                       4, 3, 3, 3, 4, 3, 3, 3, 2, 6, 1)
ext_tracks_colnames <- c("storm_id", "storm_name", "month", "day",
                          "hour", "year", "latitude", "longitude",
                          "max_wind", "min_pressure", "rad_max_wind",
                          "eye_diameter", "pressure_1", "pressure_2",
                          paste("radius_34", c("ne", "se", "sw", "nw"), sep = "_"),
                          paste("radius_50", c("ne", "se", "sw", "nw"), sep = "_"),
                          paste("radius_64", c("ne", "se", "sw", "nw"), sep = "_"),
                          "storm_type", "distance_to_land", "final")

# read_fwf (Read a table of fixed width formatted data into a data.frame)
ext_tracks <- read_fwf(ext_tracks_file, 
                       fwf_widths(ext_tracks_widths, ext_tracks_colnames),
                       na = "-99")
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double(),
##   storm_id = col_character(),
##   storm_name = col_character(),
##   month = col_character(),
##   day = col_character(),
##   hour = col_character(),
##   storm_type = col_character(),
##   final = col_character()
## )
## ℹ Use `spec()` for the full column specifications.
# Tidy
#Subset to a single observation time for a single storm
#input:  ext_tracks (output from load)
#output: 3x9 tibble - cleandata[1:3, 1:9]
ext_tracks <- as.data.frame(ext_tracks)
cleandata<- tidy_hurdata(ext_tracks)
## Warning: `select_()` was deprecated in dplyr 0.7.0.
## Please use `select()` instead.
## Warning: `mutate_()` was deprecated in dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
#Filter
#Data now shows a single storm, a single date, a single center 
#storm location, three rows for each wind intensity (34, 50, and 64 kts), and cols for
#the wind radii for a given intensity in each of the four quadrants
#input:  data, hurricane, observation
#output: long format dataframe

storm_observation <- filter_hurdata(cleandata, hurricane = 'Ike-2008', observation = '2008-09-13 12:00:00')
## Warning: `filter_()` was deprecated in dplyr 0.7.0.
## Please use `filter()` instead.
## See vignette('programming') for more help
stormdf <- as.data.frame(storm_observation)
stormdf
##   storm_id                date latitude longitude wind_speed  ne nw  se  sw
## 1 Ike-2008 2008-09-13 12:00:00     30.3     -95.2         34 125 60 180 125
## 2 Ike-2008 2008-09-13 12:00:00     30.3     -95.2         50  75 45  90  60
## 3 Ike-2008 2008-09-13 12:00:00     30.3     -95.2         64  50 20  45  30

Ike-2008 data now ready as dataframe, ‘stormdf’ in long format Next, render geom with google map layer in backgound; windplot layer in foreground

geom_hurricane(data = stormdf) +
       
       scale_color_manual(name = "Wind speed (kts)",
                          values = c("red",
                                     "orange",
                                     "yellow")) +
       
       scale_fill_manual(name = "Wind speed (kts)",
                         values = c("red",
                               "orange",
                               "yellow"))
## Warning: The `x` argument of `as_tibble()` can't be missing as of tibble 3.0.0.
## Loading required namespace: devtools
## Skipping install of 'ggmap' from a github remote, the SHA1 (2d756e5e) has not changed since last install.
##   Use `force = TRUE` to force installation
## Loading required package: usethis
## maptype = "toner-background" is only available with source = "stamen".
## resetting to source = "stamen"...
## Source : https://maps.googleapis.com/maps/api/staticmap?center=Louisiana&zoom=6&size=640x640&scale=2&maptype=terrain&key=xxx
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Louisiana&key=xxx
## Source : http://tile.stamen.com/toner-background/6/14/24.png
## Source : http://tile.stamen.com/toner-background/6/15/24.png
## Source : http://tile.stamen.com/toner-background/6/16/24.png
## Source : http://tile.stamen.com/toner-background/6/14/25.png
## Source : http://tile.stamen.com/toner-background/6/15/25.png
## Source : http://tile.stamen.com/toner-background/6/16/25.png
## Source : http://tile.stamen.com/toner-background/6/14/26.png
## Source : http://tile.stamen.com/toner-background/6/15/26.png
## Source : http://tile.stamen.com/toner-background/6/16/26.png
## Source : http://tile.stamen.com/toner-background/6/14/27.png
## Source : http://tile.stamen.com/toner-background/6/15/27.png
## Source : http://tile.stamen.com/toner-background/6/16/27.png