Tailor-Made classification of large scale predictor#

In this section, we are going to classify our tailor-made predictor (SLP fields and gradients modified using ESTELA) into a number of clusters (N=81) by reducing the dimensionality of the data by means of a principal component analysis (PCA) and then, performing a clustering of the PCs with the K-Means algorithm.

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

import os
import os.path as op
import sys
import pickle as pk

# pip
import numpy as np
import xarray as xr

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

# bluemath modules 
from bluemath.pca import dynamic_estela_predictor, PCA_EstelaPred
from bluemath.kma import kmeans_clustering_pcs

from bluemath.plotting.slp import Plot_slp_field, Plot_slp_kmeans
from bluemath.plotting.estela import Plot_estela
from bluemath.plotting.pca import Plot_EOFs_EstelaPred
from bluemath.plotting.teslakit import Plot_DWTs_Probs

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')

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


# output files:

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

# PCA 
p_pca_fit = op.join(p_est, 'PCA_fit.pkl')
p_pca_dyn = op.join(p_est, 'dynamic_predictor_PCA.nc')

# KMA
p_kma = op.join(p_deliv, 'kma')

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

# make output folders
for p in [p_kma]:
    if not op.isdir(p):
        os.makedirs(p)
        

# K-Means classification parameters
num_clusters = 81
min_data = 5

SLP data

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

PCA : Principal component analysis#

Area of SLP fields#

First, we define the area of influence, based on the wave generation area extracted from ESTELA
Here is an example for a given time, of the SLP area used for constructing the predictor.

Plot_slp_field(slp, extent = (50, 290, -67, 65), time = 4804);

PCA#

The PCA is used to reduce the high dimensionality of the original data space and thus, simplify the classification process, transforming the SLP predictor fields into a number of spatial modes with an associated temporal variability over time.

PCA projects the original data on a new space searching for the maximum variance of the sample data.

The eigenvectors (the empirical orthogonal functions, EOFs) of the data covariance matrix define the vectors of the new space.
They represent the spatial variability.

# PCA estela predictor
[pca, ipca] = PCA_EstelaPred(pca_dy, 'slp')

# store PCA fit
with open(p_pca_fit, 'wb') as fW:
    pk.dump(ipca, fW)
    
# store dynamic predictor PCA
pca.to_netcdf(p_pca_dyn)

    
# load PCA 
#pca = xr.open_dataset(p_pca_dyn)
#ipca = pk.load(open(p_pca_fit, 'rb'))


# set pred_time as a coordinate
PCs = pca.assign_coords({'time': (('time'), pca.pred_time.values)}).PCs

# select slp at PCA predictor time
slp = slp.sel(time = pca.pred_time)

Here we select the PCs explaining 90% of the variance

ev = 90

n_percent = np.cumsum((pca.variance / np.sum(pca.variance)) * 100)
n_pcs = np.where(n_percent > ev)[0][0]

print('Number of PCs explaining ' + str(ev) + '% of the EV is: ' + str(n_pcs))
print('Number of PCs explaining {0}% of the EV is: {1}'.format(ev, n_pcs))
      

Number of PCs explaining 90% of the EV is: 521
Number of PCs explaining 90% of the EV is: 521

The eigenvectors (the empirical orthogonal functions, EOFs) of the data covariance matrix define the vectors of the new space.
They represent the spatial variability.

# plot first PCA EOFs

n_pcs_plot = 3

Plot_EOFs_EstelaPred(pca, n_pcs_plot, show=True, mask_land=[]);

K-MEANS Clustering#

Daily synoptic patterns are obtained using the K-means clustering algorithm. It divides the data space into 81 clusters, a number that which must be a compromise between an easy handle characterization of the synoptic patterns and the best reproduction of the variability in the data space. Previous works with similar anlaysis confirmed that the selection of this number is adequate (Rueda et al. 2017).

Each cluster is defined by a prototype and formed by the data for which the prototype is the nearest.

Finally the best match unit (bmus) of daily clusters are reordered into a lattice following a geometric criteria, so that similar clusters are placed next to each other for a more intuitive visualization.

# KMeans clustering for slp data
kma, kma_order, bmus, sorted_bmus, centers = kmeans_clustering_pcs(
    PCs, n_pcs, 
    num_clusters, 
    min_data,
    kma_tol = 0.0001, 
    kma_n_init = 10,
)

# store kma output
pk.dump([kma, kma_order, bmus, sorted_bmus, centers], open(p_kma_slp_model, 'wb'))

# load kma 
#kma, kma_order, bmus, sorted_bmus, centers = pk.load(open(p_kma_slp_model, 'rb'))

Iteration: 422
Number of sims: 152
Minimum number of data: 5

# add KMA data to SLP
slp['bmus'] = ('time', sorted_bmus)
slp['kma_order'] = kma_order
slp['cluster'] = np.arange(num_clusters)
slp['centers'] = (('cluster', 'dims'), centers)

# add estela mask to SLP
slp['est_mask'] = (('latitude','longitude'), estela_mask.values)

# store updated dataset
slp.to_netcdf(p_kma_slp_classification)

Mean#

The mean SLP associated to each of the different clusters is shown below.

Plot_slp_kmeans(slp, extent=(50, 290, -67, 65), figsize=[23,14]);

Anomaly#

The SLP anomaly with respect to the mean and associated to each of the different clusters is shown below.

Plot_slp_kmeans(slp, extent=(50, 290, -67, 65), figsize=[23,14], anomaly=True);

Seasonality#

Probability of each KMA cluster, conditioned to each month and season

Plot_DWTs_Probs(slp.bmus+1, pca.pred_time.values, num_clusters, vmax=0.1);

Tailor-Made classification of large scale predictor

Contents