Marine Heatwaves Definition¶
This document describes the marineHeatWaves module for python which implements the Marine Heatwave (MHW) definition of Hobday et al. (manuscript in preparation). This software 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. This was developed as an interactive ipython notebook which the user can run themselves provided they have the required python modules (numpy, scipy, datetime, and matplotlib).
In [1]:
# Load required modules
import numpy as np
from datetime import date
from matplotlib import pyplot as plt
%pylab inline
In [2]:
# Load marineHeatWaves definition module
import marineHeatWaves as mhw
Marine Heatwave Detection¶
The marineHeatWaves (mhw) module consists of a number of functions for the detection and characterization of MHWs. The main function is the detection function (detect) which takes as input a time series of temperature (and a corresponding time vector) and outputs a set of detected MHWs.
Case study #1: 2011 Western Australia event¶
As an example, let's load a daily time series of SST off Western Australia (WA; 112.5$^\circ$E, 29.5$^\circ$S) over the 1982 to 2014 period, remotely-sensed from the AVHRR satellite platform:
In [3]:
# Load WA SST time series
sst = np.loadtxt('data/sst_WA.csv', delimiter=',')
# Generate time vector using datetime format (January 1 of year 1 is day 1)
t = np.arange(date(1982,1,1).toordinal(),date(2014,12,31).toordinal()+1)
dates = [date.fromordinal(tt.astype(int)) for tt in t]
Next we run the MHW detection algorithm which returns the variable mhws, consisting of the detected MHWs, and clim, consisting of the climatological (varying by day-of-year) seasonal cycle and extremes threshold:
In [4]:
mhws, clim = mhw.detect(t, sst)
This algorithm has detected the following number of MHW events:
In [5]:
mhws['n_events']
Out[5]:
The first ten events, for example, have the following maximum intensities (in $^\circ$C):
In [6]:
mhws['intensity_max'][0:10]
Out[6]:
Let's have a look at some properties of the event with the largest maximum intensity
In [7]:
ev = np.argmax(mhws['intensity_max']) # Find largest event
print 'Maximum intensity:', mhws['intensity_max'][ev], 'deg. C'
print 'Average intensity:', mhws['intensity_mean'][ev], 'deg. C'
print 'Cumulative intensity:', mhws['intensity_cumulative'][ev], 'deg. C-days'
print 'Duration:', mhws['duration'][ev], 'days'
print 'Start date:', mhws['date_start'][ev].strftime("%d %B %Y")
print 'End date:', mhws['date_end'][ev].strftime("%d %B %Y")
This turns out to be the infamous 2011 MHW off WA. Let's plot the SST time series over the full record and also have a closer look at the identified MHW event:
In [8]:
plt.figure(figsize=(14,10))
plt.subplot(2,1,1)
# Plot SST, seasonal cycle, and threshold
plt.plot(dates, sst, 'k-')
plt.plot(dates, clim['thresh'], 'g-')
plt.plot(dates, clim['seas'], 'b-')
plt.title('SST (black), seasonal climatology (blue), \
threshold (green), detected MHW events (shading)')
plt.xlim(t[0], t[-1])
plt.ylim(sst.min()-0.5, sst.max()+0.5)
plt.ylabel(r'SST [$^\circ$C]')
plt.subplot(2,1,2)
# Find indices for all ten MHWs before and after event of interest and shade accordingly
for ev0 in np.arange(ev-10, ev+11, 1):
t1 = np.where(t==mhws['time_start'][ev0])[0][0]
t2 = np.where(t==mhws['time_end'][ev0])[0][0]
plt.fill_between(dates[t1:t2+1], sst[t1:t2+1], clim['thresh'][t1:t2+1], \
color=(1,0.6,0.5))
# Find indices for MHW of interest (2011 WA event) and shade accordingly
t1 = np.where(t==mhws['time_start'][ev])[0][0]
t2 = np.where(t==mhws['time_end'][ev])[0][0]
plt.fill_between(dates[t1:t2+1], sst[t1:t2+1], clim['thresh'][t1:t2+1], \
color='r')
# Plot SST, seasonal cycle, threshold, shade MHWs with main event in red
plt.plot(dates, sst, 'k-', linewidth=2)
plt.plot(dates, clim['thresh'], 'g-', linewidth=2)
plt.plot(dates, clim['seas'], 'b-', linewidth=2)
plt.title('SST (black), seasonal climatology (blue), \
threshold (green), detected MHW events (shading)')
plt.xlim(mhws['time_start'][ev]-150, mhws['time_end'][ev]+150)
plt.ylim(clim['seas'].min() - 1, clim['seas'].max() + mhws['intensity_max'][ev] + 0.5)
plt.ylabel(r'SST [$^\circ$C]')
Out[8]:
Yep, It's certainly picked out the largest event in the series (dark red shading). This event also seems to have been preceded and succeeded by a number of shorter, weaker events (light red shading). Let's have a look at how the MHW statistics are distributed across all the detected events:
In [9]:
plt.figure(figsize=(15,7))
# Duration
plt.subplot(2,2,1)
evMax = np.argmax(mhws['duration'])
plt.bar(range(mhws['n_events']), mhws['duration'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['duration'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['duration'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.ylabel('[days]')
plt.title('Duration')
# Maximum intensity
plt.subplot(2,2,2)
evMax = np.argmax(mhws['intensity_max'])
plt.bar(range(mhws['n_events']), mhws['intensity_max'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_max'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_max'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.ylabel(r'[$^\circ$C]')
plt.title('Maximum Intensity')
# Mean intensity
plt.subplot(2,2,4)
evMax = np.argmax(mhws['intensity_mean'])
plt.bar(range(mhws['n_events']), mhws['intensity_mean'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_mean'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_mean'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.title('Mean Intensity')
plt.ylabel(r'[$^\circ$C]')
plt.xlabel('MHW event number')
# Cumulative intensity
plt.subplot(2,2,3)
evMax = np.argmax(mhws['intensity_cumulative'])
plt.bar(range(mhws['n_events']), mhws['intensity_cumulative'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_cumulative'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_cumulative'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.title(r'Cumulative Intensity')
plt.ylabel(r'[$^\circ$C$\times$days]')
plt.xlabel('MHW event number')
Out[9]:
The largest event on record by duration, maximum intensity, mean intensity, and cumulative intensity (shown by red shaded bars) happens to be the 2011 event (red outlined bars).
Case study #2: 2012 Northwest Atlantic event¶
We can also have a look at the 2012 Northwest Atlantic (NWA) event. Let's load a daily time series of SST from the NWA (67$^\circ$W, 43$^\circ$N) from the same data set as above and repeat the analysis:
In [10]:
# Load Northwest Atlantic SST time series
sst = np.loadtxt('data/sst_NW_Atl.csv', delimiter=',')
# Detect MHW events
mhws, clim = mhw.detect(t, sst)
# Find largest event
ev = np.argmax(mhws['intensity_max'])
print 'Maximum intensity:', mhws['intensity_max'][ev], 'deg. C'
print 'Average intensity:', mhws['intensity_mean'][ev], 'deg. C'
print 'Cumulative intensity:', mhws['intensity_cumulative'][ev], 'deg. C-days'
print 'Duration:', mhws['duration'][ev], 'days'
print 'Start date:', mhws['date_start'][ev].strftime("%d %B %Y")
print 'End date:', mhws['date_end'][ev].strftime("%d %B %Y")
Again, the event with largest maximum intensity turns out to be an event known in the literature: the 2012 NWA event. Let's plot it as above:
In [11]:
plt.figure(figsize=(14,10))
plt.subplot(2,1,1)
# Plot SST, seasonal cycle, and threshold
plt.plot(dates, sst, 'k-')
plt.plot(dates, clim['thresh'], 'g-')
plt.plot(dates, clim['seas'], 'b-')
plt.title('SST (black), seasonal climatology (blue), threshold (green), \
detected MHW events (shading)')
plt.xlim(t[0], t[-1])
plt.ylim(sst.min()-0.5, sst.max()+0.5)
plt.ylabel(r'SST [$^\circ$C]')
plt.subplot(2,1,2)
# Find indices for all ten MHWs before and after event of interest and shade accordingly
for ev0 in np.arange(ev-10, ev+11, 1):
t1 = np.where(t==mhws['time_start'][ev0])[0][0]
t2 = np.where(t==mhws['time_end'][ev0])[0][0]
plt.fill_between(dates[t1:t2+1], sst[t1:t2+1], clim['thresh'][t1:t2+1], \
color=(1,0.6,0.5))
# Find indices for MHW of interest (2012 NWA event) and shade accordingly
t1 = np.where(t==mhws['time_start'][ev])[0][0]
t2 = np.where(t==mhws['time_end'][ev])[0][0]
plt.fill_between(dates[t1:t2+1], sst[t1:t2+1], clim['thresh'][t1:t2+1], color='r')
# Plot SST, seasonal cycle, threshold, shade MHWs with main event in red
plt.plot(dates, sst, 'k-', linewidth=2)
plt.plot(dates, clim['thresh'], 'g-', linewidth=2)
plt.plot(dates, clim['seas'], 'b-', linewidth=2)
plt.title('SST (black), seasonal climatology (blue), threshold (green), \
detected MHW events (shading)')
plt.xlim(mhws['time_start'][ev]-150, mhws['time_end'][ev]+150)
plt.ylim(clim['seas'].min() - 1, clim['seas'].max() + \
mhws['intensity_max'][ev] + 0.5)
plt.ylabel(r'SST [$^\circ$C]')
Out[11]:
Interestingly, the 2012 NWA event doesn't appear to be as intense as the 2011 WA event. That is partly because it occured outside the summer months. It is nonetheless the largest event on record according to maximum intensity:
In [12]:
plt.figure(figsize=(15,7))
# Duration
plt.subplot(2,2,1)
evMax = np.argmax(mhws['duration'])
plt.bar(range(mhws['n_events']), mhws['duration'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['duration'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['duration'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.ylabel('[days]')
plt.title('Duration')
# Maximum intensity
plt.subplot(2,2,2)
evMax = np.argmax(mhws['intensity_max'])
plt.bar(range(mhws['n_events']), mhws['intensity_max'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_max'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_max'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.ylabel(r'[$^\circ$C]')
plt.title('Maximum Intensity')
# Mean intensity
plt.subplot(2,2,4)
evMax = np.argmax(mhws['intensity_mean'])
plt.bar(range(mhws['n_events']), mhws['intensity_mean'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_mean'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_mean'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.title('Mean Intensity')
plt.ylabel(r'[$^\circ$C]')
plt.xlabel('MHW event number')
# Cumulative intensity
plt.subplot(2,2,3)
evMax = np.argmax(mhws['intensity_cumulative'])
plt.bar(range(mhws['n_events']), mhws['intensity_cumulative'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_cumulative'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_cumulative'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.title(r'Cumulative Intensity')
plt.ylabel(r'[$^\circ$C$\times$days]')
plt.xlabel('MHW event number')
Out[12]:
However, the 2012 NWA event (red outlined bars) is not the highest ranked event (red shaded bars) when sorted according to mean intensity, cumulative intenstiy and duration:
In [13]:
ranks = mhws['n_events']- np.array(mhws['duration']).argsort().argsort()
print "The 2012 NWA event is ranked number " \
+ str(ranks[ev]) + " by duration"
ranks = mhws['n_events']- np.array(mhws['intensity_max']).argsort().argsort()
print "The 2012 NWA event is ranked number " \
+ str(ranks[ev]) + " by maximum intensity"
ranks = mhws['n_events']- np.array(mhws['intensity_mean']).argsort().argsort()
print "The 2012 NWA event is ranked number " \
+ str(ranks[ev]) + " by mean intensity"
ranks = mhws['n_events']- np.array(mhws['intensity_cumulative']).argsort().argsort()
print "The 2012 NWA event is ranked number " \
+ str(ranks[ev]) + " by cumulative intensity"
Case study #3: 2003 Mediterranean event¶
We can also have a look at the 2003 Mediterranean (Med) event. Let's load a daily time series of SST from the Med just south of France (9$^\circ$E, 43.5$^\circ$N) from the same data set as above and repeat the analysis:
In [14]:
# Load Mediterranean SST time series
sst = np.loadtxt('data/sst_Med.csv', delimiter=',')
# Detect MHW events
mhws, clim = mhw.detect(t, sst)
# Find largest event
ev = np.argmax(mhws['intensity_max'])
print 'Maximum intensity:', mhws['intensity_max'][ev], 'deg. C'
print 'Average intensity:', mhws['intensity_mean'][ev], 'deg. C'
print 'Cumulative intensity:', mhws['intensity_cumulative'][ev], 'deg. C-days'
print 'Duration:', mhws['duration'][ev], 'days'
print 'Start date:', mhws['date_start'][ev].strftime("%d %B %Y")
print 'End date:', mhws['date_end'][ev].strftime("%d %B %Y")
We are looking for the 2003 Med event but the largest event on record, according to maximum intensity, occurred in 2008. However, the 2003 Med event is the largest according to mean intensity:
In [15]:
# Find largest event (by mean intensity)
ev = np.argmax(mhws['intensity_mean'])
print 'Maximum intensity:', mhws['intensity_max'][ev], 'deg. C'
print 'Average intensity:', mhws['intensity_mean'][ev], 'deg. C'
print 'Cumulative intensity:', mhws['intensity_cumulative'][ev], 'deg. C-days'
print 'Duration:', mhws['duration'][ev], 'days'
print 'Start date:', mhws['date_start'][ev].strftime("%d %B %Y")
print 'End date:', mhws['date_end'][ev].strftime("%d %B %Y")
And a time series plot shows it clearly to be a large event:
In [16]:
plt.figure(figsize=(14,10))
plt.subplot(2,1,1)
# Plot SST, seasonal cycle, and threshold
plt.plot(dates, sst, 'k-')
plt.plot(dates, clim['thresh'], 'g-')
plt.plot(dates, clim['seas'], 'b-')
plt.title('SST (black), seasonal climatology (blue), threshold (green), \
detected MHW events (shading)')
plt.xlim(t[0], t[-1])
plt.ylim(sst.min()-0.5, sst.max()+0.5)
plt.ylabel(r'SST [$^\circ$C]')
plt.subplot(2,1,2)
# Find indices for all ten MHWs before and after event of interest and shade accordingly
for ev0 in np.arange(ev-10, ev+11, 1):
t1 = np.where(t==mhws['time_start'][ev0])[0][0]
t2 = np.where(t==mhws['time_end'][ev0])[0][0]
plt.fill_between(dates[t1:t2+1], sst[t1:t2+1], clim['thresh'][t1:t2+1], \
color=(1,0.6,0.5))
# Find indices for MHW of interest (2003 Med event) and shade accordingly
t1 = np.where(t==mhws['time_start'][ev])[0][0]
t2 = np.where(t==mhws['time_end'][ev])[0][0]
plt.fill_between(dates[t1:t2+1], sst[t1:t2+1], clim['thresh'][t1:t2+1], color='r')
# Plot SST, seasonal cycle, threshold, shade MHWs with main event in red
plt.plot(dates, sst, 'k-', linewidth=2)
plt.plot(dates, clim['thresh'], 'g-', linewidth=2)
plt.plot(dates, clim['seas'], 'b-', linewidth=2)
plt.title('SST (black), seasonal climatology (blue), threshold (green), \
detected MHW events (shading)')
plt.xlim(mhws['time_start'][ev]-150, mhws['time_end'][ev]+150)
plt.ylim(clim['seas'].min() - 1, clim['seas'].max()\
+ mhws['intensity_max'][ev] + 0.5)
plt.ylabel(r'SST [$^\circ$C]')
Out[16]:
And we can see where this event (red outlined bars) fits amongst all other detected events:
In [17]:
plt.figure(figsize=(15,7))
# Duration
plt.subplot(2,2,1)
evMax = np.argmax(mhws['duration'])
plt.bar(range(mhws['n_events']), mhws['duration'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['duration'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['duration'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.ylabel('[days]')
plt.title('Duration')
# Maximum intensity
plt.subplot(2,2,2)
evMax = np.argmax(mhws['intensity_max'])
plt.bar(range(mhws['n_events']), mhws['intensity_max'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_max'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_max'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.ylabel(r'[$^\circ$C]')
plt.title('Maximum Intensity')
# Mean intensity
plt.subplot(2,2,4)
evMax = np.argmax(mhws['intensity_mean'])
plt.bar(range(mhws['n_events']), mhws['intensity_mean'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_mean'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_mean'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.title('Mean Intensity')
plt.ylabel(r'[$^\circ$C]')
plt.xlabel('MHW event number')
# Cumulative intensity
plt.subplot(2,2,3)
evMax = np.argmax(mhws['intensity_cumulative'])
plt.bar(range(mhws['n_events']), mhws['intensity_cumulative'], width=0.6, \
color=(0.7,0.7,0.7))
plt.bar(evMax, mhws['intensity_cumulative'][evMax], width=0.6, color=(1,0.5,0.5))
plt.bar(ev, mhws['intensity_cumulative'][ev], width=0.6, edgecolor=(1,0.,0.), \
color='none')
plt.xlim(0, mhws['n_events'])
plt.title(r'Cumulative Intensity')
plt.ylabel(r'[$^\circ$C$\times$days]')
plt.xlabel('MHW event number')
Out[17]:
In [18]:
ranks = mhws['n_events']- np.array(mhws['duration']).argsort().argsort()
print "The 2003 Med event is ranked number " \
+ str(ranks[ev]) + " by duration"
ranks = mhws['n_events']- np.array(mhws['intensity_max']).argsort().argsort()
print "The 2003 Med event is ranked number " \
+ str(ranks[ev]) + " by maximum intensity"
ranks = mhws['n_events']- np.array(mhws['intensity_mean']).argsort().argsort()
print "The 2003 Med event is ranked number " \
+ str(ranks[ev]) + " by mean intensity"
ranks = mhws['n_events']- np.array(mhws['intensity_cumulative']).argsort().argsort()
print "The 2003 Med event is ranked number " \
+ str(ranks[ev]) + " by cumulative intensity"
Block-averaged Marine Heatwave properties¶
The marineHeatWaves (mhw) module also consists of functions to calculate the average of MHW properties over blocks in time (e.g., annually, decadally). The block-averaging function (blockAverage) takes as input a set of detected MHWs (i.e., the output from detect, the detection function described above) and outputs the MHW properties averaged over the specified block-length. This output can then be passed through the meanTrend function in order to calculate the time-mean and linear trend of the MHW properties over the measurement period.
Let's start by applying the block-averaging function to the Mediterranean MHWs which are stored in the variable mhws, using the default block length of 1 year (i.e., annual averages):
In [19]:
mhwBlock = mhw.blockAverage(t, mhws)
The variable mhwBlock has a set of keys which are time series of the MHW properties over the blocks. The central year of the blocks are stored in the key years_centre so we can look at, as an example, time series of MHW counts in each year and the average maximum intensity in each year:
In [20]:
plt.figure(figsize=(14,4))
plt.subplot(1,2,1)
plt.plot(mhwBlock['years_centre'], mhwBlock['count'], 'k-o')
plt.ylim(0,9)
plt.ylabel('[count]')
plt.title('Number of MHWs by year')
plt.subplot(1,2,2)
plt.plot(mhwBlock['years_centre'], mhwBlock['intensity_max'], 'k-o')
plt.ylabel(r'[$^\circ$C]')
plt.title('Average MHW maximum intensity by year')
Out[20]:
There certainly looks like a positive trend in the annual number of MHWs and possibly in the maximum intensity as well. We can calculate the mean and trend of the MHW properties using the meanTrend function:
In [22]:
mean, trend = mhw.meanTrend(mhwBlock)
print "There are on average " + str(mean['count']) + " MHWs in each year, \n \
with a linear trend of " + str(10*trend['count']) + " MHW events per decade\n"
print "The average maximum intensity is " + str(mean['intensity_max']) + " deg. C, \n \
with a linear trend of " + str(10*trend['intensity_max']) + " deg. C per decade"
Appendix: Function documentation¶
detect¶
def detect(t, sst, climatologyPeriod=[1983,2012], pctile=90, windowHalfWidth=5, smoothPercentile=True, smoothPercentileWidth=31, minDuration=5, joinAcrossGaps=True, maxGap=2):
Applies the Hobday et al. (2015) marine heat wave definition to an input time
series of sst ('sst') along with a time vector ('t'). Outputs properties of
all detected marine heat waves.
Inputs:
t Time vector, in datetime format (e.g., date(1982,1,1).toordinal())
[1D numpy array of length T]
sst Temperature vector [1D numpy array of length T]
Outputs:
mhw Detected marine heat waves (MHWs). Each key (following list) is a
list of length N where N is the number of detected MHWs:
'time_start' Start time of MHW [datetime format]
'time_end' End time of MHW [datetime format]
'time_peak' Time of MHW peak [datetime format]
'date_start' Start date of MHW [datetime format]
'date_end' End date of MHW [datetime format]
'date_peak' Date of MHW peak [datetime format]
'duration' Duration of MHW [days]
'intensity_max' Maximum (peak) intensity [deg. C]
'intensity_mean' Mean intensity [deg. C]
'intensity_var' Intensity variability [deg. C]
'intensity_cumulative' Cumulative intensity [deg. C x days]
'rate_onset' Onset rate of MHW [deg. C / days]
'rate_decline' Decline rate of MHW [deg. C / days]
'intensity_max_relThresh', 'intensity_mean_relThresh', 'intensity_var_relThresh',
and 'intensity_cumulative_relThresh' are as above except relative to the
threshold (e.g., 90th percentile) rather than the seasonal climatology
'n_events' A scalar integer (not a list) indicating the total
number of detected MHW events
clim Climatology of SST. Each key (following list) is a seasonally-varying
time series [1D numpy array of length T] of a particular measure:
'thresh' Seasonally varying threshold (e.g., 90th percentile)
'seas' Climatological seasonal cycle
Options:
climatologyPeriod Period over which climatology is calculated, specified
as list of start and end years (DEFAULT = [1983,2012])
pctile Threshold percentile (%) for detection of extreme values
(DEFAULT = 90)
windowHalfWidth Width of window (one sided) about day-of-year used for
the pooling of values and calculation of threshold percentile
(DEFAULT = 5 [days])
smoothPercentile Boolean switch indicating whether to smooth the threshold
percentile timeseries with a moving average (DEFAULT = True)
smoothPercentileWidth Width of moving average window for smoothing threshold
(DEFAULT = 31 [days])
minDuration Minimum duration for acceptance detected MHWs
(DEFAULT = 5 [days])
joinAcrossGaps Boolean switch indicating whether to join MHWs
which occur before/after a short gap (DEFAULT = True)
maxGap Maximum length of gap allowed for the joining of MHWs
(DEFAULT = 2 [days])
Notes:
1. This function assumes that the input time series consist of continuous daily values
with no missing values.
2. This function supports leap years. This is done by ignoring Feb 29s for the initial
calculation of the climatology and threshold. The value of these for Feb 29 is then
linearly interpolated from the values for Feb 28 and Mar 1.
3. The calculation of onset and decline rates assumes that the heat wave 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.)
Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Feb 2015blockAverage¶
def blockAverage(t, mhw, blockLength=1):
Calculate statistics of marine heatwave (MHW) properties averaged over blocks of
a specified length of time. Takes as input a collection of detected MHWs
(using the marineHeatWaves.detect function) and a time vector for the source
SST series.
Inputs:
t Time vector, in datetime format (e.g., date(1982,1,1).toordinal())
mhw Marine heat waves (MHWs) detected using marineHeatWaves.detect
Outputs:
mhwBlock Time series of block-averaged MHW properties. Each key (following list)
is a list of length N where N is the number of blocks:
'years_start' Start year blocks (inclusive)
'years_end' End year of blocks (inclusive)
'years_centre' Decimal year at centre of blocks
'count' Total MHW count in each block
'duration' Average MHW duration in each block [days]
'intensity_max' Average MHW "maximum (peak) intensity" in each block [deg. C]
'intensity_mean' Average MHW "mean intensity" in each block [deg. C]
'intensity_var' Average MHW "intensity variability" in each block [deg. C]
'intensity_cumulative' Average MHW "cumulative intensity" in each block [deg. C x days]
'rate_onset' Average MHW onset rate in each block [deg. C / days]
'rate_decline' Average MHW decline rate in each block [deg. C / days]
'intensity_max_relThresh', 'intensity_mean_relThresh', 'intensity_var_relThresh',
and 'intensity_cumulative_relThresh' are as above except relative to the
threshold (e.g., 90th percentile) rather than the seasonal climatology
Options:
blockLength Size of block (in years) over which to calculate the
averaged MHW properties. Must be an integer greater than
or equal to 1 (DEFAULT = 1 [year])
Notes:
This function assumes that the input time vector consists of continuous daily values.
Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Feb-Mar 2015meanTrend¶
def meanTrend(mhwBlock):
Calculates the mean and trend of marine heatwave (MHW) properties. Takes as input a
collection of block-averaged MHW properties (using the marineHeatWaves.blockAverage
function). Handles missing values (which should be specified by NaNs).
Inputs:
mhwBlock Time series of block-averaged MHW statistics calculated using the
marineHeatWaves.blockAverage function
Outputs:
mean Mean of all MHW properties over all block-averaged values
trend Linear trend of all MHW properties over all block-averaged values
Both mean and trend have the following keys, the units the trend
are the units of the property of interest per year:
'duration' Duration of MHW [days]
'intensity_max' Maximum (peak) intensity [deg. C]
'intensity_mean' Mean intensity [deg. C]
'intensity_var' Intensity variability [deg. C]
'intensity_cumulative' Cumulative intensity [deg. C x days]
'rate_onset' Onset rate of MHW [deg. C / days]
'rate_decline' Decline rate of MHW [deg. C / days]
'intensity_max_relThresh', 'intensity_mean_relThresh', 'intensity_var_relThresh',
and 'intensity_cumulative_relThresh' are as above except relative to the
threshold (e.g., 90th percentile) rather than the seasonal climatology
Notes:
This calculation performs a multiple linear regression of the form
y ~ beta * X + eps
where y is the MHW property of interest and X is a matrix of predictors. The first
column of X is all ones to estimate the mean, the second column is the time vector
which is taken as mhwBlock['years_centre'] and offset to be equal to zero at its
mid-point.
Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Feb-Mar 2015In [ ]: