Producing the ENSO index verification heatmaps

In this example we use the computed El Niño–Southern Oscillation (ENSO) indices from the previous Notebook to produce heatmaps of temporal correlations, as included on the C3S seasonal real-time forecast verification page.

Some information on ENSO impacts in Europe can be found on a page in the C3S documentation.

1. Configuration

2. Load the pre-computed ERA5 indices and prepare to compute correlations

3. Loop over the hindcast start months and compute correlations

4. Plot the correlation heatmaps

Note

Some of the python code uses a set of utility functions available in the GitHub respository which hosts these Notebooks.

Configuration

Import required modules, define the library of C3S systems, and define some parameters related to the indices previously computed.

import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import xarray as xr
import regionmask
import seaborn as sns
from matplotlib.colors import BoundaryNorm
from dateutil import relativedelta
import sys
from pathlib import Path

# Add parent directory to path so we can import the utils package
sys.path.insert(0, str(Path.cwd().parent))
from utils.tools import getForecastsystemDetails

# the data path where the indices were saved
data_path = 'data'

# define the common hindcast period
hc_period = ['1993', '2016']
hc_str = '_'.join([hc_period[0], hc_period[1]])

# select a forecast system
origin = 'ukmo'  # originating centre label for the CDS request
system = '604'  # forecast system identifier for the CDS request
prov = f'{origin}.s{system}'

# Get system details from configuration files in utils directory using format ukmo.s604
system_details = getForecastsystemDetails(prov)
print(system_details)

{'origin_details': {'MARS_origin': 'egrr', 'label': 'Met Office'}, 'system_details': {'name': 'GloSea6 (system=604)', 'isLagged': True, 'ens_size': {'products': {'fc': 50}}}}

Load the pre-computed ERA5 indices and prepare to compute correlations

# load in reanalysis
era5_file_name = '/era_5_nino_ind_1993_2016.nc'
era5_ind = xr.open_dataarray(data_path + era5_file_name)
era5_ind
# create a correlation dataset
corr = xr.Dataset()

Loop over the hindcast start months and compute correlations

Loop over the start months. For each start month do the following:

• deal with the missing Jan 1993 start if the centre is ukmo

• create a copy of the ERA5 array (to rearrange to match the valid times of the start month in question)

• read in the hindcast index for the given start, and refactor to include valid time

• sub-select qualifying ensemble members, if needed (for forecast systems where not all available ensemble members are used in the graphical products)

• refactor the ERA5 array to align with the hindcast leadtimes

• for each leadmonth compute the temporal correlation

• combine the leadmonth correlations and append to the array of correlations by start month

The result is an array of temporal correlations for each start month and leadmonth.

for stmonth in range(1, 13):
    print(f'START MONTH = {stmonth}')
    # Note that if the chosen system does not yet have all starts then this needs to be dealt with here
    # This also includes January 1993 for UKMO meaning the ERA5 data needs to be adjusted to match this
    # Met Office systems missing data - start month
    if origin == 'ukmo' and stmonth == 1:
        # edit config to deal with Jan 1993 data gap in Met Office data
        hc_period[0] = '1994'
        o = era5_ind.sel(time=slice(hc_period[0], hc_period[1])).copy()
    else:
        o = era5_ind.copy()

    # read hindcast data array for current start
    hc_file_name = f'/sst_mm_{origin}_{system}_{hc_str}_st_{stmonth}_nino_ind.nc'
    hcst = xr.open_dataarray(data_path + hc_file_name)

    #hcst = hcst.reset_coords(names="valid_time", drop=True)
    vt = xr.DataArray(dims=('start_date', 'forecastMonth'),
                      coords={'forecastMonth': hcst.forecastMonth, 'start_date': hcst.start_date})
    vt.data = [[pd.to_datetime(std) + relativedelta.relativedelta(months=fcmonth - 1) for fcmonth in vt.forecastMonth.values] for
               std
               in vt.start_date.values]
    hcst = hcst.assign_coords(valid_time=vt)

    # select qualifying ensemble members
    # Check if ensemble size in data exceeds the configured size
    ens_size_hc = system_details['system_details'].get('ens_size', {}).get('products', {}).get('hc',None)
    if ens_size_hc is not None and len(hcst.number) > ens_size_hc:
        hcst = hcst.sel(number=slice(None, ens_size_hc - 1))

    # refactor reanalysis to align with hindcast leadtimes
    o = o.rename({'time': 'start_date'}).swap_dims({'start_date': 'valid_time'})
    fcmonths = [mm + 1 if mm >= 0 else mm + 13 for mm in
                [t.month - stmonth for t in pd.to_datetime(o.valid_time.values)]]
    o = o.assign_coords(forecastMonth=('valid_time', fcmonths))
    o = o.loc[(o['forecastMonth'] <= 6)]
    
    # loop over forecast month, define the hindcast mean over the ensemble members, and compute the correlation over the hindcast years
    l_corr = list()
    l_corr_det = list()
    for this_fcmonth in hcst.forecastMonth.values:
        thishcst = (hcst.sel(forecastMonth=this_fcmonth).swap_dims({'start_date': 'valid_time'}))
        thishcst = thishcst.mean('number')
        thisobs = o.where(o.valid_time == thishcst.valid_time, drop=True)
        l_corr.append(xr.corr(thishcst, thisobs, dim='valid_time'))
        
    thiscorr = xr.concat(l_corr, dim='forecastMonth')
    thiscorr = thiscorr.assign_coords(forecastMonth=[1, 2, 3, 4, 5, 6])
    thiscorr = thiscorr.expand_dims(dim={'start_month': [stmonth]}, axis=0)
    corr = xr.merge([corr, thiscorr], compat='no_conflicts', join='outer')

Save the result, and inspect it.

file_name = f'/sst_mm_{origin}_{system}_{hc_str}_nino_ind_corr.nc'
corr.to_netcdf(data_path + file_name, mode='w')
corr

<xarray.Dataset> Size: 3kB
Dimensions:        (start_month: 12, region: 4, forecastMonth: 6)
Coordinates:
  * start_month    (start_month) int64 96B 1 2 3 4 5 6 7 8 9 10 11 12
  * region         (region) int64 32B 0 1 2 3
  * forecastMonth  (forecastMonth) int64 48B 1 2 3 4 5 6
    surface        float64 8B 0.0
    abbrevs        (region) <U3 48B '1+2' '3' '3.4' '4'
    names          (region) <U7 112B 'NINO1+2' 'NINO3' 'NINO3.4' 'NINO4'
    number         int64 8B 0
    step           timedelta64[ns] 8B 00:00:00
Data variables:
    sst            (start_month, forecastMonth, region) float64 2kB 0.9424 .....

xarray.Dataset

• Dimensions:
  • start_month: 12
  • region: 4
  • forecastMonth: 6
• Coordinates: (8)
  • start_month
    (start_month)
    int64
    1 2 3 4 5 6 7 8 9 10 11 12

    array([ 1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12])

  • region
    (region)
    int64
    0 1 2 3

    array([0, 1, 2, 3])

  • forecastMonth
    (forecastMonth)
    int64
    1 2 3 4 5 6

    array([1, 2, 3, 4, 5, 6])

  • surface
    ()
    float64
    0.0

    array(0.)

  • abbrevs
    (region)
    <U3
    '1+2' '3' '3.4' '4'

    array(['1+2', '3', '3.4', '4'], dtype='<U3')

  • names
    (region)
    <U7
    'NINO1+2' 'NINO3' 'NINO3.4' 'NINO4'

    array(['NINO1+2', 'NINO3', 'NINO3.4', 'NINO4'], dtype='<U7')

  • number
    ()
    int64
    0

    array(0)

  • step
    ()
    timedelta64[ns]
    00:00:00

    array(0, dtype='timedelta64[ns]')

• Data variables: (1)
  • sst
    (start_month, forecastMonth, region)
    float64
    0.9424 0.9897 ... 0.6788 0.8228

    array([[[0.94240909, 0.98972404, 0.99513449, 0.98433147],
            [0.65706559, 0.94456101, 0.97790239, 0.96814965],
            [0.3502185 , 0.76739183, 0.93429987, 0.95757185],
            [0.41117435, 0.6784339 , 0.94070999, 0.91308116],
            [0.69327205, 0.66371638, 0.76225142, 0.86683375],
            [0.71094977, 0.58482052, 0.63686247, 0.77263278]],
    
           [[0.84439034, 0.98039013, 0.9921852 , 0.98571335],
            [0.63437098, 0.83999722, 0.96169238, 0.98473573],
            [0.69523474, 0.7973263 , 0.95832987, 0.96960242],
            [0.79983336, 0.75746209, 0.83300311, 0.92077784],
            [0.79109931, 0.69715437, 0.68607617, 0.79037122],
            [0.70465558, 0.72059784, 0.69652643, 0.61608645]],
    
           [[0.84579243, 0.93653245, 0.97877778, 0.98839766],
            [0.82822748, 0.86311771, 0.96285862, 0.96842286],
            [0.83363579, 0.81567227, 0.89106493, 0.95272668],
            [0.80203682, 0.75730489, 0.78079115, 0.85424148],
            [0.70593984, 0.73151795, 0.75490607, 0.73192892],
            [0.71320212, 0.74761858, 0.76314281, 0.63595066]],
    ...
           [[0.97805953, 0.99650971, 0.99249697, 0.98231661],
            [0.97641444, 0.98692731, 0.98342466, 0.98762761],
            [0.96227185, 0.97686991, 0.97559375, 0.98005345],
            [0.85054648, 0.96312176, 0.96768634, 0.94825366],
            [0.56113528, 0.89069524, 0.94063967, 0.92528301],
            [0.27714915, 0.81224505, 0.92947953, 0.91472918]],
    
           [[0.98067472, 0.99470999, 0.99604183, 0.99570169],
            [0.95538153, 0.98364869, 0.9834228 , 0.98429683],
            [0.79658272, 0.95780655, 0.97128157, 0.95772124],
            [0.57149912, 0.90365301, 0.95058988, 0.93414128],
            [0.32978678, 0.8080699 , 0.91797313, 0.91215383],
            [0.38524437, 0.66461042, 0.9207957 , 0.89619838]],
    
           [[0.9732305 , 0.99369366, 0.99454826, 0.99267464],
            [0.80495406, 0.97119806, 0.98439414, 0.96698167],
            [0.50778328, 0.91225588, 0.96005467, 0.95275177],
            [0.32848753, 0.80908871, 0.93519551, 0.92937432],
            [0.42555509, 0.68735956, 0.92444252, 0.89963223],
            [0.57616062, 0.513548  , 0.67879785, 0.82275277]]])

Plot the correlation heatmaps

Using the computed correlations, construct correlation heatmaps for each NINO index computed in the previous Notebook.

Define labels (valid months or start months, origin and system) and load the correlations.

import calendar
mon_labels = [m.lower() for m in calendar.month_abbr[1:]]
orig_label = system_details.get('origin_details', {}).get('label', 'UnknownOrigin')
sys_name = system_details.get('system_details',{}).get('name','UnknownSystem')

# load the correlations
file_name = f'/sst_mm_{origin}_{system}_{hc_str}_nino_ind_corr.nc'
corr = xr.open_dataarray(data_path + file_name)

Loop over the regions for which the indices and correlations were computed:

• refactor the correlations for each region/index in a pandas dataframe (remove unnecessary dimensions and variables)

• plot correlation as a function of start month and leadtime using a seaborn heatmap

• refactor the correlations again to include valid month instead of start month

• plot correlation as a function of valid month and leadtime

# make the plot - loop over nino regions
for ind in corr.region.values:
    print(ind, corr.sel(region=ind).names.values)
    ind_name = str(corr.sel(region=ind).names.values)
    hm_cm = plt.get_cmap('RdYlBu_r')
    hm_levs = np.linspace(0, 1, 11)
    norm = BoundaryNorm(hm_levs, ncolors=hm_cm.N, clip=True)

    # HEATMAP (correlation as a function of start_month and leadtime)
    df = (corr.sel(region=ind).to_dataframe()
          .drop(['region', 'surface', 'abbrevs', 'step', 'names', 'number'], axis=1))
    df = df.unstack().droplevel(0, axis=1)
    fig = plt.figure(figsize=(10, 10))
    ax = sns.heatmap(df.T, square=True, cmap=hm_cm, vmin=hm_levs.min(), vmax=hm_levs.max(), norm=norm,
                     annot=True, fmt=".2f", cbar=False, linewidths=0.5, linecolor='white')
    ax.grid(False)
    ax.invert_yaxis()

    tit_txt1 = f'{ind_name} index (temporal correlation with ERA5) '
    tit_txt2 = f'{orig_label} {sys_name} '
    
    plt.title(tit_txt2 + tit_txt1)
    plt.xlabel('start month')
    plt.ylabel('leadtime')
    # Set ticks labels for x-axis
    ax.set_xticklabels(mon_labels)  
    ax.set_yticklabels(['1', '2', '3', '4', '5', '6'])
    plt.show()

    # HEATMAP (correlation as a function of valid_month and leadtime)
    df_vt = df.stack(future_stack=True).reset_index()
    df_vt['valid_month'] = df_vt.apply(lambda row: row['start_month']
                                       + row['forecastMonth'] - 1, axis=1).astype('int32')
    df_vt['valid_month'] = [int(mm) if mm < 13 else mm - 12 for mm in df_vt['valid_month']]
    df_vt = df_vt.drop('start_month', axis=1).set_index(
        ['valid_month', 'forecastMonth']).unstack().transpose().droplevel(0)
    fig = plt.figure(figsize=(10, 10))
    ax = sns.heatmap(df_vt, square=True, cmap=hm_cm, vmin=hm_levs.min(), vmax=hm_levs.max(), norm=norm,
                     annot=True, fmt=".2f", cbar=False, linewidths=0.5, linecolor='white')
    ax.grid(False)
    ax.invert_yaxis()
    plt.title('\n'.join([tit_txt2 + tit_txt1, '']))
    plt.xlabel('valid month')
    plt.ylabel('leadtime')
    # Set ticks labels for x-axis
    ax.set_xticklabels(mon_labels)  # , rotation='vertical', fontsize=18)
    plt.show()

0 NINO1+2

1 NINO3

2 NINO3.4

3 NINO4

Extensions to these single-system monthly correlations heatmaps would be to compute indices for multiple systems and then combine them into multi-system combinations. The approach to variance-normalisation in the C3S SST index charts is outlined on the products description page, and in the SEAS5 user guide.