Detect Marine Heat Waves and Marine Cold Spells

Given a time series of daily temperatures, the package provides tools to detect extreme thermal events, including marine heat waves, and to calculate the exceedances above or below specified threshold values. It outputs the properties of all detected events and exceedances.

CRAN_Status_Badge Travis-CI Build Status

The RmarineHeatWaves package is a translation of the original Python code written by Eric C. J. Oliver that can be found on GitHub.

The RmarineHeatWaves R package contains a number of functions which calculate and display marine heat waves according to the definition of Hobday et al. (2016). The marine cold spell option was implemented in version 0.13 (21 Nov 2015) of the Python module as a result of the preparation of Schlegel et al. (2017), wherein the cold events are introduced and briefly discussed.

This package may be found on CRAN. Alternatively, you can install it from GitHub by issuing the following command:


The functions

Function Description
detect() The main function which detects the events as per the definition of Hobday et al. (2016).
make_whole() Constructs a continuous, uninterrupted time series of temperatures.
block_average() Calculates annual means for event metrics.
event_line() Creates a line plot of marine heat waves or cold spells.
lolli_plot() Creates a timeline of selected event metrics.
exceedance() A function similar to detect() but that detects consecutive days above/below a given threshold.
geom_flame() Creates flame polygons of marine heat waves or cold spells.
geom_lolli() Creates a lolliplot timeline of selected event metric.

The package also provides data of observed SST records for three historical MHWs: the 2011 Western Australia event, the 2012 Northwest Atlantic event and the 2003 Mediterranean event.

The detection and graphing functions

The detect() function is the package's core function. Here is the detect() function applied to the Western Australian test data, which are also discussed by Hobday et al. (2016):

library(RmarineHeatWaves); library(plyr); library(dplyr); library(ggplot2)
ts <- make_whole(sst_WA)
mhw <- detect(ts, climatology_start = 1983, climatology_end = 2012)
mhw$event %>% 
  dplyr::ungroup() %>%
  dplyr::select(event_no, duration, date_start, date_peak, int_mean, int_max, int_cum) %>% 
  dplyr::arrange(-int_cum) %>% 
#> 1       22       95 1999-05-13 1999-05-22 2.498305 3.601700 237.33900
#> 2       42       60 2011-02-06 2011-02-28 3.211903 6.505969 192.71420
#> 3       49       47 2012-01-11 2012-01-27 2.225734 3.300112 104.60948
#> 4       50       46 2012-03-01 2012-04-10 1.993709 2.957609  91.71061
#> 5       41       40 2010-12-24 2011-01-28 2.157016 3.274803  86.28064

The corresponding event_line() and lolli_plot(), which represent the massive Western Australian heatwave of 2011, look like this:

event_line(mhw, spread = 200, metric = "int_cum",
           start_date = "2010-10-01", end_date = "2011-08-30")


The event_line() and lolli_plot() functions were designed to work directly on one of the list returned by detect(). If more control over the figures is required, it may be useful to create them in ggplot2 by stacking 'geoms'. We specifically created two new ggplot2 geoms to reproduce the functionality of event_line() and lolli_plot(). These functions are more general in their functionality and can be used outside of the RmarineHeatWave package too. To apply them to MHWs and MCSs, they require that we access the clim or event data frames within the list that is produced by detect(). Here is how:

mhw2 <- mhw$clim # find the climatology dataframe
mhw2 <- mhw2[10580:10690,] # identify the region of the time series of interest
ggplot(mhw2, aes(x = t, y = temp, y2 = thresh_clim_year)) +
  geom_flame() +
  geom_text(aes(x = as.Date("2011-02-01"), y = 28, label = "The MHW that launched\na thousand papers."))

ggplot(mhw$event, aes(x = date_start, y = int_max)) +
  geom_lolli(colour = "salmon", colour.n = "red", n = 3) +
  geom_text(aes(x = as.Date("2006-10-01"), y = 5,
                label = "The distribution of events\nis skewed towards the\nend of the time series."),
            colour = "black")

The default output of these function may not be to your liking. If so, not to worry. As ggplot2 geoms, they are highly maleable. For example, if we were to choose to reproduce the format of the MHWs as seen in Hobday et al. (2016), the code would look something like this:

# 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[49:110,]
ggplot(data = mhw2, aes(x = t)) +
  geom_flame(aes(y = temp, y2 = thresh_clim_year, fill = "all"), show.legend = T) +
  geom_flame(data = mhw_top, aes(y = temp, y2 = thresh_clim_year, fill = "top"), show.legend = T) +
  geom_line(aes(y = temp, colour = "temp")) +
  geom_line(aes(y = thresh_clim_year, colour = "thresh"), size = 1.0) +
  geom_line(aes(y = seas_clim_year, 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")) +
  guides(colour = guide_legend(override.aes = list(fill = NA))) +
  xlab("Date") +
  ylab(expression(paste("Temperature [", degree, "C]")))

Conversely, should we not wish to highlight any events with geom_lolli(), it would look like this:

# Note that this is accomplished by setting 'colour.n = NA', not by setting 'n = 0'.
ggplot(mhw$event, aes(x = date_start, y = int_cum)) +
  geom_lolli(colour = "salmon", n = 3, colour.n = NA)

The calculation and visualisation of marine cold spells is also accommodated within this package. Here is a cold spell detected in the OISST data for Western Australia:

mcs <- detect(ts, climatology_start = 1983, climatology_end = 2012, cold_spells = TRUE)
mcs$event %>% 
  dplyr::ungroup() %>%
  dplyr::select(event_no, duration, date_start, date_peak, int_mean, int_max, int_cum) %>% 
  dplyr::arrange(int_cum) %>% 
#>   event_no duration date_start  date_peak  int_mean   int_max    int_cum
#> 1       16       76 1990-04-13 1990-05-11 -2.538017 -3.218054 -192.88929
#> 2       54       58 2003-12-19 2004-01-23 -1.798455 -2.662320 -104.31038
#> 3       71       52 2014-04-14 2014-05-05 -1.818984 -2.565533  -94.58715
#> 4        8       38 1986-06-24 1986-07-17 -2.009802 -2.950536  -76.37248
#> 5       51       32 2003-09-08 2003-09-16 -1.560817 -2.116583  -49.94613

The plots showing the marine cold spells look like this:

# this function needs to be updated as it currently only works when 't' is the column name for date
event_line(mcs, spread = 200, metric = "int_cum",
           start_date = "1990-01-01", end_date = "1990-08-30")

Cold spell figures may be created as geoms in ggplot2, too:

mcs2 <- mcs$clim
mcs2 <- mcs2[2990:3190,]
# # Note that the plot centres on the polygons, so it may be necessary to manually zoom out a bit
ggplot(data = mcs2, aes(x = t)) +
  geom_flame(aes(y = thresh_clim_year, y2 = temp), fill = "steelblue3", show.legend = F) +
  geom_line(aes(y = temp, colour = "temp")) +
  geom_line(aes(y = thresh_clim_year, colour = "thresh"), size = 1.0) +
  geom_line(aes(y = seas_clim_year, colour = "seas"), size = 1.2) +
  scale_colour_manual(name = "Line Colour",
                      values = c("temp" = "black", "thresh" =  "forestgreen", "seas" = "grey80")) +
  scale_y_continuous(limits = c(18, 23.5)) +
  xlab("Date") +
  ylab(expression(paste("Temperature [", degree, "C]")))
ggplot(mcs$event, aes(x = date_start, y = int_cum)) +
  geom_lolli(colour = "steelblue3", colour.n = "navy", n = 7) +
  xlab("Date") +
  ylab(expression(paste("Cumulative intensity [days x ", degree, "C]")))

The exceedance function

In addition to the calculation of extreme events, consecutive days over a given static threshold may be calculated with the exceedance() function.

exc <- exceedance(ts, threshold = 25)
exc$exceedance %>%
  ungroup() %>%
  select(exceedance_no, duration, date_start, date_peak, int_mean, int_cum) %>%
  dplyr::arrange(-int_cum) %>%
#>   exceedance_no duration date_start  date_peak  int_mean   int_cum
#> 1             7       52 2011-02-08 2011-02-28 1.6740379 87.049969
#> 2             6       25 2008-04-03 2008-04-14 0.9799994 24.499985
#> 3            10       41 2012-03-03 2012-04-10 0.4385360 17.979977
#> 4             2       17 1999-05-13 1999-05-22 0.8558818 14.549990
#> 5             5       10 2000-05-03 2000-05-04 0.6969994  6.969994

The same function may be used to calculate consecutive days below a threshold, too.

exc <- exceedance(ts, threshold = 19, below = TRUE)
exc$exceedance %>%
  dplyr::ungroup() %>%
  dplyr::select(exceedance_no, duration, date_start, date_peak, int_mean, int_cum) %>%
  dplyr::arrange(int_cum) %>%
#>   exceedance_no duration date_start  date_peak   int_mean   int_cum
#> 1            17       46 2003-09-06 2003-09-16 -0.6008700 -27.64002
#> 2            16       31 2002-09-08 2002-09-25 -0.8480649 -26.29001
#> 3            13       24 1997-09-03 1997-09-15 -0.7691671 -18.46001
#> 4            20       25 2005-09-26 2005-10-12 -0.5420004 -13.55001
#> 5            12       18 1997-08-13 1997-08-22 -0.6944449 -12.50001

Working with gridded SST data

We can also load the gridded 0.25 degree Reynolds OISST data and apply the function pixel by pixel over all of the days of data. The example data used here have 93 longitude steps, 43 latitude steps, and cover 12797 days (1981 to 2016). We apply the detect() function to these data, fit a generalised linear model (GLM), and then plot the trend per decade of the marine heatwave count. In other words, have marine heatwaves become more or less frequent in recent years? Under climate change we can expect that extreme events would tend to occur more frequently and be of greater intensity. Indeed, we can clearly see in the figure below of the result of the GLM, how the Agulhas Current has been experiencing marine heat waves more frequently in recent decades. But there are two smaller areas, one along the western side of the Cape Peninsula in the Benguela Upwelling system and another around the Eastern Cape Province near Algoa Bay, where the frequency of marine heat waves seems to have actually been decreasing -- although the P-value of the decreasing trend is > 0.05, and therefore not significant.



Please read the package vignette to see how to load a netCDF file with the OISST data, apply the RmarineHeatWaves function to the whole 3D array of data, and then fit the GLM and plot the data.


Hobday, A.J. et al. (2016). A hierarchical approach to defining marine heatwaves, Progress in Oceanography, 141, pp. 227-238.

Schlegel, R. W., Oliver, E. C. J., Wernberg, T. W., Smit, A. J. (2017). Coastal and offshore co-occurrences of marine heatwaves and cold-spells. Progress in Oceanography, 151, pp. 189-205.


The Python code was written by Eric C. J. Oliver.

Contributors to the Marine Heatwaves definition and its numerical implementation include Alistair J. Hobday, Lisa V. Alexander, Sarah E. Perkins, Dan A. Smale, Sandra C. Straub, Jessica Benthuysen, Michael T. Burrows, Markus G. Donat, Ming Feng, Neil J. Holbrook, Pippa J. Moore, Hillary A. Scannell, Alex Sen Gupta, and Thomas Wernberg.

The translation from Python to R was done by A. J. Smit and the graphing functions were contributed to by Robert. W. Schlegel.


A. J. Smit Department for Biodiversity & Conservation Biology, University of the Western Cape, Private Bag X17, Bellville 7535, South Africa, E-mail: [email protected], Work tel.: +27 (0)21 959 3783


output: pdf_document


Changes in version 0.16.0 (6 January 2018)

  • The make_whole(), detect() and exceedance() functions now handle any arbitrary column names, which are specified via arguments to the function.
  • Remove dependence on magrittr.
  • More consistent use of tidyverse packages throughout.
  • Updates to examples.
  • Miscellaneous other improvements troughout the function internals.

Changes in version 0.15.9 (19 June 2017)

  • Warnings given by exceedance() when exceedances begin/end on the first/last day of the time series have been corrected --- it will now report an NA when the exceedance occurs at the start/end of the time series.
  • Some plyr routines phased out in exceedance function in favour of a dplyr approach.

Changes in version 0.15.8 (18 June 2017)

  • Further removed need for plyr functions --- use dplyr instead.

Changes in version 0.15.7 (18 June 2017)

  • Warnings given by detect() when events begin/end on the first/last day of the time series have been corrected --- it will now report an NA when the event occurs at the start/end of the time series. This bug was highlighted by Juan Bazo.
  • Documentation in detect() about the use of percentiles with cold-spells has been made more clear.
  • Dependency on some plyr functions have been started to be phased out; this is done to eventually ensure full compliance with the tidyverse suite of functions.
  • Remove some unwanted attributes that resulted from dplyr::group_by on lines 496 and 497. The list containing 'event' and 'clim' is now uncluttered by unwanted attributes.

Changes in version 0.15.6 (6 March 2017)

  • Update README.Rmd to correct figure paths.

Changes in version 0.15.5 (6 March 2017)

  • Update to 'protoFunc' helper function within 'detect'; in the previous version it returned the error "Error: not compatible with STRSXP" under rare circumstances, the cause of which remains unknown.
  • Update the calculation of proto_gaps for the same reason.

Changes in version 0.15.4 (5 March 2017)

  • The detect() function now also produces a 'climatological' variance.
  • Improve documentation in a help file.

Changes in version 0.15.3 (28 February 2017)

  • The detect() function may now be told to calculate climatologies only.

Changes in version 0.15.2 (18 February 2017)

  • Updating, README.Rmd files.
  • Updating the vignette.

Changes in version 0.15.1 (15 February 2017)

  • 'exceedence' changed to 'exceedance' throughout.

Changes in version 0.15.0 (15 February 2017)

  • Final edits to 'geom_flame' and 'geom_lolli';
    • specifically in geom_flame, small gaps in the corners of polygons are now being filled in correctly.
  • Updated examples.

Changes in version 0.14.6 (15 February 2017)

  • Edits to 'geom_flame' and 'geom_lolli'.

Changes in version 0.14.5 (14 February 2017)

  • Edits to 'geom_flame' and 'geom_lolli'.

Changes in version 0.14.4 (11 February 2017)

  • Changed 'geom_event_line' and 'geom_lolli_plot' to 'geom_flame' and 'geom_lolli'
  • These new functions are now proper ggplot geoms.

Changes in version 0.14.3 (31 January 2017)

  • Corrected plotting error in event polygons in 'event_line' function
  • Added two new functions: 'geom_event_line' and 'geom_lolli_plot' which may be used with ggplot2 directly as geoms.

Changes in version 0.14.2 (26 January 2017)

  • Fix bug which was that was brought to my attention by Mahmoud Haouari: enable the pctile option in the detect() function, which was accidentally disabled.

Changes in version 0.14.1 (10 January 2017)

  • Improved error messages in exceedance().
  • Bug fix.

Changes in version 0.14.0 (11 December 2016)

  • All changes since version 0.13.1.

Changes in version (11 December 2016)

  • Minor edits to documentation.
  • Corrected issue where the event_line graph's y-axis did not reflect the correct units for the metric selected for plotting.
  • Updated the theme, i.e. removed the default ggplot2 theme in favour of a more conventional white background with nice black axes.
  • Removed the legend that indicated the peak and secondary events from the event_line plot.
  • Repositioned the legend in the event_line plot: it is now int he bottom right corner as here it has less chance of plotting over the peak events when the graph is rescaled.

Changes in version (23 November 2016)

  • Minor update to exceedance() function to improve usability.
  • Usage of exceedance() shown in README file

Changes in version (22 October 2016)

  • The additions to the detect() function madde in the previous update have been rolled back in favour of creating an independent function to calculate exceedances.

Changes in version (22 October 2016)

  • Added 'threshold' variable to detect() function.
  • This now allows the function to detect when temperatures are over (under) a static threshold supplied by the user.
  • An example is given in the file

Changes in version (8 October 2016)

  • Expanded example included as package vignette.

Changes in version (8 October 2016)

  • Added a marine cold spell example.

Changes in version (9 September 2016)

  • Correcting a few typos.

Changes in version (9 June 2016)

  • The new development version.

Changes in version 0.13.1 (9 June 2016)

  • Edit to file.

Changes in version 0.13.0 (8 June 2016)

  • Significant changes to detection algorithm:
    • ... more extensive use of magrittr pipe operators %>% and %<>%
    • ... simplification of protoFunc helper function
    • ... removal of automatic start/end day calculation if these are not provided
    • ... climatology start/end dates must now always be explicitely set
    • ... finally enabled join_across_gaps option
    • ... fixed a bug that caused the function to fail when none of the gaps between events (i.e. proto_gaps) were less than max_gap
  • Update examples to accommodate these changes.

Changes in version 0.12.2 (6 June 2016)

  • Add file.
  • Various minor changes to documentation of several functions.

Changes in version 0.12.1 (6 June 2016)

  • Fix Windows build dependency warnings.

Changes in version 0.12.0 (2 June 2016)

  • Added lolliplot functionality -- lolli_plot().

Changes in version 0.11.2 (2 June 2016)

  • Simplify make_whole() -- it should accept dates as class POSIXct or Date without the need for unneccesary 'if' logic.
  • All example data dates (t) changed to class Date.
  • Some minor rewording to documentation.

Changes in version 0.11.1 (1 June 2016)

  • Minor edits to event_line() as per Robert Schlegel (allows broader selection of metrics for plotting.)

Changes in version 0.11.0 (1 June 2016)

  • block_average() rewritten -- it is now based upon dplyr functions so it is faster and more stream-lined.
  • Completely removed the use of reshape2 in favour of tidyr.

Changes in version 0.10.3 (1 June 2016)

  • All comments removed inside of functions.

Changes in version 0.10.2 (1 June 2016)

  • Minor changes to the version change info provided for v0.10.1.

Changes in version 0.10.1 (31 May 2016)

  • Fixed a bug that caused detect() to fail whenever it encountered fewer than two non-NAs in the period doy 59 to doy 61 when it was asked to interpolate over the non-existent day-60 during non-leap years. This is specific to versions of 'zoo' (required for na.approx()) up to 1.7-12; from 1.7-13 it works fine. A few extra lines of code were added to fix this.

Changes in version 0.10.0 (30 May 2016)

  • Added the block_average() function.

Changes in version 0.9.3 (29 May 2016)

  • Expanded documentation.

Changes in version 0.9.2 (29 May 2016)

  • Changes to event_line() documentation.
  • Renamed 'metric' options to function as 'mean', 'maximum', 'cumulative'.
  • Added more TODOs at the bottom of this file.

Changes in version 0.9.1 (28 May 2016)

  • Minor refinements to the documentation (i.e. package description).

Changes in version 0.9.0 (28 May 2016)

  • Replace 'mhw' with 'event'.
  • Add basic plotting functionality a-la Robert Schlegel in the form of the event_line() function (plus edits to make Rob's code produce a clean build process, thereby avoiding throwing 'notes' that might be frowned upon by the CRAN people.)
  • Renamed some things: 'make_whole.R' -> 'makeWhole.R'; 'marineHeatWaves-package.r' -> 'RmarineHeatWaves-package.r'; 'marineHeatWaves.R' -> 'RmarineHeatWaves.R'.

Changes in version 0.8.2 (22 May 2016)

  • Include CITATION file.
  • Update DESCRIPTION file.
  • Update author, creator and contributor roles.

Changes in version 0.8.1 (22 May 2016)

  • Add dates to file.

Changes in version 0.8.0 (22 May 2016)

  • Marine cold spell option enabled.
  • Hacky changes to suppress some notes during package compilation, needed for acceptance to CRAN, apparently.

Changes in version 0.7.0 (20 May 2016)

  • Internal variable name changes (replace more camelCase).
  • Less problematic handling of NA.
  • Improved/expanded documentation.

Changes in version 0.6.0 (19 May 2016)

  • Internal variable name changes (replace camelCase).
  • Split out make_whole() function.
  • Add make_whole() example.

Changes in version 0.5.0 (18 May 2016)

  • Added data and a functional example.
  • Slightly improved documentation.

Changes in version 0.4.0 (18 May 2016)

  • Robert Schlegel added as co-author for work on graphing functions.

Changes in version 0.3.0 (18 May 2016)

  • Changed 'eventNo' calculation to properly reflect the actual events, and not one of the proto-event types.
  • Changed ID variable in 'mhw' output to 'eventNo', thereby following the same naming convention.
  • Removed ', .id = NULL' from the event metrics 'ldply' calculations, enabling compatibility with plyr (version < 1.8.3).

Changes in version 0.2.0 (17 May 2016)

  • Add 'eventNo' to climatology output to identify each unique event and to fascilitate plotting of filled polygons using 'geom_polygon' in ggplot2, as per Robert Schlegel's suggestion.

Reference manual

It appears you don't have a PDF plugin for this browser. You can click here to download the reference manual.


0.16.0 by Albertus J. Smit, 2 months ago

Browse source code at

Authors: Albertus J. Smit [aut, cre] (R implementation.), Eric C. J. Oliver [aut] (The brain behind the Python implementation.), Robert W. Schlegel [ctb] (Graphical and data summaries.)

Documentation:   PDF Manual  

MIT + file LICENSE license

Imports tibble, ggplot2, lubridate, dplyr, stats, utils, zoo, tidyr, plyr, raster, grid, lazyeval

Suggests knitr, rmarkdown

See at CRAN