R/hurricanes.R

In cortinah/hockeystick: Download and Visualize Essential Climate Change Data

#' Download and plot essential climate data
#'
#' Retrieves Atlantic basin hurricane data since 1851 from National Oceanic and Atmospheric Administration
#' HURDAT Atlantic Hurricane Database Re-analysis Project.
#' \url{https://www.aoml.noaa.gov/hrd/hurdat/Data_Storm.html}
#'
#' @name get_hurricanes
#' @param use_cache (boolean) Return cached data if available, defaults to TRUE. Use FALSE to fetch updated data.
#' @param write_cache (boolean) Write data to cache, defaults to FALSE. Use TRUE to write data to cache for later use. Can also be set using options(hs_write_cache=TRUE)
#'
#' @return Invisibly returns a tibble with the annual HURDAT hurricane data since 1851
#'
#' @details `get_hurricanes` invisibly returns a tibble with NOAA's annual North Atlantic revised HURDAT hurricane data since 1851.
#' The returned object includes Year, and number of named storms, hurricanes, major hurricanes (S-S scale 3-5), Accumulated Cyclone Energy (ACE), and U.S. hurricane strikes.
#'
#' ACE is an index that combines the number of systems, how long they existed and how intense they became. It is calculated by squaring the maximum sustained surface wind in the system every six hours that the cyclone is a Named Storm and summing it up for the season. Please refer to above website for details.
#'
#' @import rvest
#'
#' @examples
#' \donttest{
#' # Fetch from cache if available:
#' hurricanes <- get_hurricanes()
#' #
#' # Force cache refresh:
#' hurricanes <- get_hurricanes(use_cache=FALSE)
#' #
#' # Review cache contents and last update dates:
#' hockeystick_cache_details()
#' #
#' # Plot output using package's built-in ggplot2 settings
#' plot_hurricanes(hurricanes) }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#' @references
#' HURDAT North Atlantic Hurricane Database Re-analysis Project, Hurricane Research Division, NOAA \url{https://www.aoml.noaa.gov/hrd/hurdat/Data_Storm.html}.
#'
#' Data from: \url{https://www.aoml.noaa.gov/hrd/hurdat/comparison_table.html}
#'
#'
#' \url{https://en.wikipedia.org/wiki/Accumulated_cyclone_energy}
#' @export
get_hurricanes <- function(use_cache = TRUE, write_cache = getOption("hs_write_cache")) {

hs_path <- tools::R_user_dir("hockeystick","cache")

if (use_cache) {
  if (file.exists(file.path(hs_path,'hurricanes.rds'))) return(invisible(readRDS((file.path(hs_path,'hurricanes.rds')))))
  }

file_url <- "https://www.aoml.noaa.gov/hrd/hurdat/comparison_table.html"
connected <- .isConnected(file_url)
if (!connected) {message("Retrieving remote data requires internet connectivity."); return(invisible(NULL))}

hurricanes <- read_html(file_url)
hurricanes <- html_node(hurricanes, xpath='//*[(@id = "tdcontent")]//table')
hurricanes <- html_table(hurricanes)

colnames(hurricanes) <- hurricanes[2,]
hurricanes <- hurricanes[3:nrow(hurricanes), c(1,3,5,7,9,11)]

hurricanes <- tidyr::separate(hurricanes, 6, sep=1, into=c("RevisedUShurricanes", NA), fill = "right")

hurricanes <- dplyr::mutate_all(hurricanes, as.numeric)

dir.create(hs_path, showWarnings = FALSE, recursive = TRUE)
if (write_cache) saveRDS(hurricanes, file.path(hs_path, 'hurricanes.rds'))

invisible(hurricanes)
}


#' Download and plot essential climate data
#'
#' Plots the hurricane data retrieved using `get_hurricanes()` with ggplot2. The output ggplot2 object may be further modified.
#'
#'
#' @name plot_hurricanes
#' @param dataset Name of the tibble generated by `get_hurricanes`
#' @param print (boolean) Display hurricane ggplot2 chart, defaults to TRUE. Use FALSE to not display chart.
#' @param cat (string) Select which category of hurricane to plot. May be "major", "hurricane", or "storm". Defaults to "major".
#'
#' @return Invisibly returns a ggplot2 object with hurricanes chart
#'
#' @details `plot_hurricanes` invisibly returns a ggplot2 object with a pre-defined hurricane data chart using data from `get_hurricanes`.
#' By default the chart is also displayed. Smooths using ggplot2's built-in loess smoother. Users may further modify the output ggplot2 chart.
#' Categories may be "major" (category 3-5), "hurricane" (category 1-5), or "storm" (named storm).
#'
#' @import ggplot2
#'
#' @examples
#' \donttest{
#' # Fetch hurricane data:
#' hurricanes <- get_hurricanes()
#' #
#' # Plot output using package's built-in ggplot2 defaults
#' plot_hurricanes(hurricanes)
#'
#' # Or just call plot_hurricanes(), which defaults to get_hurricanes() dataset
#' plot_hurricanes()
#'
#' p <- plot_hurricanes(hurricanes, print = FALSE)
#' # Modify plot such as: p + ggplot2::labs(title='Growing number of North Atlantic named storms') }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#'
#' @export
plot_hurricanes <- function(dataset = get_hurricanes(), cat='major', print=TRUE) {

  if(cat=='major') {dataset <- dataset[,c('Year', 'RevisedMajorHurricanes')]
                    sub_label <- 'Major Hurricanes (category 3-5)' }

  if(cat=='hurricane') {dataset <- dataset[,c('Year', 'RevisedHurricanes')]
                        sub_label <- 'Hurricanes (category 1-5)' }

  if(cat=='storm') {dataset <- dataset[,c('Year', 'RevisedNamedStorms')]
                    sub_label <- 'Named Storms' }

  colnames(dataset)[2] <- 'hurr_type'
  if (is.null(dataset) | ncol(dataset)!=2) return(invisible(NULL))

    plot <- ggplot(dataset, aes(x=Year, y=hurr_type)) +geom_line(alpha=0.75, aes(color='Hurricanes')) + theme_bw(base_size=12) +
    scale_x_continuous(name=NULL, breaks=seq(1850, 2025,25)) +
    scale_y_continuous(n.breaks = 8) +geom_smooth(size=1.1, se=F, aes(color='Loess smooth')) +
    scale_color_manual(name=NULL, values=c('dodgerblue2','firebrick1')) +theme(legend.position = c(0.15, 0.825),legend.background=element_blank()) +
    labs(title='North Atlantic Hurricane Basin', subtitle = sub_label,
        y='Number per season', caption='Source: NOAA Hurricane Research Division, HURDAT2 re-analysis\nhttps://www.aoml.noaa.gov/hrd/hurdat/Data_Storm.html')

  if (print) suppressMessages( print(plot) )
  invisible(plot)
}


#' Download and plot essential climate data
#'
#' Plots the hurricane energy data (ACE) retrieved using `get_hurricanes()` with ggplot2. The output ggplot2 object may be further modified.
#'
#'
#' @name plot_hurricane_nrg
#' @param dataset Name of the tibble generated by `get_hurricanes`
#' @param print (boolean) Display hurricane ggplot2 chart, defaults to TRUE. Use FALSE to not display chart.
#'
#' @return Invisibly returns a ggplot2 object with hurricane energy chart
#'
#' @details `plot_hurricane_nrg` invisibly returns a ggplot2 object with a pre-defined hurricane energy data chart using data from `get_hurricanes`.
#' By default the chart is also displayed. Smooths using ggplot2's built-in loess smoother. Users may further modify the output ggplot2 chart.
#'
#' ACE is an index that combines the number of systems, how long they existed and how intense they became. It is calculated by squaring the maximum sustained surface wind in the system every six hours that the cyclone is a Named Storm and summing it up for the season. Please refer to above website for details.
#'
#' @import ggplot2
#'
#' @examples
#' \donttest{
#' # Fetch hurricane data:
#' hurricanes <- get_hurricanes()
#' #
#' # Plot output using package's built-in ggplot2 defaults
#' plot_hurricane_nrg(hurricanes)
#'
#' # Or just call plot_hurricane_nrg(), which defaults to get_hurricanes() dataset
#' plot_hurricane_nrg()
#'
#' p <- plot_hurricane_nrg(hurricanes, print = FALSE)
#' # Modify plot such as: p + ggplot2::labs(title='Accumulated Cyclone Energy') }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#'
#' @export
plot_hurricane_nrg <- function(dataset = get_hurricanes(), print=TRUE) {

  if (is.null(dataset)) return(invisible(NULL))

  plot <- ggplot(dataset, aes(x=Year, y=RevisedACE)) +geom_line(alpha=0.75, aes(color='ACE')) + theme_bw(base_size=12) +
    scale_x_continuous(name=NULL, breaks=seq(1850, 2025,25)) +
    scale_y_continuous(n.breaks = 8) +geom_smooth(size=1.1, se=F, aes(color='Loess smooth')) +
    scale_color_manual(name=NULL, values=c('dodgerblue2','firebrick1')) +theme(legend.position = c(0.15, 0.875),legend.background=element_blank()) +
    labs(title='North Atlantic Hurricane Basin', subtitle='Accumulated Cyclone Energy (ACE)',
         y=expression('10'^4*' kn'^2 ), caption='Source: NOAA Hurricane Research Division, HURDAT2 re-analysis\nhttps://www.aoml.noaa.gov/hrd/hurdat/Data_Storm.html')

  if (print) suppressMessages( print(plot) )
  invisible(plot)
}