Statistical Model Validation#

For the validation of the model, we split the time into a calibration period: (1979-2015) and validation period: (2016-2021)
In the validation period, the results of the model can be compared to the real spectra, providing insights on the performance of the statistical model.

Important

The predictor and predictand classification developed in this section is only used for validating the statistical model. As this section is the repetition of the process presented in previous sections, here only some of the figures are shown for simplicity

# common
import warnings
warnings.filterwarnings('ignore')

import os
import os.path as op
import sys

# pip
import pickle as pk
import numpy as np
import xarray as xr
from sklearn.decomposition import PCA

# DEV: bluemath
sys.path.insert(0, op.join(op.abspath(''), '..', '..', '..', '..', '..'))

# bluemath modules 
from bluemath.binwaves.spectra import set_t_h, group_bins_dir_t
from bluemath.binwaves.plotting.classification import Plot_spectra_hs_pcs, Plot_kmeans_clusters_hs
from bluemath.binwaves.plotting.classification import Plot_specs_hs_inside_cluster, Plot_kmeans_clusters_probability

from bluemath.pca import dynamic_estela_predictor, PCA_EstelaPred, standardise_predictor
from bluemath.kma import kmeans_clustering_pcs, calculate_conditioned_probabilities

from bluemath.plotting.slp import Plot_slp_field, Plot_slp_kmeans
from bluemath.plotting.teslakit import Plot_DWTs_Probs, Plot_Probs_WT_WT
from bluemath.plotting.validation import Plot_validation_spec, Plot_daily_mean, Plot_entropy, calculate_entropy, Plot_prob_comparison

Database and site parameters#

# database
p_data = r'/media/administrador/HD2/SamoaTonga/data'
site = 'Samoa'
p_site = op.join(p_data, site)

# SLP
p_slp = op.join(p_data, 'SLP_hind_daily_wgrad.nc')

# deliverable folder
p_deliv = op.join(p_site, 'd04_seasonal_forecast_swells')

# spectra data
p_spec = op.join(p_deliv, 'spec')

# superpoint
p_superpoint = op.join(p_spec, 'super_point_superposition.nc')

# daily spectrum
p_sp_trans = op.join(p_spec, 'spec_daily_transformed_hs_t.nc')

# estela
p_est = op.join(p_deliv, 'estela')
p_estela = op.join(p_est, 'estela.nc')


# validation files:
p_valid = op.join(p_deliv, 'validation')
p_valid_est = op.join(p_valid, 'estela')

# dynamic predictor
p_dynp = op.join(p_valid_est, 'dynamic_predictor.nc')

# PCA 
p_pca_slp_fit = op.join(p_valid_est, 'PCA_slp_fit.pkl')
p_pca_dyn = op.join(p_valid_est, 'dynamic_predictor_PCA.nc')

p_pca_spec_fit = op.join(p_valid_est, 'PCA_spec_fit.pkl')

# KMA
p_valid_kma = op.join(p_valid, 'kma')

p_kma_slp_model = op.join(p_valid_kma, 'slp_KMA_model.pkl')                     # KMeans model SLP
p_kma_slp_classification = op.join(p_valid_kma, 'slp_KMA_classification.nc')    # classified SLP

p_kma_spec_model = op.join(p_valid_kma, 'spec_KMA_model.pkl')                   # KMeans model spectra
p_kma_spec_classification = op.join(p_valid_kma, 'spec_KMA_classification.nc')  # classified superpoint

# conditioned probabilities
p_probs_cond = op.join(p_valid, 'conditioned_prob_SLP_SPEC.nc')

# generate validation data folders
for p in [p_valid, p_valid_est, p_valid_kma]:
    if not op.isdir(p):
        os.makedirs(p)
        

# validation training time period
time_train_1 = '1979-01-01'
time_train_2 = '2015-12-31'

# validation test time period
time_test_1 = '2016-01-01'
time_test_2 = '2020-12-31'

# K-Means classification parameters - SLP
num_clusters_slp = 81
min_data_slp = 5

# K-Means classification parameters
num_clusters_spec = 49
min_data_frac_spec = 10  # to calculate minimum number of data for each cluster

Conditioned probabilities (1979-2015)#

The figure below represents the conditioned probability of the 49 predictand clusters, conditioned to the 81 predictor clusters.
We observe large variability on the probability distributions of spectral clusters associated to each SLP cluster.

c, i_sp, i_slp = np.intersect1d(sp_mod_daily.time, slp.time, return_indices=True)

slp = slp.isel(time = i_slp)
sp_mod_daily = sp_mod_daily.isel(time = i_sp)

num_clusters_slp = len(np.unique(slp.bmus))
num_clusters_sp = len(np.unique(sp_mod_daily.bmus))

f3 = Plot_Probs_WT_WT(
    slp.bmus, sp_mod_daily.bmus, 
    num_clusters_slp, num_clusters_sp,
    wt_colors = False, 
    ttl = 'Spec bmus conditioned to SLP',
    figsize = [15,15],
    vmax = 0.07,
);

# calculate conditioned probabilities
conditioned_prob = calculate_conditioned_probabilities(slp, num_clusters_slp, sp_mod_daily, num_clusters_sp)

conditioned_prob.to_netcdf(p_probs_cond)

Validation 2016-2021#

The validation period corresponds to the years from 2016-2021, which were not included in the calibration process

Obtain SLP cluster and Spec probabilities for each day#

Here, we are projecting the SLP fields from the validation period into the PCA calibration space to obtain the SLP cluster that characterizes each day

SLP fields for the validation period:

# load the slp with partitions
slp = xr.open_dataset(p_slp)

# select training time period
slp = slp.sel(time = slice(time_test_1, time_test_2))

Constructing the dynamic predictor:

slp_dy = dynamic_estela_predictor(slp, 'slp', estela_D)

Transform slp from the validation period into the calibration PCA space:

# standardise estela predictor
d_pos, xds_norm = standardise_predictor(slp_dy, 'slp')
data_norm = xds_norm.pred_norm.values

# load estela pca and get PCs
pca_fit = pk.load(open(p_pca_slp_fit, 'rb'))
PCs_2020 = pca_fit.transform(data_norm)

# load estela dynamic predictor PCA
slp_pcs = xr.open_dataset(p_pca_dyn)

# get PCs that explain 90% of the EV
#n_percent = np.cumsum((slp_pcs.variance / np.sum(slp_pcs.variance))*100)
#ev=90
#n_pcs=np.where(n_percent>ev)[0][0]
#print('Number of PCs explaining ' + str(ev) + '% of the EV is: ' + str(n_pcs))

PCs_2020 = PCs_2020[:, :n_pcs_slp]

Obtain SLP clusters where each day belongs to:

# load slp kma
slp_kma = xr.open_dataset(p_kma_slp_classification)
num_clusters_slp = len(np.unique(slp_kma.bmus))

mask = np.where(np.isnan(slp_kma.est_mask) == True)
kma_fit = pk.load(open(p_kma_slp_model, "rb"))[0]

# kma predict validation PCs
bmus_slp = kma_fit.predict(PCs_2020)

sorted_bmus = np.zeros((len(bmus_slp),),) * np.nan
for i in range(num_clusters_slp):
    posc = np.where(bmus_slp == slp_kma.kma_order[i].values)
    sorted_bmus[posc] = i

slp_dy['bmus'] = sorted_bmus

Obtain the probability of each spectral cluster at the daily scale:

spec_kma = xr.open_dataset(p_kma_spec_classification)

num_clusters_spec = len(np.unique(spec_kma.bmus))

# Intersect dates
c, i_sp, i_slp = np.intersect1d(spec_kma.time, slp_kma.time, return_indices=True)

slp_kma = slp_kma.isel(time=i_slp)
spec_kma = spec_kma.isel(time=i_sp)

# load conditioned probabilities
cond_prob = xr.open_dataset(p_probs_cond)

p_f = np.full([len(cond_prob.ir), len(cond_prob.ic), len(slp_dy.time)], np.nan)

for a in range(len(slp_dy.time)):
    p_f[:, :, a] = cond_prob.prob_sp[:, :, slp_dy.bmus[a].values.astype('int')]

slp_dy['prob_spec'] = (('ir', 'ic', 'time'), p_f)
slp_dy['bmus'] = ('time', slp_dy.bmus)

Load spectral hindcast for comparison:

# spectral hindcast
sp_mod_daily = xr.open_dataset(p_sp_trans)

Results at the daily scale#

Each day, represented by one KMA cluster, has a specific distribution of different spectra probabilities

Figure Explanation

The figures below, both at the daily and monthly scale, have 3 different panels.
Upper left panel: SLP associated to that day
Lower left panel: Real spectra from the hindcast
Right panel: This panel displays only the clusters that are probable under those SLP conditions according to the model, where the blue line surrounding each spectrum has the information of the probability.

Figure Explanation

Also, the mean hindcast spectra and anomalies are shown below for comparison.
Mean model spectra is calculated as the sum of each cluster spectra multiplied by its probability of occurrence

# select a date to plot
day = '2017-07-10'

# slp for selected date
slp_d = slp_dy.sel(time = day)

# spectra hindcast for selected date
z = 4 * np.sqrt(sp_mod_daily.sel(time = day).h.values)

extent = (50, 290, -67, 65)

Plot_validation_spec(
    slp_d, spec_kma,
    extent = extent,
    real_spec = z,
    prob_min = 0, prob_max = 0.08,
    s_min = 990, s_max = 1030,
    lw = 4,
    pt = [184.84, -21.21],
    min_prob_clust = 0.02,
    title = day,
);


# daily mean spec
num_clus_pick = 49
perct_pick = []

# mean
Plot_daily_mean(
    slp_dy, sp_mod_daily, 
    [day], 
    spec_kma, 
    num_clus_pick = num_clus_pick,
    perct_pick = perct_pick,
);


# anomaly (needs mean_spec)
mean_spec = 4 * np.sqrt(np.mean(sp_mod_daily.h.values, axis=0))

Plot_daily_mean(
    slp_dy, sp_mod_daily, 
    [day], 
    spec_kma, 
    num_clus_pick = num_clus_pick,
    perct_pick = perct_pick, 
    mean_spec = mean_spec, 
    vmax=0.5, vmin=-0.5, 
    cmap='RdBu_r',
);

Results at the monthly scale#

For the purpose of the work, which is to develop a seasonal forecast, the results are aggregated to the monthly scale
The results for two months, one in January, and one in July are shown below

# resample slp and spectra to monthly data, using monthly means

slp_dy_m = slp_dy.resample(time = '1M', skipna = True).mean()

spec_mod_m = sp_mod_daily.resample(time = '1M', skipna = True).mean()

# select a month to plot
month = ['2017-01', '2020-07']

# SLP field and cluster probabilities
min_prob_clust = 0.005

# slp for selected month
slp_m = slp_dy_m.sel(time = month[0]).isel(time = 0)

# spectra hindcast for selected month
z = 4 * np.sqrt(spec_mod_m.sel(time=month[0]).isel(time=0).h.values)

a = Plot_validation_spec(
    slp_m, spec_kma,
    extent = extent,
    real_spec = z,
    prob_min = 0, prob_max=0.08,
    s_min = 990, s_max=1030,
    lw = 4,
    pt = [184.84, -21.21],
    min_prob_clust = 0.03,
    title = month,
);

# Monthly aggregated
num_clus_pick = 49
perct_pick = []

# mean
Plot_daily_mean(
    slp_dy_m, spec_mod_m,
    month,
    spec_kma, 
    num_clus_pick = num_clus_pick,
    perct_pick = perct_pick,
);


# respect to the mean month spectra
mean_spec = 4 * np.sqrt(np.mean(spec_mod_m.h.values, axis=0))

# respect to the mean specific month
#mean_spec = 4 * np.sqrt(np.mean(spec_mod_m.isel(time=np.where(spec_mod_m.time.dt.month==slp_m.time.dt.month.values)[0]).h.values,axis=0))

# anomaly (needs mean_spec)
Plot_daily_mean(
    slp_dy_m, spec_mod_m,
    month,
    spec_kma,
    num_clus_pick = num_clus_pick,
    perct_pick = perct_pick, 
    mean_spec = mean_spec,
    vmax = 0.4, vmin = -0.4,
    cmap = 'RdBu_r',
);

Compare cluster probabilities#

To compare the capabilities of the model on correctly reproducing the ocurrence probability of the different predicant clusters, we have compared the probabilities and calculated the relative entropy (RE) of real versus modelled probabilities

# spectral hindcast
sp_mod_daily = xr.open_dataset(p_sp_trans)

# select training time period
sp_mod_daily_fit = sp_mod_daily.sel(time = slice(time_train_1, time_train_2))

# load training PCA
ipca = pk.load(open(p_pca_spec_fit, 'rb'))

# organize spatial data into lines
m = np.full((len(sp_mod_daily.time), len(sp_mod_daily.dir_bins) * len(sp_mod_daily.t_bins)), 0.000001)

for ti in range(m.shape[0]):
    m[ti, :] = np.ravel(sp_mod_daily.h_t.isel(time = ti).values)
    
# calculate PCs
PCs_val = ipca.transform(m)

# load training kma
kma_fit = pk.load(open(p_kma_spec_model, 'rb'))[0]

# predict bmus
bmus_val = kma_fit.predict(PCs_val[:, :n_pcs_spec])

# sort bmus with kma order
kma_order = spec_kma.kma_order.values

sorted_bmus_val = np.zeros((len(bmus_val),),) * np.nan
for i in range(len(spec_kma.clusters)):
    posc = np.where(bmus_val == kma_order[i])
    sorted_bmus_val[posc] = i

# add sorted bmus to spectra dataset
sp_mod_daily['bmus'] = ('time', sorted_bmus_val.astype('int'))

Here, the comparison of real versus modelled probabilities is shown for a month

month = ['2017-03']

Plot_prob_comparison(sp_mod_daily, slp_dy_m, spec_kma, month);

RE = 0.5743748637250338

# calculate entropy
re_all = calculate_entropy(sp_mod_daily, slp_dy_m, spec_kma)

print('RE mean = ' + str(np.nanmean(re_all[1:])))

RE mean = 0.9411161735077451

Here, the RE for the validation period is shown

Plot_entropy(slp_dy_m, re_all);

Statistical Model Validation

Contents