heatwaveR

https://robwschlegel.github.io/heatwaveR/
Admin
heatwaveR provides fast, research-grade tools to detect, classify, and visualise heatwaves and...
Tokens:46,815
Snippets:391
Trust Score:9.7
Update:3 months ago
Show doc for...
Tokens:
Raw
### Load required libraries

Source: https://robwschlegel.github.io/heatwaveR/articles/gridded_event_detection.html

Initial setup for data manipulation, visualization, and MHW detection.

```R
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
```

--------------------------------

### Install heatwaveR development version

Source: https://robwschlegel.github.io/heatwaveR/index.html

Use the remotes package to install the latest development version directly from GitHub.

```R
remotes::install_github("robwschlegel/heatwaveR")
```

--------------------------------

### Generating a Lollipop Plot

Source: https://robwschlegel.github.io/heatwaveR/reference/lolli_plot.html

Example workflow showing data preparation with ts2clm and detect_event followed by the creation of a default lollipop plot.

```R
ts <- ts2clm(sst_WA, climatologyPeriod = c("1983-01-01", "2012-12-31"))
res <- detect_event(ts)

library(ggplot2)

# The default output
lolli_plot(res)
```

--------------------------------

### Load required R packages for OISST data processing

Source: https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html

Installs and loads necessary libraries for data manipulation, visualization, and downloading. Uncomment the install commands if the packages are not already present in your environment.

```R
# The packages we will need
# install.packages("dplyr")
# install.packages("lubridate")
# install.packages("ggplot2")
# install.packages("tidync")
# install.packages("doParallel")
# install.packages("rerddap")
# install.packages("plyr") # Note that this library should never be loaded, only installed

# The packages we will use
library(dplyr) # A staple for modern data management in R
library(lubridate) # Useful functions for dealing with dates
library(ggplot2) # The preferred library for data visualisation
library(tidync) # For easily dealing with NetCDF data
library(rerddap) # For easily downloading subsets of data
library(doParallel) # For parallel processing
```

--------------------------------

### Install heatwaveR from CRAN

Source: https://robwschlegel.github.io/heatwaveR/index.html

Use this command to install the stable version of the package from the Comprehensive R Archive Network.

```R
install.packages("heatwaveR")
```

--------------------------------

### Visualize event metrics with geom_lolli

Source: https://robwschlegel.github.io/heatwaveR/reference/geom_lolli.html

Examples demonstrating how to use geom_lolli to visualize marine heatwave event data, including mapping aesthetics and highlighting top events.

```R
ts <- ts2clm(sst_WA, climatologyPeriod = c("1983-01-01", "2012-12-31"))
res <- detect_event(ts)
mhw <- res$event

library(ggplot2)

# Height of lollis represent event durations and their colours
# are mapped to the events' cumulative intensity:
ggplot(mhw, aes(x = date_peak, y = duration)) +
  geom_lolli(aes(colour = intensity_cumulative)) +
  scale_color_distiller(palette = "Spectral", name = "Cumulative \nintensity") +
  xlab("Date") + ylab("Event duration [days]")


# Height of lollis represent event durations and the top three (longest)
# lollis are highlighted in red:
ggplot(mhw, aes(x = date_peak, y = duration)) +
  geom_lolli(n = 3, colour_n = "red") +
  scale_color_distiller(palette = "Spectral") +
  xlab("Peak date") + ylab("Event duration [days]")


# Because this is a proper geom, any number of ill-advised things
# may be done with it:
ggplot(mhw, aes(x = event_no, y = intensity_max)) +
  geom_lolli(shape = 5, aes(colour = rate_onset), linetype = "dotted") +
  scale_color_distiller(palette = "RdYlGn", name = "Rate \nonset") +
  xlab("Event number") + ylab("Max intensity [degree C]")
```

--------------------------------

### Calculate yearly averages and visualize event counts

Source: https://robwschlegel.github.io/heatwaveR/reference/block_average.html

This example demonstrates the workflow from climatology generation and event detection to calculating yearly averages and visualizing the results with ggplot2.

```R
ts <- ts2clm(sst_WA, climatologyPeriod = c("1983-01-01", "2012-12-31"))
res <- detect_event(ts)
out <- block_average(res)
summary(glm(count ~ year, out, family = "poisson"))
#> 
#> Call:
#> glm(formula = count ~ year, family = "poisson", data = out)
#> 
#> Coefficients:
#>              Estimate Std. Error z value Pr(>|z|)   
#> (Intercept) -63.70416   20.29044  -3.140  0.00169 **
#> year          0.03209    0.01011   3.173  0.00151 **
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 116.78  on 40  degrees of freedom
#> Residual deviance: 106.27  on 39  degrees of freedom
#> AIC: 176.77
#> 
#> Number of Fisher Scoring iterations: 6
#> 

library(ggplot2)

ggplot(data = out, aes(x = year, y = count)) +
  geom_point(colour = "salmon") +
  geom_line() +
  labs(x = NULL, y = "Number of events")
```

--------------------------------

### Visualize marine heatwaves with geom_flame

Source: https://robwschlegel.github.io/heatwaveR/reference/geom_flame.html

Example demonstrating how to prepare data using ts2clm and detect_event, then visualize the results using ggplot2 and geom_flame.

```R
ts <- ts2clm(sst_WA, climatologyPeriod = c("1983-01-01", "2012-12-31"))
res <- detect_event(ts)
mhw <- res$clim
mhw <- mhw[10580:10690,]

library(ggplot2)

ggplot(mhw, aes(x = t, y = temp)) +
  geom_flame(aes(y2 = thresh)) +
  annotate(geom = "text", x = as.Date("2011-02-01"), y = 28,
           label = "That's not a heatwave.\nThis, is a heatwave.") +
  xlab("Date") + ylab(expression(paste("Temperature [", degree, "C]")))
```

--------------------------------

### Define Time Ranges for Data Download

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_metric_trends.html

Creates a data frame specifying the start and end dates for multi-year data retrieval.

```R
dl_years <- data.frame(date_index = 1:5,
                       start = as.Date(c("1982-01-01", "1990-01-01",
                                         "1998-01-01", "2006-01-01", "2014-01-01")),
                       end = as.Date(c("1989-12-31", "1997-12-31",
                                       "2005-12-31", "2013-12-31", "2021-12-31")))
```

--------------------------------

### Load required libraries for MHW analysis

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_metric_trends.html

Initializes the environment with necessary packages for data manipulation, visualization, and statistical trend analysis.

```R
library(heatwaveR)
library(lubridate)
library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(viridis)
library(trend)
library(kableExtra)
library(rerddap)
library(sf)
```

--------------------------------

### Load Required Libraries

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_to_nc.html

Initializes the necessary packages for data manipulation, NetCDF creation, and parallel processing.

```R
library(dplyr) # For basic data manipulation
library(ncdf4) # For creating NetCDF files
library(tidync) # For easily dealing with NetCDF data
library(ggplot2) # For visualising data
library(doParallel) # For parallel processing
```

--------------------------------

### make_whole_fast(data)

Source: https://robwschlegel.github.io/heatwaveR/reference/make_whole_fast.html

Takes a series of dates and temperatures, and if irregular (but ordered), inserts missing dates and fills corresponding temperatures with NAs.

```APIDOC
## make_whole_fast(data)

### Description
Constructs a continuous, uninterrupted time series of temperatures. It takes a series of dates and temperatures, and if irregular (but ordered), inserts missing dates and fills corresponding temperatures with NAs.

### Arguments
- **data** (data.frame) - Required - A data frame with columns for date (ts_x) and temperature (ts_y) data. Ordered daily data are expected.

### Value
Returns a data frame with three columns: 
- **doy** (integer) - Day-of-year (Julian day 1-366, adjusted for leap years).
- **ts_x** (Date) - A continuous series of dates.
- **y** (numeric) - The measured variable (temperature).
```

--------------------------------

### Inspect sample dataset

Source: https://robwschlegel.github.io/heatwaveR/articles/detection_and_visualisation.html

Displays the first six rows of the sst_WA dataset to verify the required date and temperature column structure.

```R
head(heatwaveR::sst_WA)
```

--------------------------------

### Create and Populate NetCDF File

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_to_nc.html

Define variables, create the NetCDF file, write data arrays, and add global metadata attributes.

```R
# Path and file name
  # NB: We are net setting time dimensions here because we are using event_no as our "time" dimension
ncpath <- "~/Desktop/"
ncname <- "MHW_results"
ncfname <- paste0(ncpath, ncname, ".nc")
# dname <- "tmp"

# Define dimensions
londim <- ncdim_def("lon", "degrees_east", lon)
latdim <- ncdim_def("lat", "degrees_north", lat)
endim <- ncdim_def("event_no", "event_number", event_no)
# timedim <- ncdim_def("time", tunits, as.double(time))

# Define variables
fillvalue <- 1e32
def_dur <- ncvar_def(name = "duration", units = "days", dim = list(endim, latdim, londim), 
                     missval = fillvalue, longname = "duration of MHW", prec = "double")
def_max_int <- ncvar_def(name = "max_int", units = "deg_c", dim = list(endim, latdim, londim), 
                         missval = fillvalue, longname = "maximum intensity during MHW", prec = "double")
def_cum_int <- ncvar_def(name = "cum_int", units = "deg_c days", dim = list(endim, latdim, londim), 
                         missval = fillvalue, longname = "cumulative intensity during MHW", prec = "double")
def_peak <- ncvar_def(name = "date_peak", units = tunits, dim = list(endim, latdim, londim),
                      missval = 0, longname = "date of peak temperature anomaly during MHW", prec = "integer")

# Create netCDF file and prepare space for our arrays
ncout <- nc_create(ncfname, list(def_peak, def_dur, def_max_int, def_cum_int), force_v4 = TRUE)

# Add the actual data
ncvar_put(ncout, def_dur, prep_dur)
ncvar_put(ncout, def_max_int, prep_max_int)
ncvar_put(ncout, def_cum_int, prep_cum_int)
ncvar_put(ncout, def_peak, prep_peak)

# Put additional attributes into dimension and data variables
ncatt_put(ncout, "lon", "axis", "X") #,verbose=FALSE) #,definemode=FALSE)
ncatt_put(ncout, "lat", "axis", "Y")
ncatt_put(ncout, "event_no", "axis", "event_no")

# Add global attributes
  # These are useful for whomever else may want to use your NetCDF file
ncatt_put(ncout, 0, "title", paste0("MHW results from lon: ",
                                    min(lon)," to ",max(lon),
                                    " and lat: ",min(lat)," to ",max(lat)))
ncatt_put(ncout, 0, "institution", "Your institution here!")
ncatt_put(ncout, 0, "source", "NOAA OISST v2.1")
ncatt_put(ncout,0, "references", "Banzon et al. (2020) J. Atmos. Oce. Tech. 37:341-349")
history <- paste0("Your name here!, ", date())
ncatt_put(ncout, 0, "history", history)
ncatt_put(ncout, 0, "Conventions", "Hobday et al. (2016)") # Assuming one has used the Hobday definition

# Get a summary of the created file:
ncout
```

--------------------------------

### Calculate Thresholds and Detect Events

Source: https://robwschlegel.github.io/heatwaveR/articles/complex_clims.html

Demonstrates calculating climatology thresholds for different temperature metrics and detecting events using one or two thresholds.

```R
# Note that because we are not using the standard column name 'temp' we must
# specify the chosen column name twice, once for ts2clm() and again for detect_event()

# First threshold based on tMin
thresh_tMin <- ts2clm(data = Algiers, y = tMin, pctile = 80, 
                      climatologyPeriod = c("1961-01-01", "1990-12-31"))

# Second threshold based on tMax
# Be careful here that you put the arguments within the correct brackets
thresh_tMax <- detect_event(ts2clm(data = Algiers, y = tMax, pctile = 90, 
                                   climatologyPeriod = c("1961-01-01", "1990-12-31")),
                            # These arguments are passed to detect_event(), not ts2clm()
                            minDuration = 3, maxGap = 0, y = tMax, protoEvents = T)

# Detect/calculate events using the two precalculated thresholds
# Because detect_event() is not able to deduce which arguments we used above,
# we must again tell it explicitly here
events_two_thresh <- detect_event(data = thresh_tMin, y = tMin, minDuration = 10, maxGap = 2,
                                  threshClim2 = thresh_tMax$event, minDuration2 = 3, maxGap2 = 0)

# Or to simply use one threshold
events_one_thresh <- detect_event(data = thresh_tMin, y = tMin, minDuration = 10, maxGap = 2)
```

--------------------------------

### Define OISST download date ranges

Source: https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html

Creates a data frame containing start and end dates for batch processing of OISST data.

```R
dl_years <- data.frame(date_index = 1:5,
                       start = c("1982-01-01", "1990-01-01", 
                                 "1998-01-01", "2006-01-01", "2014-01-01"),
                       end = c("1989-12-31", "1997-12-31", 
                               "2005-12-31", "2013-12-31", "2019-12-31"))
```

--------------------------------

### Enable interactive geom_flame with plotly

Source: https://robwschlegel.github.io/heatwaveR/articles/detection_and_visualisation.html

This snippet defines a custom geom2trace.GeomFlame function to bridge heatwaveR's geom_flame with plotly, followed by a usage example.

```R
# 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)
```

--------------------------------

### Load required libraries

Source: https://robwschlegel.github.io/heatwaveR/articles/detection_and_visualisation.html

Loads the necessary packages for data manipulation, visualization, and marine heatwave detection.

```R
library(dplyr)
library(ggplot2)
library(heatwaveR)
```

--------------------------------

### Load and Process SST Data

Source: https://robwschlegel.github.io/heatwaveR/articles/Download_SST_v2.html

Sets up the environment, defines file paths, and creates functions for slicing netCDF data and detecting marine heatwave events.

```R
# Load packages -----------------------------------------------------------

library(tidyverse)
library(tidync)
library(data.table) # for the fast fwrite() function
library(heatwaveR)
library(doParallel)
registerDoParallel(cores = 14) # not using all 16

# Define paths ------------------------------------------------------------

datadir <- "/Volumes/OceanData"
oisst_file_dir <- paste0(datadir, "/test_files")
nc_file <- paste0(oisst_file_dir, "/OISST_combined.nc")
MHW_dir <- datadir

# Define various functions ------------------------------------------------

# A load and slice function for the combined netCDF
OISST_load <- function(file_name, lon1, lon2) {
  OISST_dat <- tidync(file_name) %>%
    hyper_filter(lon = between(lon, lon1, lon2)) %>%
    hyper_tibble(select_var = "sst", force = TRUE, drop = TRUE) %>%
    select(-zlev) %>%
    dplyr::rename(t = time, temp = sst) %>%
    mutate(t = as.Date(t, origin = "1978-01-01"))
  return(OISST_dat)
  rm(OISST_dat)
}

# Rob's MHW detect function
event_only <- function(df) {
  # first calculate the climatologies
  clim <- ts2clm(data = df, climatologyPeriod = c("1991-01-01", "2020-12-31"))
  # then the events
  event <- detect_event(data = clim)
  rm(clim)
  # return only the event metric dataframe of results
  return(event$event)
  rm(event)
}
```

--------------------------------

### Retrieve package citation

Source: https://robwschlegel.github.io/heatwaveR/index.html

Run this command within an R session to display the appropriate citation information for the package.

```R
citation("heatwaveR")
```

--------------------------------

### Usage of clim_spread

Source: https://robwschlegel.github.io/heatwaveR/reference/clim_spread.html

This function is used internally to spread time series data into a wide format based on specified start and end dates and a smoothing window width.

```R
clim_spread(data, clim_start, clim_end, windowHalfWidth)
```

--------------------------------

### Visualize MHW results

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_to_nc.html

Calls the quick_grid function to plot duration and intensity metrics from gridded and NetCDF data sources.

```R
# Plot the duration results
quick_grid(MHW_res_grid, "duration")
quick_grid(MHW_res_nc, "duration")

# Cumulative intensity
quick_grid(MHW_res_grid, "intensity_cumulative")
quick_grid(MHW_res_nc, "cum_int")
```

--------------------------------

### Create a wrapper function for OISST data download

Source: https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html

This function automates the retrieval and cleaning of OISST data for a specified time range and spatial extent.

```R
# This function downloads and prepares data based on user provided start and end dates
OISST_sub_dl <- function(time_df){
  OISST_dat <- rerddap::griddap(datasetx = "ncdcOisst21Agg_LonPM180",
                                url = "https://coastwatch.pfeg.noaa.gov/erddap/", 
                                time = c(time_df$start, time_df$end), 
                                zlev = c(0, 0),
                                latitude = c(-40, -35),
                                longitude = c(15, 21),
                                fields = "sst")$data |> 
    dplyr::mutate(time = base::as.Date(stringr::str_remove(time, "T12:00:00Z"))) |> 
    dplyr::rename(t = time, temp = sst, lon = longitude, lat = latitude) |> 
    dplyr::select(lon, lat, t, temp) |> 
    stats::na.omit()
}
```

--------------------------------

### Prepare base map for visualisation

Source: https://robwschlegel.github.io/heatwaveR/articles/gridded_event_detection.html

Extracts and formats global map data for plotting.

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

--------------------------------

### Usage of make_whole_fast

Source: https://robwschlegel.github.io/heatwaveR/reference/make_whole_fast.html

Basic syntax for calling the make_whole_fast function with a data frame containing date and temperature columns.

```R
make_whole_fast(data)
```

--------------------------------

### Execute Batch Data Download

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_metric_trends.html

Iterates through the defined time ranges to download and combine the full dataset.

```R
OISST_data <- dl_years %>% 
    group_by(date_index) %>% 
    group_modify(~OISST_sub_dl(.x)) %>% 
    ungroup() %>% 
    dplyr::select(lon, lat, t, temp)
```

--------------------------------

### Generate event line plots

Source: https://robwschlegel.github.io/heatwaveR/reference/event_line.html

Creates standalone line plots visualizing temperature events. Requires output from detect_event as input.

```R
ts <- ts2clm(sst_WA, climatologyPeriod = c("1983-01-01", "2012-12-31"))
res <- detect_event(ts)

event_line(res, spread = 100, metric = duration,
start_date = "2010-12-01", end_date = "2011-06-30")


event_line(res, spread = 100, start_date = "2010-12-01",
end_date = "2011-06-30", category = TRUE)


event_line(res, spread = 100, start_date = "2010-12-01",
end_date = "2011-06-30", category = TRUE,
line_colours = c("black", "blue", "gray20", "gray20", "gray20", "gray20"))
```

--------------------------------

### Execute dplyr method

Source: https://robwschlegel.github.io/heatwaveR/articles/gridded_event_detection.html

Initiate the event detection process using the dplyr method.

```R
system.time(
```

--------------------------------

### Define OISST Data URLs and File Paths

Source: https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html

Sets the base URL and generates a data frame of file paths for a specific date range.

```R
# First we tell R where the data are on the interwebs
OISST_base_url <- "https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/v2.1/access/avhrr/"
# Note that one may go to this URL in any web browser to manually inspect the files

# Now we create a data.frame that contains all of the dates we want to download
  # NB: In order to change the dates download changes the dates in the following line
OISST_dates <- base::data.frame(t = seq(as.Date("2019-12-01"), as.Date("2019-12-31"), by = "day"))

# To finish up this step we add some text to those dates so they match the OISST file names
OISST_files <- OISST_dates |> 
  dplyr::mutate(t_day = base::gsub("-", "", t),
                t_month = base::substr(t_day, 1, 6),
                t_year = lubridate::year(t),
                file_name = base::paste0(OISST_base_url, t_month, "/", "oisst-avhrr-v02r01.", t_day ,".nc"))
```

--------------------------------

### Download OISST NetCDF Files

Source: https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html

Defines a download function and executes it in parallel to fetch daily NetCDF files.

```R
# This function will go about downloading each day of data as a NetCDF file
# Note that this will download files into a 'data/OISST' folder in the root directory
  # If this folder does not exist it will create it
  # If it does not automatically create the folder it will need to be done manually
  # The folder that is created must be a new folder with no other files in it
  # A possible bug with netCDF files in R is they won't load correctly from 
  # existing folders with other file types in them
# This function will also check if the file has been previously downloaded
  # If it has it will not download it again
OISST_url_daily_dl <- function(target_URL){
  base::dir.create("~/data/OISST", showWarnings = F)
  file_name <- base::paste0("~/data/OISST/",base::sapply(base::strsplit(target_URL, split = "/"), "[[", 10))
  if(!base::file.exists(file_name)) utils::download.file(url = target_URL, method = "libcurl", destfile = file_name)
}

# The more cores used, the faster the data may be downloaded
  # It is best practice to not use all of the cores on one's machine
  # The laptop on which I am running this code has 8 cores, so I use 7 here
doParallel::registerDoParallel(cores = 7)

# And with that we are clear for take off
base::system.time(plyr::l_ply(OISST_files$file_name, .fun = OISST_url_daily_dl, .parallel = T)) # ~15 seconds

# In roughly 15 seconds a user may have a full month of global data downloaded
# This scales well into years and decades, and is much faster with more cores
# Download speeds will also depend on the speed of the users internet connection
```

--------------------------------

### Accessing the sst_WA dataset

Source: https://robwschlegel.github.io/heatwaveR/reference/sst_WA.html

Load the sst_WA dataframe to access temperature and date information for the Western Australian region.

```R
sst_WA
```

--------------------------------

### Internal Workflow Visualization

Source: https://robwschlegel.github.io/heatwaveR/articles/internal_workflow.html

A Mermaid flowchart illustrating the data flow between internal functions like ts2clm, make_whole_fast, clim_spread, and clim_calc_cpp.

```mermaid
flowchart TD
    A(["ts2clm()"]) --> B[/INPUT:<br>date, temperature/]
    B --> C{Checks?}
    C --> |Yes| D[/OUTPUT:<br>date, temperature/]
    C --> |No| E([End])
    D --> F(["make_whole_fast()"])
    
    G(["make_whole_fast()"]) --> H[/INPUT:<br>date, temperature/]
    H --> I[/OUTPUT:<br>doy, date, temperature/]
    I --> J(["clim_spread()"])
    
    K(["clim_spread()"]) --> L[/INPUT:<br>doy, date, temperature/]
    L ---> M[spread doy as rows<br>spread year as cols<br>grow doy by windowHalfWidth]
    M --> N[/OUTPUT:<br>temperature in<br>doy x year matrix/]
    N --> O(["clim_calc_cpp()"])
    
    P(["clim_calc_cpp()"]) --> Q[/INPUT:<br>temperature in<br>doy x year matrix/]
    Q --> 
```

--------------------------------

### ts2clm(data, x, y, climatologyPeriod, ...)

Source: https://robwschlegel.github.io/heatwaveR/reference/ts2clm.html

Creates a daily climatology from a time series of daily temperatures using a user-specified sliding window for the mean and threshold calculation.

```APIDOC
## ts2clm

### Description
Creates a daily climatology from a time series of daily temperatures using a user-specified sliding window for the mean and threshold calculation, followed by an optional moving average smoother.

### Signature
`ts2clm(data, x = t, y = temp, climatologyPeriod, maxPadLength = FALSE, windowHalfWidth = 5, pctile = 90, smoothPercentile = TRUE, smoothPercentileWidth = 31, clmOnly = FALSE, var = FALSE, roundClm = 4, returnDF = TRUE)`

### Parameters
- **data** (data.frame) - Required - A data frame with two columns (default headers 't' and 'temp').
- **x** (column name) - Optional - Name of the column with dates.
- **y** (column name) - Optional - Name of the measurement variable column.
- **climatologyPeriod** (vector) - Required - Two values representing the start and end date of the climatology period.
- **maxPadLength** (integer/FALSE) - Optional - Maximum length of days to interpolate missing data.
- **windowHalfWidth** (integer) - Optional - Width of sliding window about day-of-year. Default is 5.
- **pctile** (numeric) - Optional - Threshold percentile for detection. Default is 90.
- **smoothPercentile** (boolean) - Optional - Whether to smooth the climatology and threshold. Default is TRUE.
- **smoothPercentileWidth** (integer) - Optional - Full width of moving average window. Default is 31.
- **clmOnly** (boolean) - Optional - If TRUE, returns only the climatologies. Default is FALSE.
- **var** (boolean) - Optional - If TRUE, calculates the variance of the seasonal signal. Default is FALSE.
- **roundClm** (numeric/FALSE) - Optional - Decimal places for rounding outputs. Default is 4.
- **returnDF** (boolean) - Optional - If TRUE, returns a data.frame; if FALSE, returns a data.table. Default is TRUE.
```

--------------------------------

### Download OISST data in batches

Source: https://robwschlegel.github.io/heatwaveR/articles/OISST_preparation.html

Uses tidyverse functions to download and process OISST data in batches, utilizing hoardr for local caching.

```R
# Download all of the data with one nested request
# The time this takes will vary greatly based on connection speed
base::system.time(
  OISST_data <- dl_years |> 
    dplyr::group_by(date_index) |> 
    dplyr::group_modify(~OISST_sub_dl(.x)) |> 
    dplyr::ungroup() |> 
    dplyr::select(lon, lat, t, temp)
) # 518 seconds, ~100 seconds per batch
```

--------------------------------

### Detect events with percentile thresholds

Source: https://robwschlegel.github.io/heatwaveR/reference/detect_event3.html

Configure event detection using two different percentile thresholds by extracting the climatology from the detect_event3 output.

```R
thresh_95 <- detect_event3(ts2clm3(sst_Med, pctile = 95,
                                   climatologyPeriod = c("1982-01-01", "2011-12-31")), 
                           minDuration = 2, maxGap = 0)$climatology
```

--------------------------------

### Define File Attributes and Dimensions

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_to_nc.html

Calculate spatial and event steps and define the dimensions required for the NetCDF structure.

```R
lon_step <- mean(diff(sort(unique(MHW_res_grid$lon))))
lat_step <- mean(diff(sort(unique(MHW_res_grid$lat))))
lon <- seq(min(MHW_res_grid$lon), max(MHW_res_grid$lon), by = lon_step)
lat <- seq(min(MHW_res_grid$lat), max(MHW_res_grid$lat), by = lat_step)
event_no <- seq(min(MHW_res_grid$event_no), max(MHW_res_grid$event_no), by = 1)
tunits <- "days since 1970-01-01"

# Length of each attribute
nlon <- length(lon)
nlat <- length(lat)
nen <- length(event_no)
```

--------------------------------

### Summarize MHW metrics by realm in R

Source: https://robwschlegel.github.io/heatwaveR/articles/MHW_metric_trends.html

Joins pixel counts, calculates mean metrics per ecoregion, and fills missing years with zero values.

```R
MHW_summary_NZ_realms <- left_join(MHW_summary_NZ, MHW_summary_NZ_npix_realms, by = 'ECOREGION') %>% 
  pivot_wider(names_from = metrics,values_from = values) %>% 
  group_by(year, ECOREGION, season, npix) %>% 
  summarize(Number_MHW_days = sum(nMHWdays),
            Nevents = sum(nevents), 
            Mean_Intensity = mean(int_mean),
            Maximum_Intensity = mean(int_max),
            Cumulative_Intensity = sum(int_cumulative), .groups = "drop") %>% 
  mutate(Number_MHW_days = Number_MHW_days/npix,
         Nevents = Nevents/npix, 
         Cumulative_Intensity = Cumulative_Intensity/npix) %>% 
  group_by(ECOREGION, season) %>% 
  complete(year = 1982:2021) %>% 
  dplyr::select(-npix) %>% 
  pivot_longer(-c(year, ECOREGION, season), names_to = 'Metrics', values_to = 'values') %>% 
  mutate(Metrics = factor(Metrics, levels = c("Number_MHW_days", "Nevents", "Mean_Intensity", 
                                              "Maximum_Intensity", "Cumulative_Intensity"))) %>% 
  replace(is.na(.), 0)
```

--------------------------------

### detect_event3(data, x, y, seasClim, threshClim, ...)

Source: https://robwschlegel.github.io/heatwaveR/reference/detect_event3.html

Applies the Hobday et al. (2016) marine heat wave definition to an input time series using pre-calculated climatologies.

```APIDOC
## detect_event3

### Description
Identifies marine heatwaves and cold-spells based on an input time series, date vector, and pre-calculated seasonal and threshold climatologies.

### Usage
```r
detect_event3(data, x = t, y = temp, seasClim = seas, threshClim = thresh, threshClim2 = NA, minDuration = 5, minDuration2 = minDuration, joinAcrossGaps = TRUE, maxGap = 2, maxGap2 = maxGap, coldSpells = FALSE, protoEvents = FALSE, categories = FALSE, roundRes = 4, ...)
```

### Parameters
- **data** (data.table) - Required - The input time series data.
- **x** (column) - Optional - The date vector column (default: t).
- **y** (column) - Optional - The value vector column (default: temp).
- **seasClim** (column) - Optional - The seasonal climatology column (default: seas).
- **threshClim** (column) - Optional - The threshold climatology column (default: thresh).
- **threshClim2** (column) - Optional - The secondary threshold climatology column (default: NA).
- **minDuration** (numeric) - Optional - Minimum duration for heatwaves (default: 5).
- **minDuration2** (numeric) - Optional - Minimum duration for cold-spells (default: minDuration).
- **joinAcrossGaps** (boolean) - Optional - Whether to join events across gaps (default: TRUE).
- **maxGap** (numeric) - Optional - Maximum gap size for heatwaves (default: 2).
- **maxGap2** (numeric) - Optional - Maximum gap size for cold-spells (default: maxGap).
- **coldSpells** (boolean) - Optional - Whether to detect cold-spells (default: FALSE).
- **protoEvents** (boolean) - Optional - Whether to return proto-events (default: FALSE).
- **categories** (boolean) - Optional - Whether to calculate event categories (default: FALSE).
- **roundRes** (numeric) - Optional - Rounding resolution (default: 4).
```

--------------------------------

### Calculate and Inspect Threshold Exceedances

Source: https://robwschlegel.github.io/heatwaveR/reference/exceedance.html

Calculates exceedances for a given dataset and threshold, then displays the resulting daily data and exceedance metrics.

```R
res <- exceedance(sst_WA, threshold = 25)
# show first ten days of daily data:
res$threshold[1:10, ]
#>             t  temp thresh threshCriterion durationCriterion exceedance
#> 1  1982-01-01 20.94     25           FALSE             FALSE      FALSE
#> 2  1982-01-02 21.25     25           FALSE             FALSE      FALSE
#> 3  1982-01-03 21.38     25           FALSE             FALSE      FALSE
#> 4  1982-01-04 21.16     25           FALSE             FALSE      FALSE
#> 5  1982-01-05 21.26     25           FALSE             FALSE      FALSE
#> 6  1982-01-06 21.61     25           FALSE             FALSE      FALSE
#> 7  1982-01-07 21.74     25           FALSE             FALSE      FALSE
#> 8  1982-01-08 21.50     25           FALSE             FALSE      FALSE
#> 9  1982-01-09 21.40     25           FALSE             FALSE      FALSE
#> 10 1982-01-10 21.36     25           FALSE             FALSE      FALSE
#>    exceedance_no
#> 1             NA
#> 2             NA
#> 3             NA
#> 4             NA
#> 5             NA
#> 6             NA
#> 7             NA
#> 8             NA
#> 9             NA
#> 10            NA
# show first five exceedances:
res$exceedance[1:5, ]
#>   exceedance_no index_start index_peak index_end duration date_start  date_peak
#> 1             1        2682       2683      2686        5 1989-05-05 1989-05-06
#> 2             2        6342       6351      6358       17 1999-05-13 1999-05-22
#> 3             3        6362       6363      6368        7 1999-06-02 1999-06-03
#> 4             4        6686       6688      6691        6 2000-04-21 2000-04-23
#> 5             5        6698       6699      6707       10 2000-05-03 2000-05-04
#>     date_end intensity_mean intensity_max intensity_var intensity_cumulative
#> 1 1989-05-09         0.2860          0.36        0.0462                 1.43
#> 2 1999-05-29         0.8559          1.40        0.3859                14.55
#> 3 1999-06-08         0.2071          0.27        0.0512                 1.45
#> 4 2000-04-26         0.1550          0.41        0.1803                 0.93
#> 5 2000-05-12         0.6970          1.01        0.2747                 6.97
#>   intensity_mean_abs intensity_max_abs intensity_var_abs intensity_cum_abs
#> 1            25.2860             25.36            0.0462            126.43
#> 2            25.8559             26.40            0.3859            439.55
#> 3            25.2071             25.27            0.0512            176.45
#> 4            25.1550             25.41            0.1803            150.93
#> 5            25.6970             26.01            0.2747            256.97
#>   rate_onset rate_decline
#> 1     0.1900       0.1171
#> 2     0.1826       0.1787
#> 3     0.1400       0.0555
#> 4     0.2540       0.1214
#> 5     0.5600       0.1371
```

--------------------------------

### Visualizing heatwave results

Source: https://robwschlegel.github.io/heatwaveR/articles/complex_clims.html

Create a bubble plot and multi-panel visualization for heatwave event metrics using ggplot2 and ggarrange.

```R
# The code to create a bubble plot for the heatwave results
bubble_plot <- ggplot(data = events$event, aes(x = date_peak, y = intensity_max)) +
  geom_point(aes(size = intensity_cumulative), shape = 21, fill = "salmon", alpha = 0.8) +
  labs(x = NULL, y = "Maximum Intensity [°C] ", size = "Cumulative Intensity [°C x days]") +
  scale_size_continuous(range = c(1, 10), 
                        guide = guide_legend(title.position = "top", direction = "horizontal")) +
  theme_bw() +
  theme(legend.position = c(0.3, 0.12),
        legend.box.background = element_rect(colour = "black"))

# Don't forget to set 'event_line(y = tMax)'
ggarrange(event_line(events, y = tMax, metric = intensity_max),
          event_line(events, y = tMax, metric = intensity_max, category = T),
          lolli_plot(events),
          bubble_plot,
          ncol = 2, nrow = 2, align = "hv")
```

--------------------------------

### Load Algiers dataset

Source: https://robwschlegel.github.io/heatwaveR/articles/complex_clims.html

Explicitly assigns the built-in Algiers dataset to the environment for clarity.

```R
Algiers <- Algiers
```
Tokens:
Raw