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 verification page.

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

Configuration

Import required modules 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 dateutils import relativedelta

# dictionary of systems, versions and attributes
providers = {
    'ecmf.s51': {'cds_name': 'ecmwf', 'plot_name': 'ECMWF', 'plot_system': 'SEAS5', 'cds_system': '51', 'hcsize': 25,
                 'fcsize': 51, 'lagged': False},
    'lfpw.s8': {'cds_name': 'meteo_france', 'plot_name': 'Météo-France', 'plot_system': 'System 8', 'cds_system': '8',
                'hcsize': 25, 'fcsize': 51, 'lagged': False},
    'egrr.s602': {'cds_name': 'ukmo', 'plot_name': 'Met Office', 'plot_system': 'GloSea6', 'cds_system': '602',
                  'hcsize': 28, 'fcsize': 60, 'lagged': True},
    # 'egrr.s603': {'cds_name': 'ukmo', 'plot_name': 'Met Office', 'plot_system': 'GloSea6', 'cds_system': '603',
    #               'hcsize': 28, 'fcsize': 60, 'lagged': True},
    'edzw.s21': {'cds_name': 'dwd', 'plot_name': 'DWD', 'plot_system': 'GCFS2.1', 'cds_system': '21', 'hcsize': 30,
                 'lagged': False},
    'cmcc.s35': {'cds_name': 'cmcc', 'plot_name': 'CMCC', 'plot_system': 'SPS3.5', 'cds_system': '35', 'hcsize': 40,
                 'fcsize': 50, 'lagged': False},
   'kwbc.s2': {'cds_name': 'ncep', 'plot_name': 'NCEP', 'plot_system': 'CFSv2', 'cds_system': '2', 'hcsize': 20,
               'lagged': True},
    'rjtd.s3': {'cds_name': 'jma', 'plot_name': 'JMA', 'plot_system': 'CPS3', 'cds_system': '3', 'hcsize': 10,
                'fcsize': 150, 'lagged': True},
    'cwao.s2': {'cds_name': 'eccc', 'plot_name': 'ECCC', 'plot_system': 'CanCM4i', 'cds_system': '2', 'hcsize': 10,
                'fcsize': 10, 'lagged': False},
    'cwao.s3': {'cds_name': 'eccc', 'plot_name': 'ECCC', 'plot_system': 'GEM5-NEMO', 'cds_system': '3', 'hcsize': 10,
                'fcsize': 10, 'lagged': False},
    # 'cwao.s4': {'cds_name': 'eccc', 'plot_name': 'ECCC', 'plot_system': 'CanESM5.1p1bc', 'cds_system': '4', 'hcsize': 10,
    #             'fcsize': 10, 'lagged': False},
    # 'cwao.s5': {'cds_name': 'eccc', 'plot_name': 'ECCC', 'plot_system': 'GEM5.2-NEMO', 'cds_system': '5', 'hcsize': 10,
    #             'fcsize': 10, 'lagged': False},
}

# select a system and version
prov = 'lfpw.s8'
centre = providers[prov]['cds_name']
version = providers[prov]['cds_system']
# define the common hindcast period
hc_period = ['1993', '2016']
hc_str = '_'.join([hc_period[0], hc_period[1]])
# the data path where the indices were saved
data_path = '/data'

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 array
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 lagged start systems where not all dates are used in the graphical products)

• refactor the ERA5 array to allign with the hindcast leadtimes

• for each leadmonth compute the temporal correlation

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

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

for stmonth in range(1, 13):
    print('START MONTH = {}'.format(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 centre == '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 HCST DATAFRAME FOR GIVEN START
    hc_file_name = '/{}_mm_{}_{}_{}_st_{}_nino_ind.nc'.format('sst', centre, version, hc_str, stmonth)
    hcst = xr.open_dataarray(data_path + hc_file_name)
    hcst = hcst.rename({'time': 'start_date'})

    #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(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
    if prov == 'kwbc.s2':  # Only case where member sub-selection is needed from the CDS data
        hcst = hcst.sel(number=slice(None, 19))

    # refactor reanalysis to allign 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])

START MONTH = 1
START MONTH = 2
START MONTH = 3
START MONTH = 4
START MONTH = 5
START MONTH = 6
START MONTH = 7
START MONTH = 8
START MONTH = 9
START MONTH = 10
START MONTH = 11
START MONTH = 12

Save the result, and inspect it.

file_name = '/{}_mm_{}_{}_{}_nino_ind_corr.nc'.format('sst', centre, version, hc_str)
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.8843 .....

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.8843 0.9749 ... 0.547 0.8424

    array([[[0.88428325, 0.97486296, 0.98754286, 0.98760322],
            [0.60822438, 0.92039486, 0.95650449, 0.96280133],
            [0.19136528, 0.75826997, 0.91534092, 0.95686024],
            [0.22341532, 0.53335769, 0.8603799 , 0.91650688],
            [0.43586997, 0.44712255, 0.69001741, 0.85897938],
            [0.42898064, 0.48539336, 0.63148798, 0.7681749 ]],
    
           [[0.82171896, 0.94152061, 0.96517206, 0.97698454],
            [0.38277638, 0.80579707, 0.92722522, 0.96092377],
            [0.29266329, 0.63716356, 0.90934797, 0.93186019],
            [0.43387371, 0.44154574, 0.73777226, 0.9031279 ],
            [0.44157277, 0.59686254, 0.78593459, 0.86665536],
            [0.3013047 , 0.52985381, 0.70890726, 0.80963241]],
    
           [[0.86952867, 0.90835389, 0.9640901 , 0.98631497],
            [0.61194553, 0.72242112, 0.94430708, 0.97354831],
            [0.59348479, 0.52977095, 0.76074186, 0.92107745],
            [0.4810353 , 0.64739365, 0.80564077, 0.86959181],
            [0.36748945, 0.61026551, 0.77919187, 0.83636535],
            [0.43502273, 0.6254185 , 0.76445873, 0.76085631]],
    ...
           [[0.95246379, 0.98891461, 0.99319801, 0.97057888],
            [0.93208014, 0.97824024, 0.98414882, 0.97243545],
            [0.90395994, 0.96593881, 0.97044046, 0.96185535],
            [0.81947906, 0.96969538, 0.95983992, 0.92237681],
            [0.54816711, 0.89253996, 0.93440922, 0.91007183],
            [0.18854009, 0.74906627, 0.88863746, 0.9052694 ]],
    
           [[0.94709726, 0.99171239, 0.99365991, 0.98334209],
            [0.91997758, 0.98570181, 0.98515185, 0.97319554],
            [0.77776976, 0.96923018, 0.97782602, 0.9466315 ],
            [0.44821512, 0.89108526, 0.95258034, 0.93355036],
            [0.10712365, 0.79338086, 0.91479902, 0.91786239],
            [0.26010704, 0.53916141, 0.85040554, 0.87164005]],
    
           [[0.9523891 , 0.99577448, 0.99399865, 0.98244374],
            [0.82479727, 0.98747117, 0.99257037, 0.97964003],
            [0.50650936, 0.92990577, 0.97110609, 0.96516766],
            [0.09619206, 0.77348797, 0.92682981, 0.94847549],
            [0.23877698, 0.51166995, 0.84941693, 0.89640995],
            [0.35026087, 0.21924205, 0.54699072, 0.84240341]]])

• Indexes: (3)
  • start_month
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], dtype='int64', name='start_month'))

  • region
    PandasIndex

    PandasIndex(Index([0, 1, 2, 3], dtype='int64', name='region'))

  • forecastMonth
    PandasIndex

    PandasIndex(Index([1, 2, 3, 4, 5, 6], dtype='int64', name='forecastMonth'))

• Attributes: (0)

Plot the correlation heatmaps

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

Define labels for the valid months or start months, and load the correlations.

# month labels for the plots
mon_labels = ['jan', 'feb', 'mar', 'apr', 'may', 'jun', 'jul', 'aug', 'sep', 'oct', 'nov', 'dec']

# load the correlations
file_name = '/{}_mm_{}_{}_{}_nino_ind_corr.nc'.format('sst', centre, version, hc_str)
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,
                     linewidth=0.5,
                     annot=True, fmt=".2f", cbar=False)
    ax.invert_yaxis()

    tit_txt1 = '{} index (temporal correlation with ERA5) '.format(ind_name)
    if centre == 'ukmo':
        tit_txt2 = '{} {} (system={}), '.format(providers[prov]['plot_name'], providers[prov]['plot_system'],
                                                providers[prov]['cds_system'])
    if prov == 'ecmwf.51':
        tit_txt2 = '{} {} (C3Sv{}), '.format(providers[prov]['plot_name'], providers[prov]['plot_system'],
                                             providers[prov]['cds_system'])
    else:
        tit_txt2 = '{} {} , '.format(providers[prov]['plot_name'], providers[prov]['plot_system'])
    plt.title('\n'.join([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,
                     linewidth=0.5,
                     annot=True, fmt=".2f", cbar=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.