inst/doc/gridded_event_detection.R

In heatwaveR: Detect Heatwaves and Cold-Spells

## ----global_options, include = FALSE------------------------------------------
knitr::opts_chunk$set(fig.width = 4, fig.align = 'center',
                      echo = TRUE, warning = FALSE, message = FALSE, 
                      eval = FALSE, tidy = FALSE)

## ----load-pkg-----------------------------------------------------------------
#  library(dplyr) # For basic data manipulation
#  library(ggplot2) # For visualising data
#  library(heatwaveR) # For detecting MHWs
#  library(tidync) # For easily dealing with NetCDF data
#  library(doParallel) # For parallel processing

## ----load-data----------------------------------------------------------------
#  OISST <- readRDS("~/Desktop/OISST_vignette.Rds")

## ----detect-func--------------------------------------------------------------
#  event_only <- function(df){
#    # First calculate the climatologies
#    clim <- ts2clm(data = df, climatologyPeriod = c("1982-01-01", "2011-01-01"))
#    # Then the events
#    event <- detect_event(data = clim)
#    # Return only the event metric dataframe of results
#    return(event$event)
#  }

## ----detect-purrr-------------------------------------------------------------
#  system.time(
#  # First we start by choosing the 'OISST' dataframe
#  MHW_dplyr <- OISST %>%
#    # Then we group the data by the 'lon' and 'lat' columns
#    group_by(lon, lat) %>%
#    # Then we run our MHW detecting function on each group
#    group_modify(~event_only(.x))
#  ) # ~123 seconds

## ----detect-plyr--------------------------------------------------------------
#  # NB: One should never use ALL available cores, save at least 1 for other essential tasks
#  # The computer I'm writing this vignette on has 8 cores, so I use 7 here
#  registerDoParallel(cores = 7)
#  
#  # Detect events
#  system.time(
#  MHW_plyr <- plyr::ddply(.data = OISST, .variables = c("lon", "lat"), .fun = event_only, .parallel = TRUE)
#  ) # 33 seconds

## ----lon-files----------------------------------------------------------------
#  for(i in 1:length(unique(OISST$lon))){
#    OISST_sub <- OISST %>%
#      filter(lon == unique(lon)[i])
#    saveRDS(object = OISST_sub, file = paste0("~/Desktop/OISST_lon_",i,".Rds"))
#  }

## ----detect-both--------------------------------------------------------------
#  # The 'dplyr' wrapper function to pass to 'plyr'
#  dplyr_wraper <- function(file_name){
#    MHW_dplyr <- readRDS(file_name) %>%
#      group_by(lon, lat) %>%
#      group_modify(~event_only(.x))
#  }
#  # Create a vector of the files we want to use
#  OISST_files <- dir("~/Desktop", pattern = "OISST_lon_*", full.names = T)
#  
#  # Use 'plyr' technique to run 'dplyr' technique with multiple cores
#  system.time(
#  MHW_result <- plyr::ldply(OISST_files, .fun = dplyr_wraper, .parallel = T)
#  ) # 31 seconds
#  
#  # Save for later use as desired
#  saveRDS(MHW_result, "~/Desktop/MHW_result.Rds")

## ----event-tally--------------------------------------------------------------
#  # summarise the number of unique longitude, latitude and year combination:
#  OISST_n <- MHW_result %>%
#    mutate(year = lubridate::year(date_start)) %>%
#    group_by(lon, lat, year) %>%
#    summarise(n = n(), .groups = "drop") %>%
#    group_by(lon, lat) %>%
#    tidyr::complete(year = c(1982:2019)) %>% # Note that these dates may differ
#    mutate(n = ifelse(is.na(n), 0, n))
#  head(OISST_n)

## ----trend-fun----------------------------------------------------------------
#  lin_fun <- function(ev) {
#    mod1 <- glm(n ~ year, family = poisson(link = "log"), data = ev)
#    # extract slope coefficient and its p-value
#    tr <- data.frame(slope = summary(mod1)$coefficients[2,1],
#                     p = summary(mod1)$coefficients[2,4])
#    return(tr)
#  }

## ----apply-trend-fun-plyr-----------------------------------------------------
#  OISST_nTrend <- plyr::ddply(OISST_n, c("lon", "lat"), lin_fun, .parallel = T)
#  OISST_nTrend$pval <- cut(OISST_nTrend$p, breaks = c(0, 0.001, 0.01, 0.05, 1))
#  head(OISST_nTrend)

## ----prep-geography-----------------------------------------------------------
#  # The base map
#  map_base <- ggplot2::fortify(maps::map(fill = TRUE, plot = FALSE)) %>%
#    dplyr::rename(lon = long)

## ----the-figure---------------------------------------------------------------
#  map_slope <- ggplot(OISST_nTrend, aes(x = lon, y = lat)) +
#    geom_rect(size = 0.2, fill = NA,
#         aes(xmin = lon - 0.1, xmax = lon + 0.1, ymin = lat - 0.1, ymax = lat + 0.1,
#             colour = pval)) +
#    geom_raster(aes(fill = slope), interpolate = FALSE, alpha = 0.9) +
#    scale_fill_gradient2(name = "count/year (slope)", high = "red", mid = "white",
#                         low = "darkblue", midpoint = 0,
#                         guide = guide_colourbar(direction = "horizontal",
#                                                 title.position = "top")) +
#    scale_colour_manual(breaks = c("(0,0.001]", "(0.001,0.01]", "(0.01,0.05]", "(0.05,1]"),
#                        values = c("firebrick1", "firebrick2", "firebrick3", "white"),
#                        name = "p-value", guide = FALSE) +
#    geom_polygon(data = map_base, aes(group = group),
#                 colour = NA, fill = "grey80") +
#    coord_fixed(ratio = 1, xlim = c(13.0, 23.0), ylim = c(-33, -42), expand = TRUE) +
#    labs(x = "", y = "") +
#    theme_bw() +
#    theme(legend.position = "bottom")
#  
#  map_p <- ggplot(OISST_nTrend, aes(x = lon, y = lat)) +
#    geom_raster(aes(fill = pval), interpolate = FALSE) +
#    scale_fill_manual(breaks = c("(0,0.001]", "(0.001,0.01]", "(0.01,0.05]",
#                                 "(0.05,0.1]", "(0.1,0.5]", "(0.5,1]"),
#                      values = c("black", "grey20", "grey40",
#                                 "grey80", "grey90", "white"),
#                      name = "p-value",
#                      guide = guide_legend(direction = "horizontal",
#                                                 title.position = "top")) +
#    geom_polygon(data = map_base, aes(group = group),
#                 colour = NA, fill = "grey80") +
#    coord_fixed(ratio = 1, xlim = c(13.0, 23.0), ylim = c(-33, -42), expand = TRUE) +
#    labs(x = "", y = "") +
#    theme_bw() +
#    theme(legend.position = "bottom")
#  
#  map_both <- ggpubr::ggarrange(map_slope, map_p, align = "hv")
#  map_both