''' A set of functions which implement the Marine Heat Wave (MHW) definition of Hobday et al. (2016) ''' import numpy as np import scipy as sp from scipy import linalg from scipy import stats import scipy.ndimage as ndimage from datetime import date def detect(t, temp, climatologyPeriod=[None,None], pctile=90, windowHalfWidth=5, smoothPercentile=True, smoothPercentileWidth=31, minDuration=5, joinAcrossGaps=True, maxGap=2, maxPadLength=False, coldSpells=False, alternateClimatology=False, Ly=False): ''' Applies the Hobday et al. (2016) marine heat wave definition to an input time series of temp ('temp') 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] temp 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] 'index_start' Start index of MHW 'index_end' End index of MHW 'index_peak' Index of MHW peak '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 'intensity_max_abs', 'intensity_mean_abs', 'intensity_var_abs', and 'intensity_cumulative_abs' are as above except as absolute magnitudes rather than relative to the seasonal climatology or threshold 'category' is an integer category system (1, 2, 3, 4) based on the maximum intensity in multiples of threshold exceedances, i.e., a value of 1 indicates the MHW intensity (relative to the climatology) was >=1 times the value of the threshold (but less than 2 times; relative to climatology, i.e., threshold - climatology). Category types are defined as 1=moderate, 2=strong, 3=severe, 4=extreme. More details in Hobday et al. (in prep., Oceanography). Also supplied are the duration of each of these categories for each event. '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 'missing' A vector of TRUE/FALSE indicating which elements in temp were missing values for the MHWs detection Options: climatologyPeriod Period over which climatology is calculated, specified as list of start and end years. Default is to calculate over the full range of years in the supplied time series. Alternate periods suppled as a list e.g. [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]) maxPadLength Specifies the maximum length [days] over which to interpolate (pad) missing data (specified as nans) in input temp time series. i.e., any consecutive blocks of NaNs with length greater than maxPadLength will be left as NaN. Set as an integer. (DEFAULT = False, interpolates over all missing values). coldSpells Specifies if the code should detect cold events instead of heat events. (DEFAULT = False) alternateClimatology Specifies an alternate temperature time series to use for the calculation of the climatology. Format is as a list of numpy arrays: (1) the first element of the list is a time vector, in datetime format (e.g., date(1982,1,1).toordinal()) [1D numpy array of length TClim] and (2) the second element of the list is a temperature vector [1D numpy array of length TClim]. (DEFAULT = False) Ly Specifies if the length of the year is < 365/366 days (e.g. a 360 day year from a climate model). This affects the calculation of the climatology. (DEFAULT = False) Notes: 1. This function assumes that the input time series consist of continuous daily values with few missing values. Time ranges which start and end part-way through the calendar year are supported. 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.) 4. For the purposes of MHW detection, any missing temp values not interpolated over (through optional maxPadLLength) will be set equal to the seasonal climatology. This means they will trigger the end/start of any adjacent temp values which satisfy the MHW criteria. 5. 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. Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Feb 2015 ''' # # Initialize MHW output variable # mhw = {} mhw['time_start'] = [] # datetime format mhw['time_end'] = [] # datetime format mhw['time_peak'] = [] # datetime format mhw['date_start'] = [] # datetime format mhw['date_end'] = [] # datetime format mhw['date_peak'] = [] # datetime format mhw['index_start'] = [] mhw['index_end'] = [] mhw['index_peak'] = [] mhw['duration'] = [] # [days] mhw['duration_moderate'] = [] # [days] mhw['duration_strong'] = [] # [days] mhw['duration_severe'] = [] # [days] mhw['duration_extreme'] = [] # [days] mhw['intensity_max'] = [] # [deg C] mhw['intensity_mean'] = [] # [deg C] mhw['intensity_var'] = [] # [deg C] mhw['intensity_cumulative'] = [] # [deg C] mhw['intensity_max_relThresh'] = [] # [deg C] mhw['intensity_mean_relThresh'] = [] # [deg C] mhw['intensity_var_relThresh'] = [] # [deg C] mhw['intensity_cumulative_relThresh'] = [] # [deg C] mhw['intensity_max_abs'] = [] # [deg C] mhw['intensity_mean_abs'] = [] # [deg C] mhw['intensity_var_abs'] = [] # [deg C] mhw['intensity_cumulative_abs'] = [] # [deg C] mhw['category'] = [] mhw['rate_onset'] = [] # [deg C / day] mhw['rate_decline'] = [] # [deg C / day] # # Time and dates vectors # # Generate vectors for year, month, day-of-month, and day-of-year T = len(t) year = np.zeros((T)) month = np.zeros((T)) day = np.zeros((T)) doy = np.zeros((T)) for i in range(T): year[i] = date.fromordinal(t[i]).year month[i] = date.fromordinal(t[i]).month day[i] = date.fromordinal(t[i]).day # Leap-year baseline for defining day-of-year values year_leapYear = 2012 # This year was a leap-year and therefore doy in range of 1 to 366 t_leapYear = np.arange(date(year_leapYear, 1, 1).toordinal(),date(year_leapYear, 12, 31).toordinal()+1) dates_leapYear = [date.fromordinal(tt.astype(int)) for tt in t_leapYear] month_leapYear = np.zeros((len(t_leapYear))) day_leapYear = np.zeros((len(t_leapYear))) doy_leapYear = np.zeros((len(t_leapYear))) for tt in range(len(t_leapYear)): month_leapYear[tt] = date.fromordinal(t_leapYear[tt]).month day_leapYear[tt] = date.fromordinal(t_leapYear[tt]).day doy_leapYear[tt] = t_leapYear[tt] - date(date.fromordinal(t_leapYear[tt]).year,1,1).toordinal() + 1 # Calculate day-of-year values for tt in range(T): doy[tt] = doy_leapYear[(month_leapYear == month[tt]) * (day_leapYear == day[tt])] # Constants (doy values for Feb-28 and Feb-29) for handling leap-years feb28 = 59 feb29 = 60 # Set climatology period, if unset use full range of available data if (climatologyPeriod[0] is None) or (climatologyPeriod[1] is None): climatologyPeriod[0] = year[0] climatologyPeriod[1] = year[-1] # # Calculate threshold and seasonal climatology (varying with day-of-year) # # if alternate temperature time series is supplied for the calculation of the climatology if alternateClimatology: tClim = alternateClimatology[0] tempClim = alternateClimatology[1] TClim = len(tClim) yearClim = np.zeros((TClim)) monthClim = np.zeros((TClim)) dayClim = np.zeros((TClim)) doyClim = np.zeros((TClim)) for i in range(TClim): yearClim[i] = date.fromordinal(tClim[i]).year monthClim[i] = date.fromordinal(tClim[i]).month dayClim[i] = date.fromordinal(tClim[i]).day doyClim[i] = doy_leapYear[(month_leapYear == monthClim[i]) * (day_leapYear == dayClim[i])] else: tempClim = temp.copy() TClim = np.array([T]).copy()[0] yearClim = year.copy() monthClim = month.copy() dayClim = day.copy() doyClim = doy.copy() # Flip temp time series if detecting cold spells if coldSpells: temp = -1.*temp tempClim = -1.*tempClim # Pad missing values for all consecutive missing blocks of length <= maxPadLength if maxPadLength: temp = pad(temp, maxPadLength=maxPadLength) tempClim = pad(tempClim, maxPadLength=maxPadLength) # Length of climatological year lenClimYear = 366 # Start and end indices clim_start = np.where(yearClim == climatologyPeriod[0])[0][0] clim_end = np.where(yearClim == climatologyPeriod[1])[0][-1] # Inialize arrays thresh_climYear = np.NaN*np.zeros(lenClimYear) seas_climYear = np.NaN*np.zeros(lenClimYear) clim = {} clim['thresh'] = np.NaN*np.zeros(TClim) clim['seas'] = np.NaN*np.zeros(TClim) # Loop over all day-of-year values, and calculate threshold and seasonal climatology across years for d in range(1,lenClimYear+1): # Special case for Feb 29 if d == feb29: continue # find all indices for each day of the year +/- windowHalfWidth and from them calculate the threshold tt0 = np.where(doyClim[clim_start:clim_end+1] == d)[0] # If this doy value does not exist (i.e. in 360-day calendars) then skip it if len(tt0) == 0: continue tt = np.array([]) for w in range(-windowHalfWidth, windowHalfWidth+1): tt = np.append(tt, clim_start+tt0 + w) tt = tt[tt>=0] # Reject indices "before" the first element tt = tt[tt= 365-day year else: thresh_climYear = runavg(thresh_climYear, smoothPercentileWidth) seas_climYear = runavg(seas_climYear, smoothPercentileWidth) # Generate threshold for full time series clim['thresh'] = thresh_climYear[doy.astype(int)-1] clim['seas'] = seas_climYear[doy.astype(int)-1] # Save vector indicating which points in temp are missing values clim['missing'] = np.isnan(temp) # Set all remaining missing temp values equal to the climatology temp[np.isnan(temp)] = clim['seas'][np.isnan(temp)] # # Find MHWs as exceedances above the threshold # # Time series of "True" when threshold is exceeded, "False" otherwise exceed_bool = temp - clim['thresh'] exceed_bool[exceed_bool<=0] = False exceed_bool[exceed_bool>0] = True # Fix issue where missing temp vaues (nan) are counted as True exceed_bool[np.isnan(exceed_bool)] = False # Find contiguous regions of exceed_bool = True events, n_events = ndimage.label(exceed_bool) # Find all MHW events of duration >= minDuration for ev in range(1,n_events+1): event_duration = (events == ev).sum() if event_duration < minDuration: continue mhw['time_start'].append(t[np.where(events == ev)[0][0]]) mhw['time_end'].append(t[np.where(events == ev)[0][-1]]) # Link heat waves that occur before and after a short gap (gap must be no longer than maxGap) if joinAcrossGaps: # Calculate gap length for each consecutive pair of events gaps = np.array(mhw['time_start'][1:]) - np.array(mhw['time_end'][0:-1]) - 1 if len(gaps) > 0: while gaps.min() <= maxGap: # Find first short gap ev = np.where(gaps <= maxGap)[0][0] # Extend first MHW to encompass second MHW (including gap) mhw['time_end'][ev] = mhw['time_end'][ev+1] # Remove second event from record del mhw['time_start'][ev+1] del mhw['time_end'][ev+1] # Calculate gap length for each consecutive pair of events gaps = np.array(mhw['time_start'][1:]) - np.array(mhw['time_end'][0:-1]) - 1 if len(gaps) == 0: break # Calculate marine heat wave properties mhw['n_events'] = len(mhw['time_start']) categories = np.array(['Moderate', 'Strong', 'Severe', 'Extreme']) for ev in range(mhw['n_events']): mhw['date_start'].append(date.fromordinal(mhw['time_start'][ev])) mhw['date_end'].append(date.fromordinal(mhw['time_end'][ev])) # Get SST series during MHW event, relative to both threshold and to seasonal climatology tt_start = np.where(t==mhw['time_start'][ev])[0][0] tt_end = np.where(t==mhw['time_end'][ev])[0][0] mhw['index_start'].append(tt_start) mhw['index_end'].append(tt_end) temp_mhw = temp[tt_start:tt_end+1] thresh_mhw = clim['thresh'][tt_start:tt_end+1] seas_mhw = clim['seas'][tt_start:tt_end+1] mhw_relSeas = temp_mhw - seas_mhw mhw_relThresh = temp_mhw - thresh_mhw mhw_relThreshNorm = (temp_mhw - thresh_mhw) / (thresh_mhw - seas_mhw) mhw_abs = temp_mhw # Find peak tt_peak = np.argmax(mhw_relSeas) mhw['time_peak'].append(mhw['time_start'][ev] + tt_peak) mhw['date_peak'].append(date.fromordinal(mhw['time_start'][ev] + tt_peak)) mhw['index_peak'].append(tt_start + tt_peak) # MHW Duration mhw['duration'].append(len(mhw_relSeas)) # MHW Intensity metrics mhw['intensity_max'].append(mhw_relSeas[tt_peak]) mhw['intensity_mean'].append(mhw_relSeas.mean()) mhw['intensity_var'].append(np.sqrt(mhw_relSeas.var())) mhw['intensity_cumulative'].append(mhw_relSeas.sum()) mhw['intensity_max_relThresh'].append(mhw_relThresh[tt_peak]) mhw['intensity_mean_relThresh'].append(mhw_relThresh.mean()) mhw['intensity_var_relThresh'].append(np.sqrt(mhw_relThresh.var())) mhw['intensity_cumulative_relThresh'].append(mhw_relThresh.sum()) mhw['intensity_max_abs'].append(mhw_abs[tt_peak]) mhw['intensity_mean_abs'].append(mhw_abs.mean()) mhw['intensity_var_abs'].append(np.sqrt(mhw_abs.var())) mhw['intensity_cumulative_abs'].append(mhw_abs.sum()) # Fix categories tt_peakCat = np.argmax(mhw_relThreshNorm) cats = np.floor(1. + mhw_relThreshNorm) mhw['category'].append(categories[np.min([cats[tt_peakCat], 4]).astype(int) - 1]) mhw['duration_moderate'].append(np.sum(cats == 1.)) mhw['duration_strong'].append(np.sum(cats == 2.)) mhw['duration_severe'].append(np.sum(cats == 3.)) mhw['duration_extreme'].append(np.sum(cats >= 4.)) # Rates of onset and decline # Requires getting MHW strength at "start" and "end" of event (continuous: assume start/end half-day before/after first/last point) if tt_start > 0: mhw_relSeas_start = 0.5*(mhw_relSeas[0] + temp[tt_start-1] - clim['seas'][tt_start-1]) mhw['rate_onset'].append((mhw_relSeas[tt_peak] - mhw_relSeas_start) / (tt_peak+0.5)) else: # MHW starts at beginning of time series if tt_peak == 0: # Peak is also at begining of time series, assume onset time = 1 day mhw['rate_onset'].append((mhw_relSeas[tt_peak] - mhw_relSeas[0]) / 1.) else: mhw['rate_onset'].append((mhw_relSeas[tt_peak] - mhw_relSeas[0]) / tt_peak) if tt_end < T-1: mhw_relSeas_end = 0.5*(mhw_relSeas[-1] + temp[tt_end+1] - clim['seas'][tt_end+1]) mhw['rate_decline'].append((mhw_relSeas[tt_peak] - mhw_relSeas_end) / (tt_end-tt_start-tt_peak+0.5)) else: # MHW finishes at end of time series if tt_peak == T-1: # Peak is also at end of time series, assume decline time = 1 day mhw['rate_decline'].append((mhw_relSeas[tt_peak] - mhw_relSeas[-1]) / 1.) else: mhw['rate_decline'].append((mhw_relSeas[tt_peak] - mhw_relSeas[-1]) / (tt_end-tt_start-tt_peak)) # Flip climatology and intensties in case of cold spell detection if coldSpells: clim['seas'] = -1.*clim['seas'] clim['thresh'] = -1.*clim['thresh'] for ev in range(len(mhw['intensity_max'])): mhw['intensity_max'][ev] = -1.*mhw['intensity_max'][ev] mhw['intensity_mean'][ev] = -1.*mhw['intensity_mean'][ev] mhw['intensity_cumulative'][ev] = -1.*mhw['intensity_cumulative'][ev] mhw['intensity_max_relThresh'][ev] = -1.*mhw['intensity_max_relThresh'][ev] mhw['intensity_mean_relThresh'][ev] = -1.*mhw['intensity_mean_relThresh'][ev] mhw['intensity_cumulative_relThresh'][ev] = -1.*mhw['intensity_cumulative_relThresh'][ev] mhw['intensity_max_abs'][ev] = -1.*mhw['intensity_max_abs'][ev] mhw['intensity_mean_abs'][ev] = -1.*mhw['intensity_mean_abs'][ev] mhw['intensity_cumulative_abs'][ev] = -1.*mhw['intensity_cumulative_abs'][ev] return mhw, clim def blockAverage(t, mhw, clim=None, blockLength=1, removeMissing=False, temp=None): ''' 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_max_max' Maximum 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] 'total_days' Total number of MHW days in each block [days] 'total_icum' Total cumulative intensity over all MHWs in each block [deg. C x 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 'intensity_max_abs', 'intensity_mean_abs', 'intensity_var_abs', and 'intensity_cumulative_abs' are as above except as absolute magnitudes rather than relative to the seasonal climatology or threshold 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]) removeMissing Boolean switch indicating whether to remove (set = NaN) statistics for any blocks in which there were missing temperature values (DEFAULT = FALSE) clim The temperature climatology (including missing value information) as output by marineHeatWaves.detect (required if removeMissing = TRUE) temp Temperature time series. If included mhwBlock will output block averages of mean, max, and min temperature (DEFAULT = NONE) If both clim and temp are provided, this will output annual counts of moderate, strong, severe, and extreme days. Notes: This function assumes that the input time vector consists of continuous daily values. Note that in the case of time ranges which start and end part-way through the calendar year, the block averages at the endpoints, for which there is less than a block length of data, will need to be interpreted with care. Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Feb-Mar 2015 ''' # # Time and dates vectors, and calculate block timing # # Generate vectors for year, month, day-of-month, and day-of-year T = len(t) year = np.zeros((T)) month = np.zeros((T)) day = np.zeros((T)) for i in range(T): year[i] = date.fromordinal(t[i]).year month[i] = date.fromordinal(t[i]).month day[i] = date.fromordinal(t[i]).day # Number of blocks, round up to include partial blocks at end years = np.unique(year) nBlocks = np.ceil((years.max() - years.min() + 1) / blockLength).astype(int) # # Temperature time series included? # sw_temp = None sw_cats = None if temp is not None: sw_temp = True if clim is not None: sw_cats = True else: sw_cats = False else: sw_temp = False # # Initialize MHW output variable # mhwBlock = {} mhwBlock['count'] = np.zeros(nBlocks) mhwBlock['count'] = np.zeros(nBlocks) mhwBlock['duration'] = np.zeros(nBlocks) mhwBlock['intensity_max'] = np.zeros(nBlocks) mhwBlock['intensity_max_max'] = np.zeros(nBlocks) mhwBlock['intensity_mean'] = np.zeros(nBlocks) mhwBlock['intensity_cumulative'] = np.zeros(nBlocks) mhwBlock['intensity_var'] = np.zeros(nBlocks) mhwBlock['intensity_max_relThresh'] = np.zeros(nBlocks) mhwBlock['intensity_mean_relThresh'] = np.zeros(nBlocks) mhwBlock['intensity_cumulative_relThresh'] = np.zeros(nBlocks) mhwBlock['intensity_var_relThresh'] = np.zeros(nBlocks) mhwBlock['intensity_max_abs'] = np.zeros(nBlocks) mhwBlock['intensity_mean_abs'] = np.zeros(nBlocks) mhwBlock['intensity_cumulative_abs'] = np.zeros(nBlocks) mhwBlock['intensity_var_abs'] = np.zeros(nBlocks) mhwBlock['rate_onset'] = np.zeros(nBlocks) mhwBlock['rate_decline'] = np.zeros(nBlocks) mhwBlock['total_days'] = np.zeros(nBlocks) mhwBlock['total_icum'] = np.zeros(nBlocks) if sw_temp: mhwBlock['temp_mean'] = np.zeros(nBlocks) mhwBlock['temp_max'] = np.zeros(nBlocks) mhwBlock['temp_min'] = np.zeros(nBlocks) # Calculate category days if sw_cats: mhwBlock['moderate_days'] = np.zeros(nBlocks) mhwBlock['strong_days'] = np.zeros(nBlocks) mhwBlock['severe_days'] = np.zeros(nBlocks) mhwBlock['extreme_days'] = np.zeros(nBlocks) cats = np.floor(1 + (temp - clim['thresh']) / (clim['thresh'] - clim['seas'])) mhwIndex = np.zeros(t.shape) for ev in range(mhw['n_events']): mhwIndex[mhw['index_start'][ev]:mhw['index_end'][ev]+1] = 1. # Start, end, and centre years for all blocks mhwBlock['years_start'] = years[range(0, len(years), blockLength)] mhwBlock['years_end'] = mhwBlock['years_start'] + blockLength - 1 mhwBlock['years_centre'] = 0.5*(mhwBlock['years_start'] + mhwBlock['years_end']) # # Calculate block averages # for i in range(mhw['n_events']): # Block index for year of each MHW (MHW year defined by start year) iBlock = np.where((mhwBlock['years_start'] <= mhw['date_start'][i].year) * (mhwBlock['years_end'] >= mhw['date_start'][i].year))[0][0] # Add MHW properties to block count mhwBlock['count'][iBlock] += 1 mhwBlock['duration'][iBlock] += mhw['duration'][i] mhwBlock['intensity_max'][iBlock] += mhw['intensity_max'][i] mhwBlock['intensity_max_max'][iBlock] = np.max([mhwBlock['intensity_max_max'][iBlock], mhw['intensity_max'][i]]) mhwBlock['intensity_mean'][iBlock] += mhw['intensity_mean'][i] mhwBlock['intensity_cumulative'][iBlock] += mhw['intensity_cumulative'][i] mhwBlock['intensity_var'][iBlock] += mhw['intensity_var'][i] mhwBlock['intensity_max_relThresh'][iBlock] += mhw['intensity_max_relThresh'][i] mhwBlock['intensity_mean_relThresh'][iBlock] += mhw['intensity_mean_relThresh'][i] mhwBlock['intensity_cumulative_relThresh'][iBlock] += mhw['intensity_cumulative_relThresh'][i] mhwBlock['intensity_var_relThresh'][iBlock] += mhw['intensity_var_relThresh'][i] mhwBlock['intensity_max_abs'][iBlock] += mhw['intensity_max_abs'][i] mhwBlock['intensity_mean_abs'][iBlock] += mhw['intensity_mean_abs'][i] mhwBlock['intensity_cumulative_abs'][iBlock] += mhw['intensity_cumulative_abs'][i] mhwBlock['intensity_var_abs'][iBlock] += mhw['intensity_var_abs'][i] mhwBlock['rate_onset'][iBlock] += mhw['rate_onset'][i] mhwBlock['rate_decline'][iBlock] += mhw['rate_decline'][i] if mhw['date_start'][i].year == mhw['date_end'][i].year: # MHW in single year mhwBlock['total_days'][iBlock] += mhw['duration'][i] else: # MHW spans multiple years year_mhw = year[mhw['index_start'][i]:mhw['index_end'][i]+1] for yr_mhw in np.unique(year_mhw): iBlock = np.where((mhwBlock['years_start'] <= yr_mhw) * (mhwBlock['years_end'] >= yr_mhw))[0][0] mhwBlock['total_days'][iBlock] += np.sum(year_mhw == yr_mhw) # NOTE: icum for a MHW is assigned to its start year, even if it spans mult. years mhwBlock['total_icum'][iBlock] += mhw['intensity_cumulative'][i] # Calculation of category days if sw_cats: for i in range(int(nBlocks)): mhwBlock['moderate_days'][i] = ((year >= mhwBlock['years_start'][i]) * (year <= mhwBlock['years_end'][i]) * mhwIndex * (cats == 1)).astype(int).sum() mhwBlock['strong_days'][i] = ((year >= mhwBlock['years_start'][i]) * (year <= mhwBlock['years_end'][i]) * mhwIndex * (cats == 2)).astype(int).sum() mhwBlock['severe_days'][i] = ((year >= mhwBlock['years_start'][i]) * (year <= mhwBlock['years_end'][i]) * mhwIndex * (cats == 3)).astype(int).sum() mhwBlock['extreme_days'][i] = ((year >= mhwBlock['years_start'][i]) * (year <= mhwBlock['years_end'][i]) * mhwIndex * (cats >= 4)).astype(int).sum() # Calculate averages count = 1.*mhwBlock['count'] count[count==0] = np.nan mhwBlock['duration'] = mhwBlock['duration'] / count mhwBlock['intensity_max'] = mhwBlock['intensity_max'] / count mhwBlock['intensity_mean'] = mhwBlock['intensity_mean'] / count mhwBlock['intensity_cumulative'] = mhwBlock['intensity_cumulative'] / count mhwBlock['intensity_var'] = mhwBlock['intensity_var'] / count mhwBlock['intensity_max_relThresh'] = mhwBlock['intensity_max_relThresh'] / count mhwBlock['intensity_mean_relThresh'] = mhwBlock['intensity_mean_relThresh'] / count mhwBlock['intensity_cumulative_relThresh'] = mhwBlock['intensity_cumulative_relThresh'] / count mhwBlock['intensity_var_relThresh'] = mhwBlock['intensity_var_relThresh'] / count mhwBlock['intensity_max_abs'] = mhwBlock['intensity_max_abs'] / count mhwBlock['intensity_mean_abs'] = mhwBlock['intensity_mean_abs'] / count mhwBlock['intensity_cumulative_abs'] = mhwBlock['intensity_cumulative_abs'] / count mhwBlock['intensity_var_abs'] = mhwBlock['intensity_var_abs'] / count mhwBlock['rate_onset'] = mhwBlock['rate_onset'] / count mhwBlock['rate_decline'] = mhwBlock['rate_decline'] / count # Replace empty years in intensity_max_max mhwBlock['intensity_max_max'][np.isnan(mhwBlock['intensity_max'])] = np.nan # Temperature series if sw_temp: for i in range(int(nBlocks)): tt = (year >= mhwBlock['years_start'][i]) * (year <= mhwBlock['years_end'][i]) mhwBlock['temp_mean'][i] = np.nanmean(temp[tt]) mhwBlock['temp_max'][i] = np.nanmax(temp[tt]) mhwBlock['temp_min'][i] = np.nanmin(temp[tt]) # # Remove years with missing values # if removeMissing: missingYears = np.unique(year[np.where(clim['missing'])[0]]) for y in range(len(missingYears)): iMissing = np.where((mhwBlock['years_start'] <= missingYears[y]) * (mhwBlock['years_end'] >= missingYears[y]))[0][0] mhwBlock['count'][iMissing] = np.nan mhwBlock['duration'][iMissing] = np.nan mhwBlock['intensity_max'][iMissing] = np.nan mhwBlock['intensity_max_max'][iMissing] = np.nan mhwBlock['intensity_mean'][iMissing] = np.nan mhwBlock['intensity_cumulative'][iMissing] = np.nan mhwBlock['intensity_var'][iMissing] = np.nan mhwBlock['intensity_max_relThresh'][iMissing] = np.nan mhwBlock['intensity_mean_relThresh'][iMissing] = np.nan mhwBlock['intensity_cumulative_relThresh'][iMissing] = np.nan mhwBlock['intensity_var_relThresh'][iMissing] = np.nan mhwBlock['intensity_max_abs'][iMissing] = np.nan mhwBlock['intensity_mean_abs'][iMissing] = np.nan mhwBlock['intensity_cumulative_abs'][iMissing] = np.nan mhwBlock['intensity_var_abs'][iMissing] = np.nan mhwBlock['rate_onset'][iMissing] = np.nan mhwBlock['rate_decline'][iMissing] = np.nan mhwBlock['total_days'][iMissing] = np.nan if sw_cats: mhwBlock['moderate_days'][iMissing] = np.nan mhwBlock['strong_days'][iMissing] = np.nan mhwBlock['severe_days'][iMissing] = np.nan mhwBlock['extreme_days'][iMissing] = np.nan mhwBlock['total_icum'][iMissing] = np.nan return mhwBlock def meanTrend(mhwBlock, alpha=0.05): ''' 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 alpha Significance level for estimate of confidence limits on trend, e.g., alpha = 0.05 for 5% significance (or 95% confidence) (DEFAULT = 0.05) Outputs: mean Mean of all MHW properties over all block-averaged values trend Linear trend of all MHW properties over all block-averaged values dtrend One-sided width of (1-alpha)% confidence intevfal on linear trend, i.e., trend lies within (trend-dtrend, trend+dtrend) with specified level of confidence. 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 'intensity_max_abs', 'intensity_mean_abs', 'intensity_var_abs', and 'intensity_cumulative_abs' are as above except as absolute magnitudes rather than relative to the seasonal climatology or threshold 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 2015 ''' # Initialize mean and trend dictionaries mean = {} trend = {} dtrend = {} # Construct matrix of predictors, first column is all ones to estimate the mean, # second column is the time vector, equal to zero at mid-point. t = mhwBlock['years_centre'] X = np.array([np.ones(t.shape), t-t.mean()]).T # Loop over all keys in mhwBlock for key in mhwBlock.keys(): # Skip time-vector keys of mhwBlock if (key == 'years_centre') + (key == 'years_end') + (key == 'years_start'): continue # Predictand (MHW property of interest) y = mhwBlock[key] valid = ~np.isnan(y) # non-NaN indices # Perform linear regression over valid indices if np.isinf(nonans(y).sum()): # If contains Inf values beta = [np.nan, np.nan] elif np.sum(~np.isnan(y)) > 0: # If at least one non-NaN value beta = linalg.lstsq(X[valid,:], y[valid])[0] else: beta = [np.nan, np.nan] # Insert regression coefficients into mean and trend dictionaries mean[key] = beta[0] trend[key] = beta[1] # Confidence limits on trend yhat = np.sum(beta*X, axis=1) t_stat = stats.t.isf(alpha/2, len(t[valid])-2) s = np.sqrt(np.sum((y[valid] - yhat[valid])**2) / (len(t[valid])-2)) Sxx = np.sum(X[valid,1]**2) - (np.sum(X[valid,1])**2)/len(t[valid]) # np.var(X, axis=1)[1] dbeta1 = t_stat * s / np.sqrt(Sxx) dtrend[key] = dbeta1 # Return mean, trend return mean, trend, dtrend def rank(t, mhw): ''' Calculate the rank and return periods of marine heatwaves (MHWs) according to each metric. 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: rank The rank of each MHW according to each MHW property. A rank of 1 is the largest, 2 is the 2nd largest, etc. Each key (listed below) is a list of length N where N is the number of MHWs. returnPeriod The return period (in years) of each MHW according to each MHW property. The return period signifies, statistically, the recurrence interval for an event at least as large/long as the event in quetion. Each key (listed below) is a list of length N where N is the number of MHWs. '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] 'total_days' Total number of MHW days in each block [days] 'total_icum' Total cumulative intensity over all MHWs in each block [deg. C x 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 'intensity_max_abs', 'intensity_mean_abs', 'intensity_var_abs', and 'intensity_cumulative_abs' are as above except as absolute magnitudes rather than relative to the seasonal climatology or threshold Notes: This function assumes that the MHWs were calculated over a suitably long record that return periods make sense. If the record length is a few years or less than this becomes meaningless. Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Sep 2015 ''' # Initialize rank and return period dictionaries rank = {} returnPeriod = {} # Number of years on record nYears = len(t)/365.25 # Loop over all keys in mhw for key in mhw.keys(): # Skip irrelevant keys of mhw, only calculate rank/returns for MHW properties if (key == 'date_end') + (key == 'date_peak') + (key == 'date_start') + (key == 'date_end') + (key == 'date_peak') + (key == 'date_start') + (key == 'index_end') + (key == 'index_peak') + (key == 'index_start') + (key == 'n_events'): continue # Calculate ranks rank[key] = mhw['n_events'] - np.array(mhw[key]).argsort().argsort() # Calculate return period as (# years on record + 1) / (# of occurrences of event) # Return period is for events of at least the event magnitude/duration returnPeriod[key] = (nYears + 1) / rank[key] # Return rank, return return rank, returnPeriod def runavg(ts, w): ''' Performs a running average of an input time series using uniform window of width w. This function assumes that the input time series is periodic. Inputs: ts Time series [1D numpy array] w Integer length (must be odd) of running average window Outputs: ts_smooth Smoothed time series Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Feb-Mar 2015 ''' # Original length of ts N = len(ts) # make ts three-fold periodic ts = np.append(ts, np.append(ts, ts)) # smooth by convolution with a window of equal weights ts_smooth = np.convolve(ts, np.ones(w)/w, mode='same') # Only output central section, of length equal to the original length of ts ts = ts_smooth[N:2*N] return ts def pad(data, maxPadLength=False): ''' Linearly interpolate over missing data (NaNs) in a time series. Inputs: data Time series [1D numpy array] maxPadLength Specifies the maximum length over which to interpolate, i.e., any consecutive blocks of NaNs with length greater than maxPadLength will be left as NaN. Set as an integer. maxPadLength=False (default) interpolates over all NaNs. Written by Eric Oliver, Institue for Marine and Antarctic Studies, University of Tasmania, Jun 2015 ''' data_padded = data.copy() bad_indexes = np.isnan(data) good_indexes = np.logical_not(bad_indexes) good_data = data[good_indexes] interpolated = np.interp(bad_indexes.nonzero()[0], good_indexes.nonzero()[0], good_data) data_padded[bad_indexes] = interpolated if maxPadLength: blocks, n_blocks = ndimage.label(np.isnan(data)) for bl in range(1, n_blocks+1): if (blocks==bl).sum() > maxPadLength: data_padded[blocks==bl] = np.nan return data_padded def nonans(array): ''' Return input array [1D numpy array] with all nan values removed ''' return array[~np.isnan(array)]