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