inst/doc/detection_and_visualisation.R

In heatwaveR: Detect Heatwaves and Cold-Spells

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

## ----file-structure-----------------------------------------------------------
head(heatwaveR::sst_WA)

## ----detect-example1----------------------------------------------------------
library(dplyr)
library(ggplot2)
library(heatwaveR)

# Detect the events in a time series
ts <- ts2clm(sst_WA, climatologyPeriod = c("1982-01-01", "2011-12-31"))
mhw <- detect_event(ts)

# View just a few metrics
mhw$event %>% 
  dplyr::ungroup() %>%
  dplyr::select(event_no, duration, date_start, date_peak, intensity_max, intensity_cumulative) %>% 
  dplyr::arrange(-intensity_max) %>% 
  head(5)

## ----fig-example1, echo = TRUE, eval = TRUE-----------------------------------
event_line(mhw, spread = 180, metric = "intensity_max", 
           start_date = "1982-01-01", end_date = "2014-12-31")

lolli_plot(mhw, metric = "intensity_max")

## ----fig-example2, echo = TRUE, eval = TRUE-----------------------------------
# Select the region of the time series of interest
mhw2 <- mhw$climatology %>% 
  slice(10580:10720)

ggplot(mhw2, aes(x = t, y = temp, y2 = thresh)) +
  geom_flame() +
  geom_text(aes(x = as.Date("2011-02-25"), y = 25.8, label = "the Destroyer\nof Kelps"))

ggplot(mhw$event, aes(x = date_start, y = intensity_max)) +
  geom_lolli(colour = "salmon", colour_n = "red", n = 3) +
  geom_text(colour = "black", aes(x = as.Date("2006-08-01"), y = 5,
                label = "The marine heatwaves\nTend to be left skewed in a\nGiven time series")) +
  labs(y = expression(paste("Max. intensity [", degree, "C]")), x = NULL)

## ----fig-example3, echo = TRUE, eval = TRUE-----------------------------------
# It is necessary to give geom_flame() at least one row on either side of 
# the event in order to calculate the polygon corners smoothly
mhw_top <- mhw2 %>% 
  slice(5:111)

ggplot(data = mhw2, aes(x = t)) +
  geom_flame(aes(y = temp, y2 = thresh, fill = "all"), show.legend = T) +
  geom_flame(data = mhw_top, aes(y = temp, y2 = thresh, fill = "top"),  show.legend = T) +
  geom_line(aes(y = temp, colour = "temp")) +
  geom_line(aes(y = thresh, colour = "thresh"), size = 1.0) +
  geom_line(aes(y = seas, colour = "seas"), size = 1.2) +
  scale_colour_manual(name = "Line Colour",
                      values = c("temp" = "black", 
                                 "thresh" =  "forestgreen", 
                                 "seas" = "grey80")) +
  scale_fill_manual(name = "Event Colour", 
                    values = c("all" = "salmon", 
                               "top" = "red")) +
  scale_x_date(date_labels = "%b %Y") +
  guides(colour = guide_legend(override.aes = list(fill = NA))) +
  labs(y = expression(paste("Temperature [", degree, "C]")), x = NULL)

## ----flame_gaps---------------------------------------------------------------
mhw3 <- mhw$climatology %>% 
  slice(850:950)

ggplot(mhw3, aes(x = t, y = temp, y2 = thresh)) +
  geom_flame(fill = "black", alpha = 0.5) +
  # Note the use of n = 5 and n_gap = 2 below
  geom_flame(n = 5, n_gap = 2, fill = "red", alpha = 0.5) +
  ylim(c(22, 25)) +
    geom_text(colour = "black", aes(x = as.Date("1984-05-16"), y = 24.5,
                label = "heat\n\n\n\n\nspike"))

## ----fig-example4, echo = TRUE, eval = TRUE-----------------------------------
ggplot(mhw$event, aes(x = date_peak, y = intensity_max)) +
  geom_lolli(colour = "firebrick") +
  labs(x = "Peak Date", 
       y = expression(paste("Max. intensity [", degree, "C]")), x = NULL) +
  theme_linedraw()

## ----detect-example2----------------------------------------------------------
# First calculate the cold-spells
ts_10th <- ts2clm(sst_WA, climatologyPeriod = c("1982-01-01", "2011-12-31"), pctile = 10)
mcs <- detect_event(ts_10th, coldSpells = TRUE)

# Then look at the top few events
mcs$event %>% 
  dplyr::ungroup() %>%
  dplyr::select(event_no, duration, date_start,
                date_peak, intensity_mean, intensity_max, intensity_cumulative) %>%
  dplyr::arrange(intensity_cumulative) %>% 
  head(5)

## ----fig-example5, echo = TRUE, eval = TRUE-----------------------------------
event_line(mcs, spread = 200, metric = "intensity_cumulative",
           start_date = "1982-01-01", end_date = "2014-12-31")

lolli_plot(mcs, metric = "intensity_cumulative", xaxis = "event_no")

## ----fig-example6, echo = TRUE, eval = TRUE-----------------------------------
# Select the region of the time series of interest
mcs2 <- mcs$climatology %>% 
  slice(2900:3190)

# Note that one must specify a colour other than the default 'salmon'
ggplot(mcs2, aes(x = t, y = thresh, y2 = temp)) +
  geom_flame(fill = "steelblue3")

ggplot(mcs$event, aes(x = date_start, y = intensity_max)) +
  geom_lolli(colour = "steelblue3", colour_n = "navy", n = 3) +
  labs(x = "Start Date",
       y = expression(paste("Max. intensity [", degree, "C]")))

## ----fig-example7, echo = TRUE, eval = TRUE-----------------------------------
mcs_top <- mcs2 %>% 
  slice(125:202)

ggplot(data = mcs2, aes(x = t)) +
  geom_flame(aes(y = thresh, y2 = temp, fill = "all"), show.legend = T) +
  geom_flame(data = mcs_top, aes(y = thresh, y2 = temp, fill = "top"), show.legend = T) +
  geom_line(aes(y = temp, colour = "temp")) +
  geom_line(aes(y = thresh, colour = "thresh"), size = 1.0) +
  geom_line(aes(y = seas, colour = "seas"), size = 1.2) +
  scale_colour_manual(name = "Line Colour",
                      values = c("temp" = "black", "thresh" =  "forestgreen", "seas" = "grey80")) +
  scale_fill_manual(name = "Event Colour", values = c("all" = "steelblue3", "top" = "navy")) +
  scale_x_date(date_labels = "%b %Y") +
  guides(colour = guide_legend(override.aes = list(fill = NA))) +
  labs(y = expression(paste("Temperature [", degree, "C]")), x = NULL)

ggplot(mcs$event, aes(x = date_start, y = intensity_cumulative)) +
  geom_lolli(colour = "steelblue3", colour_n = "navy", n = 7) +
  labs( x = "Start Date", y = expression(paste("Cumulative intensity [days x ", degree, "C]")))

## ----plotly-example, eval=FALSE-----------------------------------------------
#  # Must load plotly library first
#  library(plotly)
#  
#  # Function needed for making geom_flame() work with plotly
#  geom2trace.GeomFlame <- function (data,
#                                    params,
#                                    p) {
#  
#    x <- y <- y2 <- NULL
#  
#    # Create data.frame for ease of use
#    data1 <- data.frame(x = data[["x"]],
#                        y = data[["y"]],
#                        y2 = data[["y2"]])
#  
#    # Grab parameters
#    n <- params[["n"]]
#    n_gap <- params[["n_gap"]]
#  
#    # Find events that meet minimum length requirement
#    data_event <- heatwaveR::detect_event(data1, x = x, y = y,
#                                          seasClim = y,
#                                          threshClim = y2,
#                                          minDuration = n,
#                                          maxGap = n_gap,
#                                          protoEvents = T)
#  
#    # Detect spikes
#    data_event$screen <- base::ifelse(data_event$threshCriterion == FALSE, FALSE,
#                                      ifelse(data_event$event == FALSE, TRUE, FALSE))
#  
#    # Screen out spikes
#    data1 <- data1[data_event$screen != TRUE,]
#  
#    # Prepare to find the polygon corners
#    x1 <- data1$y
#    x2 <- data1$y2
#  
#    # # Find points where x1 is above x2.
#    above <- x1 > x2
#    above[above == TRUE] <- 1
#    above[is.na(above)] <- 0
#  
#    # Points always intersect when above=TRUE, then FALSE or reverse
#    intersect.points <- which(diff(above) != 0)
#  
#    # Find the slopes for each line segment.
#    x1.slopes <- x1[intersect.points + 1] - x1[intersect.points]
#    x2.slopes <- x2[intersect.points + 1] - x2[intersect.points]
#  
#    # # Find the intersection for each segment.
#    x.points <- intersect.points + ((x2[intersect.points] - x1[intersect.points]) / (x1.slopes - x2.slopes))
#    y.points <- x1[intersect.points] + (x1.slopes * (x.points - intersect.points))
#  
#    # Coerce x.points to the same scale as x
#    x_gap <- data1$x[2] - data1$x[1]
#    x.points <- data1$x[intersect.points] + (x_gap*(x.points - intersect.points))
#  
#    # Create new data frame and merge to introduce new rows of data
#    data2 <- data.frame(y = c(data1$y, y.points), x = c(data1$x, x.points))
#    data2 <- data2[order(data2$x),]
#    data3 <- base::merge(data1, data2, by = c("x","y"), all.y = T)
#    data3$y2[is.na(data3$y2)] <- data3$y[is.na(data3$y2)]
#  
#    # Remove missing values for better plotting
#    data3$y[data3$y < data3$y2] <- NA
#    missing_pos <- !stats::complete.cases(data3[c("x", "y", "y2")])
#    ids <- cumsum(missing_pos) + 1
#    ids[missing_pos] <- NA
#  
#    # Get the correct positions
#    positions <- data.frame(x = c(data3$x, rev(data3$x)),
#                            y = c(data3$y, rev(data3$y2)),
#                            ids = c(ids, rev(ids)))
#  
#    # Convert to a format geom2trace is happy with
#    positions <- plotly::group2NA(positions, groupNames = "ids")
#    positions <- positions[stats::complete.cases(positions$ids),]
#    positions <- dplyr::left_join(positions, data[,-c(2,3)], by = "x")
#    if(length(stats::complete.cases(positions$PANEL)) > 1)
#      positions$PANEL <- positions$PANEL[stats::complete.cases(positions$PANEL)][1]
#    if(length(stats::complete.cases(positions$group)) > 1)
#      positions$group <- positions$group[stats::complete.cases(positions$group)][1]
#  
#    # Run the plotly polygon code
#    if(length(unique(positions$PANEL)) == 1){
#      getFromNamespace("geom2trace.GeomPolygon", asNamespace("plotly"))(positions)
#    } else{
#      return()
#    }
#  }
#  
#  # Time series
#  ts_res <- heatwaveR::ts2clm(data = heatwaveR::sst_WA,
#                              climatologyPeriod = c("1982-01-01", "2011-12-31"))
#  ts_res_sub <- ts_res[10500:10800,]
#  
#  # Flame Figure
#  p <- ggplot(data = ts_res_sub, aes(x = t, y = temp)) +
#    heatwaveR::geom_flame(aes(y2 = thresh), n = 5, n_gap = 2) +
#    geom_line(aes(y = temp)) +
#    geom_line(aes(y = seas), colour = "green") +
#    geom_line(aes(y = thresh), colour = "red") +
#    labs(x = "", y = "Temperature (°C)")
#  
#  # Create interactive visuals
#  ggplotly(p)