R/instrumental_temp.R

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

#' Download and plot essential climate data
#'
#' Retrieves the combined global land- and sea-surface temperature anomaly (Land-Ocean Temperature Index, LOTI).
#' Source is NASA/GISS Surface Temperature Analysis (GISTEMP v4), an estimate of global surface temperature change.
#' \url{https://data.giss.nasa.gov/gistemp/}
#'
#' @name get_temp
#' @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 mean and monthly Combined Land-Surface Air and Sea-Surface Water Temperature Anomalies.
#'
#' `get_temp` invisibly returns a tibble with the NASA/GISS annual mean and monthly global temperature anomaly. Data are global from 1880 to present, and represent the deviations from the 1951-1980 mean.
#'
#' The returned object includes monthly and annual average anomalies, as well as seasonal anomalies.
#' GISS Surface Temperature Analysis (GISTEMP v4) is an estimate of global surface temperature change.
#'
#' Data are updated around the middle of every month using current data files from NOAA GHCN v4 (meteorological stations) and ERSST v5 (ocean areas).
#' Station data are combined as described in Hansen et al. (2010) \url{https://data.giss.nasa.gov/gistemp/references.html} and Lenssen et al. (2019) \url{https://pubs.giss.nasa.gov/abs/le05800h.html}
#'
#' @importFrom readr read_csv
#' @importFrom lubridate ymd
#' @importFrom utils download.file
#'
#' @examples
#' \donttest{
#' # Fetch temp anomaly from cache if available:
#' anomaly <- get_temp()
#' #
#' # Force cache refresh:
#' anomaly <- get_temp(use_cache=FALSE)
#' #
#' # Review cache contents and last update dates:
#' hockeystick_cache_details()
#' #
#' # Plot output using package's built-in ggplot2 settings
#' plot_temp(anomaly) }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#' @references
#' \itemize{
#' \item GISS Surface Temperature Analysis (GISTEMP v4): \url{https://data.giss.nasa.gov/gistemp/}
#' \item GISTEMP Team, 2020: \emph{GISS Surface Temperature Analysis (GISTEMP), version 4.} NASA Goddard Institute for Space Studies.
#' \item Lenssen, N., G. Schmidt, J. Hansen, M. Menne, A. Persin, R. Ruedy, and D. Zyss, 2019:  Improvements in the GISTEMP uncertainty model.
#'  \emph{J. Geophys. Res. Atmos.}, 124, no. 12, 6307-6326, doi:10.1029/2018JD029522. \url{https://pubs.giss.nasa.gov/abs/le05800h.html}
#'  }
#'
#' @export
get_temp <- 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,'gisstemp.rds'))) return(invisible(readRDS((file.path(hs_path,'gisstemp.rds')))))
  }

file_url <- 'https://data.giss.nasa.gov/gistemp/tabledata_v4/GLB.Ts+dSST.csv'
connected <- .isConnected(file_url)
if (!connected) {message("Retrieving remote data requires internet connectivity."); return(invisible(NULL))}

dl <- tempfile()
download.file(file_url, dl)
gisstemp <- suppressMessages( readr::read_csv(dl, skip=1, na='***') )
# Compute year-to-date mean for current year
gisstemp[nrow(gisstemp), 'J-D'] <- mean(as.numeric(gisstemp[nrow(gisstemp), 2:13]), na.rm = T)
gisstemp$Year <- ymd(paste(gisstemp$Year,'12','31',sep='-'))
dir.create(hs_path, showWarnings = FALSE, recursive = TRUE)
if (write_cache) saveRDS(gisstemp, file.path(hs_path, 'gisstemp.rds'))

invisible(gisstemp)
}



#' Download and plot essential climate data
#'
#' Plots the global annual mean temperature anomaly retrieved using `get_temp()` with ggplot2. The output ggplot2 object may be further modified.
#'
#'
#' @name plot_temp
#' @param dataset Name of the tibble generated by `get_temp`
#' @param print (boolean) Display temperature anomaly ggplot2 chart, defaults to TRUE. Use FALSE to not display chart.
#'
#' @return Invisibly returns a ggplot2 object with temperature anomaly chart
#'
#' @details `plot_temp` invisibly returns a ggplot2 object with a pre-defined temperature annual mean anomaly chart using data from `get_temp`.
#' By default the chart is also displayed. Smooths using ggplot2's built-in loess smoother. Users may further modify the output ggplot2 chart.
#'
#' @importFrom ggplot2 ggplot
#'
#' @examples
#' \donttest{
#' # Fetch temperature anomaly:
#' anomaly <- get_temp()
#' #
#' # Plot output using package's built-in ggplot2 defaults
#' plot_temp(anomaly)
#'
#' # Or just call plot_temp(), which defaults to get_temp() dataset
#' plot_temp()
#'
#' p <- plot_temp(anomaly, print = FALSE)
#' # Modify plot such as: p + ggplot2::labs(title='The Signature of Climate Change') }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#'
#' @export

plot_temp <- function(dataset = get_temp(), print=TRUE) {

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

plot <- ggplot(dataset, aes(x=Year, y=`J-D`)) +geom_line(alpha=0.75, aes(color='Annual mean')) + theme_bw(base_size=12) +
  scale_x_date(name=NULL, limits=c(as.Date('1878-01-01'), ymd(max(dataset$Year))), date_breaks='15 years', date_labels='%Y') +
  scale_y_continuous(n.breaks = 9) +geom_smooth(size=1.1, se=F, span=0.3, aes(color='Loess smoothing')) +
  labs(title='Global Land-Ocean Temperature Index (LOTI)', subtitle='Global surface temperature relative to 1951-1980 mean',
       y='Temperature Anomaly (C\U00B0)', caption='Source: NASA Goddard Institute for Space Studies\nhttps://data.giss.nasa.gov/gistemp/') +
  scale_color_manual(name=NULL, values=c('dodgerblue2','firebrick1')) +theme(legend.position = c(0.175, 0.825),legend.background=element_blank())

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



#' Download and plot essential climate data
#'
#' Plots the monthly mean temperature anomaly retrieved using `get_temp()` with ggplot2. The output ggplot2 object may be further modified.
#'
#'
#' @name plot_temp_monthly
#' @param dataset Name of the tibble generated by `get_temp`
#' @param print (boolean) Display temperature anomaly ggplot2 chart, defaults to TRUE. Use FALSE to not display chart.
#'
#' @return Invisibly returns a ggplot2 object with temperature anomaly chart
#'
#' @details `plot_temp_monthly` invisibly returns a ggplot2 object with a pre-defined temperature monthly mean anomaly chart using data from `get_temp`.
#' By default the chart is also displayed. Smooths using ggplot2's built-in loess smoother. Users may further modify the output ggplot2 chart.
#'
#' @importFrom ggplot2 ggplot
#' @importFrom dplyr mutate select
#' @importFrom lubridate ceiling_date days year
#'
#' @examples
#' \donttest{
#' # Fetch temperature anomaly:
#' anomaly <- get_temp()
#' #
#' # Plot output using package's built-in ggplot2 defaults
#' plot_temp_monthly(anomaly)
#'
#' # Or just call plot_temp_monthly(), which defaults to get_temp() dataset
#' plot_temp_monthly()
#'
#' p <- plot_temp_monthly(anomaly, print = FALSE)
#' # Modify plot such as: p + ggplot2::labs(title='The Signature of Climate Change') }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#'
#' @export

plot_temp_monthly <- function(dataset = get_temp(), print=TRUE) {

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

  dataset <- dataset |> select(1:13) |> tidyr::pivot_longer(cols = 2:13, names_to = 'Month', values_to = 'Anomaly') |>
    filter(!is.na(Anomaly))

  dataset <- dataset |> mutate(Date = ceiling_date(ym(paste(year(Year), Month)), 'month') - days(1))

  plot <- ggplot(dataset, aes(x=Date, y=Anomaly)) +geom_line(alpha=0.75, aes(color='Monthly mean')) + theme_bw(base_size=12) +
    scale_x_date(name=NULL, limits=c(as.Date('1878-01-01'), ymd(max(dataset$Year))), date_breaks='15 years', date_labels='%Y') +
    scale_y_continuous(n.breaks = 10) +geom_smooth(size=1.1, se=F, span=0.3, aes(color='Loess smoothing')) +
    labs(title='Global Land-Ocean Temperature Index (LOTI)', subtitle='Global surface temperature relative to 1951-1980 mean',
         y='Temperature Anomaly (C\U00B0)', caption='Source: NASA Goddard Institute for Space Studies\nhttps://data.giss.nasa.gov/gistemp/') +
    scale_color_manual(name=NULL, values=c('firebrick1','dodgerblue2')) +
    theme(legend.position = c(0.175, 0.825), legend.background=element_blank())

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



#' Download and plot essential climate data
#'
#' Plots scatter of monthly temperature anomaly retrieved using `get_temp()` with ggplot2. The output ggplot2 object may be further modified.
#'
#'
#' @name plot_temp_scatter
#' @param dataset Name of the tibble generated by `get_temp`
#' @param print (boolean) Display temperature anomaly ggplot2 chart, defaults to TRUE. Use FALSE to not display chart.
#' @param labellatest (boolean) Display date of latest value, default = TRUE
#' @param labelmax (boolean) Display date of max value, default = FALSE
#'
#' @return Invisibly returns a ggplot2 object with temperature anomaly chart
#'
#' @details `plot_temp_scatter` invisibly returns a ggplot2 object with a pre-defined temperature monthly mean anomaly chart using data from `get_temp`.
#' By default the chart is also displayed. Smooths using ggplot2's built-in loess smoother. Users may further modify the output ggplot2 chart.
#'
#' @importFrom ggplot2 ggplot
#' @importFrom dplyr mutate select
#' @importFrom lubridate ceiling_date days ym
#'
#' @examples
#' \donttest{
#' # Fetch temperature anomaly:
#' anomaly <- get_temp()
#' #
#' # Plot output using package's built-in ggplot2 defaults
#' plot_temp_scatter(anomaly)
#'
#' # Or just call plot_temp_scatter(), which defaults to get_temp() dataset
#' plot_temp_scatter()
#'
#' p <- plot_temp_scatter(anomaly, print = FALSE)
#' # Modify plot such as: p + ggplot2::labs(title='The Signature of Climate Change') }
#'
#' @author Hernando Cortina, \email{hch@@alum.mit.edu}
#'
#' @export

plot_temp_scatter <- function(dataset = get_temp(), print=TRUE, labelmax=FALSE, labellatest=TRUE) {

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

  dataset <- dataset |> select(1:13) |> tidyr::pivot_longer(cols = 2:13, names_to = 'Month', values_to = 'Anomaly') |>
    filter(!is.na(Anomaly))

  dataset <- dataset |> mutate(Date = ceiling_date(ym(paste(year(Year), Month)), 'month') - days(1))

  col_strip <- RColorBrewer::brewer.pal(11, "RdBu")[-(5:7)]

  plot <- ggplot(dataset, aes(x=Year, y=Anomaly, color=Anomaly)) + geom_point(size=1) + theme_bw(base_size=12) +
    scale_x_date(name=NULL, limits=c(as.Date('1878-01-01'), ymd(max(dataset$Year))), date_breaks='15 years', date_labels='%Y') +
    scale_y_continuous(n.breaks = 8)  +
    labs(title='Global Land-Ocean Temperature Index (LOTI)', subtitle='Global surface temperature relative to 1951-1980 mean', color='',
         y='Monthly Temperature Anomaly (C\U00B0)', caption='Source: NASA Goddard Institute for Space Studies\nhttps://data.giss.nasa.gov/gistemp/') +
    theme(legend.position='inside', legend.position.inside = c(0.09, 0.85),legend.background=element_blank() ) +
    scale_color_gradientn(colors = rev(col_strip))

  if (labelmax) {
    dataset <- dataset |> mutate(max = ifelse(Anomaly == max(dataset$Anomaly), T, F))
    dataset <- dataset |> mutate(name = ifelse(max == T, substr(as.character(Date),1,7), ''))
    plot <- plot + geom_text(aes(y = Anomaly+0.1, label=dataset$name), size=3) }


  if (labellatest) {
    if ( !exists("name", dataset) ) dataset$name <- ''
    dataset[nrow(dataset),'name'] <- substr(as.character(pull(dataset[nrow(dataset),'Date'])), 1, 7)
    plot <- plot + geom_text(aes(y = Anomaly+0.2, label=dataset$name), size=3) }

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