Title: | Detect Marine Heat Waves and Marine Cold Spells |
---|---|
Description: | 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. |
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.) |
Maintainer: | Albertus J. Smit <[email protected]> |
License: | MIT + file LICENSE |
Version: | 0.17.1 |
Built: | 2025-02-18 04:57:00 UTC |
Source: | https://github.com/ajsmit/rmarineheatwaves |
Calculate Yearly Means for Event Metrics.
block_average(data, x = t, y = temp, report = "full")
block_average(data, x = t, y = temp, report = "full")
data |
Accepts the data returned by the |
x |
This column is expected to contain a vector of dates as per the
specification of |
y |
This is a column containing the measurement variable. If the column
name differs from the default (i.e. |
report |
Specify either |
This function needs to be provided with the full output from the detect
function. Note that the yearly averages are calculted only for complete years
(i.e. years that start/end part-way through the year at the beginning or end
of the original time series are removed from the calculations).
This function differs from the python implementation of the function of the
same name (i.e., blockAverage
, see https://github.com/ecjoliver/marineHeatWaves)
in that we only provide the ability to calculate the average (or aggregate)
event metrics in 'blocks' of one year, while the python version allows
arbitrary (integer) block sizes.
The function will return a data frame of the averaged (or aggregate) metrics. It includes the following:
year |
The year over which the metrics were averaged. |
temp_mean |
Seawater temperature for the specified year [deg. C]. |
temp_min |
The minimum temperature for the specified year [deg. C]. |
temp_max |
The maximum temperature for the specified year [deg. C]. |
count |
The number of events per year. |
duration |
The average duration of events per year [days]. |
int_mean |
The average event "mean intensity" in each year [deg. C]. |
int_max |
The average event "maximum (peak) intensity" in each year [deg. C]. |
int_var |
The average event "intensity variability" in each year [deg. C]. |
int_cum |
The average event "cumulative intensity" in each year [deg. C x days]. |
rate_onset |
Average event onset rate in each year [deg. C / days]. |
rate_decline |
Average event decline rate in each year [deg. C / days]. |
total_days |
Total number of events days in each year [days]. |
total_icum |
Total cumulative intensity over all events in each year [deg. C x days]. |
int_max_rel_thresh
, int_mean_rel_thresh
,
int_var_rel_thresh
, and int_cum_rel_thresh
are as above except relative to the threshold (e.g., 90th percentile) rather
than the seasonal climatology.
int_max_abs
, int_mean_abs
, int_var_abs
, and
int_cum_abs
are as above except as absolute magnitudes
rather than relative to the seasonal climatology or threshold.
int_max_norm
and int_mean_norm
are as above except
units are in multiples of threshold exceedances, i.e., a value of 1.5
indicates the event intensity (relative to the climatology) was 1.5 times the
value of the threshold (relative to climatology,
i.e., threshold - climatology.)
Albertus J. Smit, Eric C. J. Oliver
Hobday, A.J. et al. (2016), A hierarchical approach to defining marine heatwaves, Progress in Oceanography, 141, pp. 227-238, doi: 10.1016/j.pocean.2015.12.014
# ts_dat <- make_whole(sst_Med) # res <- detect(ts_dat, climatology_start = "1983-01-01", # climatology_end = "2012-12-31") # out <- block_average(res) # summary(glm(count ~ year, out, family = "poisson")) ## Not run: plot(out$year, out$count, col = "salmon", pch = 16, xlab = "Year", ylab = "Number of events") lines(out$year, out$count) ## End(Not run)
# ts_dat <- make_whole(sst_Med) # res <- detect(ts_dat, climatology_start = "1983-01-01", # climatology_end = "2012-12-31") # out <- block_average(res) # summary(glm(count ~ year, out, family = "poisson")) ## Not run: plot(out$year, out$count, col = "salmon", pch = 16, xlab = "Year", ylab = "Number of events") lines(out$year, out$count) ## End(Not run)
Applies the Hobday et al. (2016) marine heat wave definition to an input time series of temperature along with a daily date vector.
detect(data, doy = doy, x = t, y = temp, climatology_start, climatology_end, pctile = 90, window_half_width = 5, smooth_percentile = TRUE, smooth_percentile_width = 31, clim_only = FALSE, min_duration = 5, join_across_gaps = TRUE, max_gap = 2, max_pad_length = 3, cold_spells = FALSE)
detect(data, doy = doy, x = t, y = temp, climatology_start, climatology_end, pctile = 90, window_half_width = 5, smooth_percentile = TRUE, smooth_percentile_width = 31, clim_only = FALSE, min_duration = 5, join_across_gaps = TRUE, max_gap = 2, max_pad_length = 3, cold_spells = FALSE)
data |
A data frame with three columns. In the default setting (i.e. ommitting
the arguments |
doy |
If a column headed |
x |
This column is expected to contain a vector of dates as per the
specification of |
y |
This is a column containing the measurement variable. If the column
name differs from the default (i.e. |
climatology_start |
Required. The start date for the period across which the (varying by day-of-year) seasonal cycle and extremes threshold are calculated. |
climatology_end |
Required. The end date for the period across which the (varying by day-of-year) seasonal cycle and extremes threshold are calculated. |
pctile |
Threshold percentile (%) for detection of extreme values.
Default is |
window_half_width |
Width of sliding window about day-of-year (to one
side of the center day-of-year) used for the pooling of values and
calculation of climatology and threshold percentile. Default is |
smooth_percentile |
Boolean switch selecting whether to smooth the
climatology and threshold percentile timeseries with a moving average of
width |
smooth_percentile_width |
Full width of moving average window for smoothing
climatology and threshold. Default is |
clim_only |
Choose to calculate only the climatologies and not the
events. Default is |
min_duration |
Minimum duration for acceptance of detected MHWs.
Default is |
join_across_gaps |
Boolean switch indicating whether to join MHWs which
occur before/after a short gap as specified by |
max_gap |
Maximum length of gap allowed for the joining of MHWs. Default
is |
max_pad_length |
Specifies the maximum length of days over which to
interpolate (pad) missing data (specified as |
cold_spells |
Boolean specifying if the code should detect cold events
instead of heat events. Default is |
This function assumes that the input time series consists of continuous
daily values with few missing values. Time ranges which start and end
part-way through the calendar year are supported. The accompanying function
make_whole
aids in the preparation of a time series that is
suitable for use with detect
, although this may also be accomplished
'by hand' as long as the criteria are met as discussed in the documentation
to make_whole
.
It is recommended that a climatology period of at least 30 years is specified in order to capture decadal thermal periodicities. It is further advised that full the start and end dates for the climatology period result in full years, e.g. "1982-01-01" to "2011-12-31" or "1982-07-01" to "2012-06-30"; if not, this may result in an unequal weighting of data belonging with certain months within a time series.
This function supports leap years. This is done by ignoring Feb 29s for the initial calculation of the climatology and threshold. The values for Feb 29 are then linearly interpolated from the values for Feb 28 and Mar 1.
The calculation of onset and decline rates assumes that the events
started a half-day before the start day and ended a half-day after the
end-day. This is consistent with the duration definition as implemented,
which assumes duration = end day - start day + 1. As of version 0.15.7, an
event that is already present at the beginning of a time series, or an event
that is still present at the end of a time series, will report the rate of
onset or the rate of decline as NA
, as it is impossible to know what
the temperature half a day before or after the start or end of the event is.
This may be a departure from the python marineHeatWaves function.
For the purposes of event detection, any missing temperature values not
interpolated over (through optional max_pad_length
) will be set equal
to the seasonal climatology. This means they will trigger the end/start of
any adjacent temperature values which satisfy the event definition criteria.
If the code is used to detect cold events (coldSpells
= TRUE),
then it works just as for heat waves except that events are detected as
deviations below the (100 - pctile)th percentile (e.g., the 10th instead of
90th) for at least 5 days. Intensities are reported as negative values and
represent the temperature anomaly below climatology.
If only the climatology for the time series is required, and not the
events themselves, this may be done by setting clim_only
= TRUE.
The original Python algorithm was written by Eric Oliver, Institute for Marine and Antarctic Studies, University of Tasmania, Feb 2015, and is documented by 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 our preparation of Schlegel et al. (submitted), wherein the cold events receive a brief overview.
The function will return a list of two tibbles (see the tidyverse
),
clim
and event
, which are the climatology and events,
respectively. The climatology contains the full time series of daily temperatures,
as well as the the seasonal climatology, the threshold and various aspects of the
events that were detected. The software was designed for detecting extreme
thermal events, and the units specified below reflect that intended purpose.
However, the various other kinds of extreme events may be detected according
to the 'marine heat wave' specifications, and if that is the case, the appropriate
units need to be determined by the user.
doy |
Julian day (day-of-year). For non-leap years it runs 1...59 and
61...366, while leap years run 1...366. This column will be named differently if
another name was specified to the |
t |
The date of the temperature measurement. This column will be
named differently if another name was specified to the |
temp |
If the software was used for the purpose for which it was designed,
seawater temperature [deg. C] on the specified date will be returned. This
column will of course be named differently if another kind of measurement was
specified to the |
seas_clim_year |
Climatological seasonal cycle [deg. C]. |
thresh_clim_year |
Seasonally varying threshold (e.g., 90th percentile) [deg. C]. |
var_clim_year |
Seasonally varying variance (standard deviation) [deg. C]. |
thresh_criterion |
Boolean indicating if |
duration_criterion |
Boolean indicating whether periods of consecutive
|
event |
Boolean indicating if all criteria that define a MHW or MCS are met. |
event_no |
A sequential number indicating the ID and order of occurence of the MHWs or MCSs. |
The events are summarised using a range of event metrics:
index_start |
Start index of event. |
index_stop |
Stop index of event. |
event_no |
A sequential number indicating the ID and order of the events. |
duration |
Duration of event [days]. |
date_start |
Start date of event [date]. |
date_stop |
Stop date of event [date]. |
date_peak |
Date of event peak [date]. |
int_mean |
Mean intensity [deg. C]. |
int_max |
Maximum (peak) intensity [deg. C]. |
int_var |
Intensity variability (standard deviation) [deg. C]. |
int_cum |
Cumulative intensity [deg. C x days]. |
rate_onset |
Onset rate of event [deg. C / day]. |
rate_decline |
Decline rate of event [deg. C / day]. |
int_max_rel_thresh
, int_mean_rel_thresh
,
int_var_rel_thresh
, and int_cum_rel_thresh
are as above except relative to the threshold (e.g., 90th percentile) rather
than the seasonal climatology.
int_max_abs
, int_mean_abs
, int_var_abs
, and
int_cum_abs
are as above except as absolute magnitudes
rather than relative to the seasonal climatology or threshold.
int_max_norm
and int_mean_norm
are as above except
units are in multiples of threshold exceedances, i.e., a value of 1.5
indicates the event intensity (relative to the climatology) was 1.5 times the
value of the threshold (relative to climatology,
i.e., threshold - climatology.)
Note that rate_onset
and rate_decline
will return NA
when the event begins/ends on the first/last day of the time series. This
may be particularly evident when the function is applied to large gridded
data sets. Although the other metrics do not contain any errors and
provide sensible values, please take this into account in its
interpretation.
Albertus J. Smit, Robert W. Schlegel, Eric C. J. Oliver
Hobday, A.J. et al. (2016). A hierarchical approach to defining marine heatwaves, Progress in Oceanography, 141, pp. 227-238, doi:10.1016/j.pocean.2015.12.014
Schlegel, R. W., Oliver, 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, doi:10.1016/j.pocean.2017.01.004
ts_dat <- make_whole(sst_WA) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") # show a portion of the climatology: res$clim[1:10, ] # show some of the heat waves: res$event[1:5, 1:10]
ts_dat <- make_whole(sst_WA) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") # show a portion of the climatology: res$clim[1:10, ] # show some of the heat waves: res$event[1:5, 1:10]
Creates a graph of warm or cold events as per the second row of Figure 3 in Hobday et al. (2016).
event_line(data, x = t, y = temp, min_duration = 5, spread = 150, metric = "int_cum", start_date, end_date)
event_line(data, x = t, y = temp, min_duration = 5, spread = 150, metric = "int_cum", start_date, end_date)
data |
The function receives the output from the |
x |
This column is expected to contain a vector of dates as per the
specification of |
y |
This is a column containing the measurement variable. If the column
name differs from the default (i.e. |
min_duration |
The minimum duration that an event has to for it to qualify as a marine heat wave or marine cold spell. |
spread |
The the number of days leading and trailing the largest event
(as per |
metric |
One of the following options: |
start_date |
The start date of a period of time within which the largest
event (as per |
end_date |
The end date of a period of time within which the largest
event (as per |
The function will return a line plot indicating the climatology,
threshold and temperature, with the hot or cold events that meet the
specifications of Hobday et al. (2016) shaded in as appropriate. The plotting
of hot or cold events depends on which option is specified in detect
.
The top event detect during the selected time period will be visible in a
brighter colour. This function differs in use from geom_flame
in that it creates a stand alone figure. The benefit of this being
that one must not have any prior knowledge of ggplot2 to create the figure.
Robert W. Schlegel
Hobday, A.J. et al. (2016), A hierarchical approach to defining marine heatwaves, Progress in Oceanography, 141, pp. 227-238, doi: 10.1016/j.pocean.2015.12.014
ts_dat <- make_whole(sst_WA) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") ## Not run: event_line(res, spread = 200, metric = "int_cum", start_date = "2010-10-01", end_date = "2011-08-30") ## End(Not run)
ts_dat <- make_whole(sst_WA) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") ## Not run: event_line(res, spread = 200, metric = "int_cum", start_date = "2010-10-01", end_date = "2011-08-30") ## End(Not run)
Detect consecutive days in exceedance of a given threshold.
exceedance(data, x = t, y = temp, threshold = 20, below = FALSE, min_duration = 5, join_across_gaps = TRUE, max_gap = 2, max_pad_length = 3)
exceedance(data, x = t, y = temp, threshold = 20, below = FALSE, min_duration = 5, join_across_gaps = TRUE, max_gap = 2, max_pad_length = 3)
data |
A data frame with at least the two following columns:
a |
x |
This column is expected to contain a vector of dates as per the
specification of |
y |
This is a column containing the measurement variable. If the column
name differs from the default (i.e. |
threshold |
The static threshold used to determine how many consecutive
days are in exceedance of the temperature of interest. Default is
|
below |
Default is |
min_duration |
Minimum duration that temperatures must be in exceedance
of the |
join_across_gaps |
A TRUE/FALSE statement that indicates whether
or not to join consecutive days of temperatures in exceedance of the
|
max_gap |
The maximum length of the gap across which to connect
consecutive days in exceedance of the |
max_pad_length |
Specifies the maximum length of days over which to
interpolate (pad) missing data (specified as |
This function assumes that the input time series consists of continuous
daily temperatures, with few missing values. The accompanying function
make_whole
aids in the preparation of a time series that is
suitable for use with exceedance
, although this may also be accomplished
'by hand' as long as the criteria are met as discussed in the documentation
to make_whole
.
Future versions seek to accomodate monthly and annual time series, too.
The calculation of onset and decline rates assumes that exceedance of the
threshold
started a half-day before the start day and ended a half-day
after the end-day. This is consistent with the duration definition as implemented,
which assumes duration = end day - start day + 1.
For the purposes of exceedance detection, any missing temperature values not
interpolated over (through optional max_pad_length
) will remain as
NA
. This means they will trigger the end of an exceedance if the adjacent
temperature values are in exceedance of the threshold
.
If the function is used to detect consecutive days of temperature under
the given theshold
, these temperatures are then taken as being in
exceedance below the threshold
as there is no antonym in the English
language for 'exceedance'.
This function is based largely on the detect
function found in this
package, which was ported from the Python algorithm that was written by Eric
Oliver, Institute for Marine and Antarctic Studies, University of Tasmania,
Feb 2015, and is documented by Hobday et al. (2016).
The function will return a list of two components. The first being
threshold
, which shows the daily temperatures and on which specific days
the given threshold
was exceeded. The second component of the list is
exceedance
, which shows a medley of statistics for each discrete group
of days in exceedance of the given threshold
. Note that any additional
columns left in the data frame given to this function will be output in the
threshold
component of the output. For example, if one uses
make_whole
to prepare a time series for analysis and leaves
in the doy
column, this column will appear in the output.
The information shown in the threshold
component is:
t |
The date of the temperature measurement. This variable may named
differently if an alternative name is supplied to the function's |
temp |
Temperature on the specified date [deg. C]. This variable may
named differently if an alternative name is supplied to the function's |
thresh |
The static |
thresh_criterion |
Boolean indicating if |
duration_criterion |
Boolean indicating whether periods of consecutive
|
exceedance |
Boolean indicting if all criteria that define a discrete
group in exceedance of the |
exceedance_no |
A sequential number indicating the ID and order of occurence of exceedances. |
The individual exceedances are summarised using the following metrics:
index_start |
Row number on which exceedance starts. |
index_stop |
Row number on which exceedance stops. |
exceedance_no |
The same sequential number indicating the ID and
order of the exceedance as found in the |
duration |
Duration of exceedance [days]. |
date_start |
Start date of exceedance [date]. |
date_stop |
Stop date of exceedance [date]. |
date_peak |
Date of exceedance peak [date]. |
int_mean |
Mean intensity [deg. C]. |
int_max |
Maximum (peak) intensity [deg. C]. |
int_var |
Intensity variability (standard deviation) [deg. C]. |
int_cum |
Cumulative intensity [deg. C x days]. |
rate_onset |
Onset rate of exceedance [deg. C / day]. |
rate_decline |
Decline rate of exceedance [deg. C / day]. |
int_max_abs
, int_mean_abs
, int_var_abs
, and
int_cum_abs
are as above except as absolute magnitudes
rather than relative to the threshold.
Robert W. Schlegel, Albertus J. Smit
ts_dat <- make_whole(sst_WA) res <- exceedance(ts_dat, threshold = 25) # show first ten days of daily data: res$threshold[1:10, ] # show first five exceedances: res$exceedance[1:5, ]
ts_dat <- make_whole(sst_WA) res <- exceedance(ts_dat, threshold = 25) # show first ten days of daily data: res$threshold[1:10, ] # show first five exceedances: res$exceedance[1:5, ]
This function will create polygons between two lines. If given a
temperature and theshold time series, like that produced by detect
,
the output will meet the specifications of Hobday et al. (2016) shown as
'flame polygons.' If one wishes to plot polygons below a given threshold, and not
above, switch the values being fed to the y
and y2
aesthetics. This function differs in use from event_line
in that it must be created as a ggplot
'geom' object. The benefit
of this being that one may add additional information to the figure as geom
layers to ggplot2 graphs as may be necessary.
geom_flame(mapping = NULL, data = NULL, stat = "identity", position = "identity", ..., na.rm = FALSE, show.legend = NA, inherit.aes = TRUE)
geom_flame(mapping = NULL, data = NULL, stat = "identity", position = "identity", ..., na.rm = FALSE, show.legend = NA, inherit.aes = TRUE)
mapping |
Set of aesthetic mappings created by |
data |
The data to be displayed in this layer. There are three options: If NULL, the default, the data is inherited from the plot data as specified
in the call to A data.frame, or other object, will override the plot data. All objects will
be fortified to produce a data frame. See A function will be called with a single argument, the plot data. The return
value must be a |
stat |
The statistical transformation to use on the data for this layer, as a string. |
position |
Position adjustment, either as a string, or the result of a call to a position adjustment function. |
... |
other arguments passed on to |
na.rm |
If |
show.legend |
Logical. Should this layer be included in the legends? |
inherit.aes |
If |
geom_flame
understands the following aesthetics (required aesthetics
are in bold):
x
y
y2
colour
fill
size
alpha
linetype
Robert W. Schlegel
Hobday, A.J. et al. (2016), A hierarchical approach to defining marine heatwaves, Progress in Oceanography, 141, pp. 227-238, doi: 10.1016/j.pocean.2015.12.014
event_line
for a non-ggplot2 based flame function.
ts_dat <- make_whole(sst_WA) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") mhw <- res$clim mhw <- mhw[10580:10690,] ## Not run: require(ggplot2) ggplot(mhw, aes(x = t, y = temp)) + geom_flame(aes(y2 = thresh_clim_year)) + geom_text(aes(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]"))) ## End(Not run)
ts_dat <- make_whole(sst_WA) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") mhw <- res$clim mhw <- mhw[10580:10690,] ## Not run: require(ggplot2) ggplot(mhw, aes(x = t, y = temp)) + geom_flame(aes(y2 = thresh_clim_year)) + geom_text(aes(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]"))) ## End(Not run)
The function will return a graph of the intensity of the selected
metric along the *y*-axis versus a time variable along the *x*-axis.
The number of top events (n
) from the chosen metric may be highlighted
in a brighter colour with the aesthetic value colour.n
.
This function differs in use from lolli_plot
in that it must be created as a ggplot2 'geom' object. The benefit of this being
that one may add additional information layer by layer to the figure as
geoms as necessary.
geom_lolli(mapping = NULL, data = NULL, ..., n = 1, na.rm = FALSE, show.legend = NA, inherit.aes = TRUE)
geom_lolli(mapping = NULL, data = NULL, ..., n = 1, na.rm = FALSE, show.legend = NA, inherit.aes = TRUE)
mapping |
Set of aesthetic mappings created by |
data |
The data to be displayed in this layer. There are three options: If NULL, the default, the data is inherited from the plot data as specified
in the call to A data.frame, or other object, will override the plot data. All objects will
be fortified to produce a data frame. See A function will be called with a single argument, the plot data. The return
value must be a |
... |
other arguments passed on to |
n |
The number of top events to highlight. Default is 1. This parameter
has no effect if |
na.rm |
If |
show.legend |
Logical. Should this layer be included in the legends? |
inherit.aes |
If |
geom_lolli
understands the following aesthetics (required aesthetics
are in bold):
x
y
alpha
color
linetype
size
shape
stroke
fill
colour.n
: While this value may be used as an aesthetic, it also
works as a parameter for this function. If one chooses not to highlight
any events, use colour.n = NA
outside of aes()
. One may
also provide a non-static value to colour.na
but remember that
one may not provide multiple continuous or discrete scales to a single
ggplot2 object. Therefore, if one provides a continuous value to
aes(colour)
, the values supplied to colour.n
must be
discrete. ggplot2 will attempt to do this automatically.
Robert W. Schlegel
lolli_plot
for a non-geom based lolliplot function.
ts_dat <- make_whole(sst_NW_Atl) # with defaults: res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") mhw <- res$event ## Not run: require(lubridate) # Height of lollis represent event durations and their colours # are mapped to the events' cumulative intensity: ggplot(mhw, aes(x = mhw$date_peak, y = mhw$duration)) + geom_lolli(n = 0, shape = 20, aes(colour = mhw$int_cum), colour.n = NA) + 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 = mhw$date_peak, y = mhw$duration)) + geom_lolli(n = 3, shape = 20, colour.n = "red") + scale_color_distiller(palette = "Spectral", name = "Cumulative \nintensity") + xlab("Date") + ylab("Event duration [days]") ## End(Not run)
ts_dat <- make_whole(sst_NW_Atl) # with defaults: res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") mhw <- res$event ## Not run: require(lubridate) # Height of lollis represent event durations and their colours # are mapped to the events' cumulative intensity: ggplot(mhw, aes(x = mhw$date_peak, y = mhw$duration)) + geom_lolli(n = 0, shape = 20, aes(colour = mhw$int_cum), colour.n = NA) + 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 = mhw$date_peak, y = mhw$duration)) + geom_lolli(n = 3, shape = 20, colour.n = "red") + scale_color_distiller(palette = "Spectral", name = "Cumulative \nintensity") + xlab("Date") + ylab("Event duration [days]") ## End(Not run)
Visualise a timeline of several event metrics as 'lollipop' graphs.
lolli_plot(data, metric = "int_max", event_count = 3, xaxis = "date_start")
lolli_plot(data, metric = "int_max", event_count = 3, xaxis = "date_start")
data |
Output from the |
metric |
One of |
event_count |
The number of top events to highlight. Default is 3. |
xaxis |
One of |
The function will return a graph of the intensity of the selected
metric along the y-axis versus either t
or event_no
.
The number of top events as per event_count
will be highlighted
in a brighter colour. This function differs in use from geom_lolli
in that it creates a stand alone figure. The benefit of this being
that one must not have any prior knowledge of ggplot2 to create the figure.
Albertus J. Smit and Robert W. Schlegel
ts_dat <- make_whole(sst_NW_Atl) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") ## Not run: lolli_plot(res, metric = "int_cum", event_count = 3, xaxis = "date_peak") ## End(Not run)
ts_dat <- make_whole(sst_NW_Atl) res <- detect(ts_dat, climatology_start = "1983-01-01", climatology_end = "2012-12-31") ## Not run: lolli_plot(res, metric = "int_cum", event_count = 3, xaxis = "date_peak") ## End(Not run)
Takes a series of dates and temperatures, and if irregular (but ordered), inserts missing dates and fills correpsonding temperatures with NAs.
make_whole(data, x = t, y = temp)
make_whole(data, x = t, y = temp)
data |
A data frame with columns for date and temperature data. Ordered daily data are expected, and although missing values (NA) can be accommodated, the function is only recommended when NAs occur infrequently, preferably at no more than 3 consecutive days. |
x |
A column with the daily time vector (see details). For backwards
compatibility, the column is named |
y |
A column with the response vector. RmarineHeatWaves version <= 0.15.9
assumed that this would be daily seawater temperatures, but as of version 0.16.0
it may be any arbitrary measurement taken at a daily frequency. The default
remains temperature, and the default column name is therefore |
Upon import, the package uses 'zoo' and 'lubridate' to process the input
date and temperature data. It reads in daily data with the time vector
specified as either POSIXct
or Date
(e.g. "1982-01-01 02:00:00" or
"1982-01-01"). The data may be an irregular time series, but date must be
ordered. The function constructs a complete time series from the start date
to the end date, and fills in the regions in the time series where temperature
data are missing, with NAs in the temperature vector. There must only be one
temperature value per day otherwise the function will fail. It is up to the
user to calculate daily data from sub-daily measurements. Leap years are
automatically accommodated by 'zoo'.
This function can handle some of missing days, but this is not a
licence to actually use these data for the detection of anomalous thermal
events. Hobday et al. (2016) recommend gaps of no more than 3 days, which
may be adjusted by setting the max_pad_length
argument of the
detect
function. The longer and more frequent the gaps become
the lower the fidelity of the annual climatology and threshold that can be
calculated, which will not only have repercussions for the accuracy at which
the event metrics can be determined, but also for the number of events that
can be detected.
It is recommended that a climatology period of at least 30 years is specified in order to capture any decadal thermal periodicities.
The function will return a data frame with three columns. The column
headed doy
(day-of-year) is the Julian day running from 1 to 366, but
modified so that the day-of-year series for non-leap-years runs 1...59 and
then 61...366. For leap years the 60th day is February 29. See the example,
below. The other two columns take the names of x
and y
, if supplied,
or it will be t
and temp
in case the default values were used.
The x
(or t
) column is a series of dates of class Date
,
while y
(or temp
) is the measured variable. This time series will
be uninterrupted and continuous daily values between the first and last dates
of the input data.
Smit, A. J.
require(dplyr); require(tidyr); require(lubridate) ts_dat <- make_whole(sst_WA) # default columns "t" and "temp", in that order clim_start <- "1983-01-01" clim_end <- "2012-12-31" ts_dat %>% filter(t >= clim_start & t <= clim_end) %>% mutate(t = year(t)) %>% spread(t, temp) %>% filter(doy >= 55 & doy <= 65)
require(dplyr); require(tidyr); require(lubridate) ts_dat <- make_whole(sst_WA) # default columns "t" and "temp", in that order clim_start <- "1983-01-01" clim_end <- "2012-12-31" ts_dat %>% filter(t >= clim_start & t <= clim_end) %>% mutate(t = year(t)) %>% spread(t, temp) %>% filter(doy >= 55 & doy <= 65)
This package is an R implementation of the python script marineHeatWaves
(https://github.com/ecjoliver/marineHeatWaves) written by Eric C. J.
Oliver as part of the marine heat waves definition by Hobday et al. (2016).
Although the title of the package refers to marine heat waves (MHW), it is equally capable of detecting marine cold spells (MCS). This functionality to detect cold events is also present in the python package, where it was implemented as a result of the publication Schlegel et al. (2017) that discusses the quantification and detection of anomalously cold events. As of release 0.16.0, the detect function may also be applied to time series of other natural phenomena which one might want to express in terms of the summary metrics outlined in the paper by Hobday et al. (2016).
The main function is the detection function detect
which takes as
input a time series of temperature (and a corresponding series of dates)
and outputs a set of detected MHWs or MCS, as well as the climatological
(varying by day-of-year) seasonal cycle and extremes threshold. There are
various helper functions to fascilitate developing an uninterrupted time
series of temperatures (e.g. make_whole
) and some options to produce
graphical summaries and representations of the detected events such as
event_line
and lolli_plot
, or the ggplot2
equivalents, geom_flame
and geom_lolli
.
This package is demonstrated by applying the MHW definition to observed SST records and showing how it identifies three historical MHWs: the 2011 Western Australia event, the 2012 Northwest Atlantic event and the 2003 Mediterranean event. These data are included herewith.
One may also use the exceedance
function to calculate consecutive
days above or below a given static threshold. The output of this function is
similar to detect
.
Albertus J. Smit <[email protected]>, Robert W. Schlegel, Eric C. J. Oliver
Hobday, A. J. et al. (2016), A hierarchical approach to defining marine heatwaves. Progress in Oceanography, 141, pp. 227-238, <DOI:10.1016/j.pocean.2015.12.014> (official citation for this package).
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, <DOI:10.1016/j.pocean.2017.01.004>
A dataset containing the sea surface temperature (in degrees Celsius) and date for the Mediterranean region from 1982-01-01 to 2014-12-31.
sst_Med
sst_Med
A data frame with 12053 rows and 2 variables:
date, as POSIXct
SST, in degrees Celsius
...
https://www.ncdc.noaa.gov/oisst
A dataset containing the sea surface temperature (in degrees Celsius) and date for the Northwest Atlantic region from 1982-01-01 to 2014-12-31.
sst_NW_Atl
sst_NW_Atl
A data frame with 12053 rows and 2 variables:
date, as POSIXct
SST, in degrees Celsius
...
https://www.ncdc.noaa.gov/oisst
A dataset containing the sea surface temperature temperature (in degrees Celsius) and date for the Western Australian region for the period 1982-01-01 to 2014-12-31.
sst_WA
sst_WA
A data frame with 12053 rows and 2 variables:
date, as POSIXct
SST, in degrees Celsius
...
https://www.ncdc.noaa.gov/oisst