R/icecurves.R

In hockeystick: Download and Visualize Essential Climate Change Data

  #' Download and plot essential climate data
#'
#' Retrieves Arctic or Antarctic monthly Sea Ice Index time series (in million square km).
#' Source is the National Snow and Ice Data Center, defaults to Arctic (Northern Hemisphere) monthly sea ice extent since 1979.
#' \url{https://nsidc.org/data/explore-data} \cr \cr
#'
#' @name get_icecurves
#' @param pole 'N' for Arctic or 'S' for Antarctic
#' @param measure Must be 'extent' or 'area', defaults to 'extent'. Please see terminology link in references for details.
#' @param use_cache (boolean) Return cached data if available, defaults to TRUE. Use FALSE to fetch updated data, or to change pole or month in cache.
#' @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 series of monthly Sea Ice Index since 1979 (in million square km).
#'
#' `get_icecurves` invisibly returns a tibble with time series of monthly Sea Ice Index since 1979 (in million square km).
#'
#' User may select monthly Arctic or Antarctic sea ice extent or area (in millions of square kilometers).
#' Defaults to Arctic sea ice extent.  \url{https://nsidc.org/sea-ice-today/about-data#area_extent}
#'
#' @importFrom utils download.file read.csv
#' @importFrom tibble tibble
#'
#' @examples
#' \donttest{
#' # Fetch monthly sea ice history from cache if available:
#' icecurves <- get_icecurves()
#' #
#' # Force cache refresh:
#' icecurves <- get_icecurves(use_cache = FALSE)
#' # change region
#' south_icecurves <- get_icecurves(pole='S',  use_cache = FALSE)
#' #
#' # Review cache contents and last update dates:
#' hockeystick_cache_details()
#' #
#' # Plot output using package's built-in ggplot2 settings
#' plot_icecurves(icecurves) }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#' @references
#' \itemize{
#' \item NSIDC Data Archive: \url{https://nsidc.org/data/explore-data}
#' \item All About Sea Ice: \url{https://nsidc.org/learn/parts-cryosphere/sea-ice}
#' \item Sea Ice terminology: \url{https://nsidc.org/learn}
#'  }
#'
#' @export

get_icecurves <- function(pole='N', measure='extent', use_cache = TRUE, write_cache = getOption("hs_write_cache")) {

  if (pole != 'S' & pole != 'N') stop("pole must be 'N' or 'S'")
    if (measure != 'extent' & measure != 'area') stop("measure must be 'extent' or 'area'")

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

  if (use_cache) {
    if (file.exists(file.path(hs_path, 'icecurves.rds')))

      {cached_ice <- readRDS((file.path(hs_path,'icecurves.rds')))
       cached_type <- attr(cached_ice, 'hs_ice_type')

      if (all(c(pole, measure) == cached_type)) return(invisible(cached_ice)) }
  }

  connected <- .isConnected('ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/north/monthly/data/')
  if (!connected) {message("Retrieving remote data requires internet connectivity."); return(invisible(NULL))}

  if (pole=='N') file_url <- 'ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/north/monthly/data/'
  if (pole=='S') file_url <- 'ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/south/monthly/data/'

  month <- c('01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12')

  curve <- function(month) {
    filename <- paste0(pole, '_', month, '_extent_v3.0.csv')

    dl <- tempfile()
    download.file(file.path(file_url, filename), dl)

    seaice <- utils::read.csv(dl, header = TRUE)

    if (measure == 'extent') seaice <- seaice[,c('year', 'mo', 'extent')] else seaice <- seaice[,c('year', 'mo', 'area')]
    return(seaice)
  }

  icecurves <- lapply(month, curve)
  icecurves <- do.call("rbind", icecurves)
  icecurves <- tibble::tibble(icecurves)

    attr(icecurves, "hs_ice_type") <- c(pole, measure)

    if (write_cache) saveRDS(icecurves, file.path(hs_path, 'icecurves.rds'))

  invisible(icecurves) }

#' Download and plot essential climate data
#'
#' Plots the monthly Sea Ice Index data retrieved using `get_icecurves()` with ggplot2. The output ggplot2 object may be further modified.
#'
#'
#' @name plot_icecurves
#' @param dataset Name of the tibble generated by `get_icecurves`, defaults to calling `get_icecurves`
#' @param print (boolean) Display sea ice ggplot2 chart, defaults to TRUE. Use FALSE to not display chart.
#' @return Invisibly returns a ggplot2 object with Sea Ice Index chart
#'
#' @details `plot_icecurves` invisibly returns a ggplot2 object with a pre-defined Sea Ice Index chart using data from `get_icecurves`.
#' By default the chart is also displayed. Users may further modify the output ggplot2 chart.
#' Chart consists of one line per year showing monthly sea ice from January through December. Current year is highlighted.
#'
#' @import ggplot2
#'
#' @examples
#' \donttest{
#' # Fetch historic monthly sea ice data since 1979:
#' icecurves <- get_icecurves()
#' #
#' # Plot output using package's built-in ggplot2 defaults
#' plot_icecurves(icecurves)
#'
#' # Or just call plot_icecurves(), which defaults to get_icecurves() dataset
#' plot_icecurves()
#'
#' p <- plot_icecurves(icecurves, print = FALSE)
#' # Modify plot such as: p + ggplot2::labs(title='Shrinking Arctic Sea Ice') }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#'
#' @export

plot_icecurves <- function(dataset = get_icecurves(), print = TRUE) {

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

  measure_name <- colnames(dataset)[3]
  measure_name <- paste(toupper(substr(measure_name, 1, 1)), substr(measure_name, 2, nchar(measure_name)), sep="")

  ice_region <- if (attr(dataset, "hs_ice_type")[1] == 'S') 'Antarctic' else 'Arctic'

  title <- paste0(ice_region, ' Sea Ice ', measure_name, ' by Year and Month')

  colnames(dataset) <- c('year', 'month', 'measure')
  current_year <- max(dataset$year)
  dataset$linew <- ifelse(dataset$year == current_year, current_year, 'Previous')
  dataset <- dataset[dataset[,3] != '-9999', ]

  plot <- ggplot(dataset, aes(x=month, y=measure, alpha=as.factor(linew), linewidth=as.factor(linew), color=as.factor(year))) +
    geom_line() +scale_y_continuous(limits=c(0, 20)) +theme_bw() +
    scale_x_continuous(breaks = seq(1, 12, 1), minor_breaks = NULL) +
    scale_color_manual(guide = 'none', values = c(rep('darkgrey', length(unique(dataset$year))-1), 'dodgerblue')) +
    scale_linewidth_manual(labels = c(paste0('1979-', current_year-1), current_year), values = c(0.7, 2), breaks = c('Previous', as.character(current_year))) +
    labs(title = title, x = 'Month', y = expression("Million km"^2),
         caption = 'Source: National Snow & Ice Data Center\nhttps://nsidc.org/data/seaice_index',
         linewidth = NULL) +
    theme(legend.position = c(0.41, 0.91), legend.background = element_blank()) +
    scale_alpha_manual(guide = 'none', labels = c('Previous', as.character(current_year)), values = c(0.4, 1), breaks = c('Previous', '2023')) +
    guides(linewidth = guide_legend(override.aes=list(color = c('darkgrey', "dodgerblue"))))


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