HyBeat metamodel: selection of cases, generating a library of pre-run cases and PCA analysis

HyBeat is a metamodel based on data mining techniques and a reduced number of 2D hydrodynaimc model Xbeach (Roelvink et al., 2015) numerical simulations, that allows to obtain the mean setup and the infragravity wave at a very high resolution (~ 5-20m) at the coast.

These are the steps followed in the HyBeat methodology:

• Selection of 100 design cases based on selected parameters (Hs,Tp, θ, 𝜂) using Latin Hybercube Sampling (LHS) algorithm and Maximum Dissimilarity Algorithm (MDA) selection technique
  

• Numerically run the 100 design cases per each computational grid with the 2D model Xbeach at a spatial resolution of 5 to 20 meters
  

• Post process the 100 simulations per computational domain by means of Principal Component Analysis (PCA)
  

• Spatial and temporal reconstruction of mean wave setup and infragravity wave using a non-linear interpolation technique based on Radial Basis Functions (RBFs)
  

HyBeat

In this notebook we are going to focus on showing significant intermediate results of HyBeat metamodel to be able to reconstruct any given event

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

# pip
from scipy.stats import qmc
import pandas as pd
import numpy as np
from sklearn.cluster import KMeans


# bluemath 
sys.path.insert(0, os.path.join(os.path.abspath(''),  '..'))
from Codes_v3.Codes.LHS import *
from Codes.mda import MaxDiss_Simplified_NoThreshold
from Codes_v3.Codes.plots import *
sys.path.append('/home/administrador/PROYECTOS/SamoaTonga/Xbeach/teslakit-master')
from teslakit.rbf import RBF_Reconstruction

LHS + MDA

A selection of cases based on selected parameters (Hs,Tp, θ, 𝜂) and its ranges is done through Latin Hypercube Sampling algorithm (LHS, McKay et al., 1979, Sheikholeslami and Razavi, 2017) and a selection algorithm, Maximum Dissimilarity Algorithm (MDA, Camus et al., 2011a). LHS cuts each dimension space into n (10,000) sections, where n is the number of sampling points, and only one combination in each section is kept. It produces samples that reflect the true underlying distribution, and it tends to require much smaller sample sizes than simple random sampling. Then Hs,Tp, θ and 𝜂 variables are chosen to apply MDA, with the aim of selecting a representative subset of 100 cases which best captures the variability of the space.

HyBeat

The following intermediate results are for Nukualofa domain in Tongatapu island

#Ocean parameters
n_dims = 4 #Number of dimensions
name_dims=['Hs', 'Steepness', 'Dir', 'SWL']
lower_bounds = [1, 0.005, -45+30, -0.7] 
upper_bounds = [10, 0.05, 45+30, 2.5]
n_samples=10000 #Number of combinations to extract from LHS

#LHS params and selection
sampler = qmc.LatinHypercube(d=n_dims,seed=0)
dataset = sampler.random(n=n_samples)
dataset = qmc.scale(dataset, lower_bounds, upper_bounds)
dataset_pd = pd.DataFrame(dataset, columns=name_dims).to_xarray()
dataset_pd['Tp'] = (('index'), np.sqrt((dataset_pd.Hs.values*2*np.pi)/(9.81*dataset_pd.Steepness.values)))

#MDA selection
n_subset = 100 #Subset for MDA
ix_scalar = [0,1,3]
ix_directional = [2]

subset_mda = MaxDiss_Simplified_NoThreshold(dataset, n_subset, ix_scalar, ix_directional)
df = pd.DataFrame(subset_mda, columns=name_dims).to_xarray()
df['Tp'] = (('index'), np.sqrt((df.Hs.values*2*np.pi)/(9.81*df.Steepness.values)))

MaxDiss dataset: 10000 --> 100

   MDA centroids: 100/100

n=100
Plot_MDA_Data(pd.DataFrame(dataset), pd.DataFrame(subset_mda[:n,:]), d_lab=name_dims, figsize=[15,12])

dataset_plot = np.vstack([dataset_pd['Hs'].values, dataset_pd['Tp'].values, dataset_pd['Dir'].values,dataset_pd['SWL'].values]).T
mda_plot = np.vstack([df['Hs'].values, df['Tp'].values, df['Dir'].values,df['SWL'].values]).T
name_plot=['Hs', 'Tp', 'Dir', 'SWL']

Plot_MDA_Data(pd.DataFrame(dataset_plot), pd.DataFrame(mda_plot[:n,:]), d_lab=name_plot, figsize=[15,12])

Library of Xbeach pre-run cases

A set of computational grids are carefully design according to the HR topo-bathymetry data (5 m) so that the wave setup results along the entire coastline of the islands is well defined without instabilities due to lateral boundary conditions. The following figure shows the domain of each computational grid for Upolu island, as well as the overlapping areas between them. These computational grids are generated using Delft3D (https://oss.deltares.nl/web/delft3d) meshing tools RFGRID and QUICKIN, which enables the generation of rectilinear and curvilinear grids from higher (~ 80m) to lower resolution (~ 5-20m) near the coast. Finally, each model run simulates in surfbeat mode a one-hour event according the 100 design cases, which Hydraulic Boundary Conditions (HBCs) are defined by a parametric Jownsap spectra and a constant water level (Hs,Tp, θ, 𝜂).

HyBeat

The following intermediate results are for Nukualofa domain in Tongatpu island

p_results = '/media/administrador/HD/SamoaTonga/xbeach/Cluster/Results/Postprocess_results/Tongatapu/PCA_results/TX02'
subset_mda = xr.open_dataset(os.path.join(p_results, 'LHS_mda_sel_XBeach_TX02.nc'))
C_R_zs_mean = xr.open_dataset(os.path.join(p_results, 'PCA_zs_mean_ncases_100.nc'))
C_R_zs_var = xr.open_dataset(os.path.join(p_results, 'PCA_zs_var_ncases_100.nc'))
C_R_H_mean = xr.open_dataset(os.path.join(p_results, 'PCA_H_mean_ncases_100.nc'))

fig = plt.figure(figsize=[10,10])
gs = gridspec.GridSpec(1,1, hspace = 0.005, wspace = 0.005)
ax = fig.add_subplot(gs[0,0])
ax.set_title('Topo-bathymetry', fontsize=20, color='navy')
var_name = 'zb'
Plot_Maps_XBeach(C_R_zs_mean, var_name,p_coast =1, ax=ax, mask=0)

ncases = len(C_R_zs_mean.case)
for i in range(ncases):
    
    
    textstr = '\n'.join((
    r'$\mathrm{Hs}=%.2f \mathrm{\; m}$' % (round(subset_mda.Hs.values[i],2), ),
    r'$\mathrm{Tp}=%.2f \mathrm{\; s}$' % (round(subset_mda.Tp.values[i],2), ),
    r'$\mathrm{Dir}=%.2f$' % (round(subset_mda.Dir.values[i],2), ),
    r'$\eta=%.2f \mathrm{\; m}$' % (round(subset_mda.SWL.values[i],2), )))
    
    if i==0 or i==5 or i==10 or i==21 or i==30 or i==34 or i==44 or i==49 or i==67 or i==90  or i==27:
        continue
    
    props = dict(boxstyle='round', facecolor='white', alpha=0.5)
    # props = dict(facecolor='white', alpha=0.5)
    
    
    fig = plt.figure(figsize=[30,5])
    gs = gridspec.GridSpec(1,3, hspace = 0.005, wspace = 0.005)
    
    ax = fig.add_subplot(gs[0])
    ax.set_title('Hs mean. Case '+str(i), fontsize=15, color='navy')
    Plot_Maps_XBeach(C_R_H_mean.isel(case=i)*1.4, vmin=0, vmax=10, var_name = 'H_mean',p_coast =1, ax=ax)
    ax.text(0.79, 1, textstr, transform=ax.transAxes, fontsize=14,verticalalignment='top', bbox=props)
    
    
    ax = fig.add_subplot(gs[1])
    ax.set_title('Zs mean. Case '+str(i), fontsize=15, color='navy')
    Plot_Maps_XBeach(C_R_zs_mean.isel(case=i), vmin=-1.5, vmax=1.5, var_name = 'zs_mean',p_coast =1, ax=ax)
    # ax.text(0.82, 0.95, textstr, transform=ax.transAxes, fontsize=14,verticalalignment='top', bbox=props)

    
    ax = fig.add_subplot(gs[2])
    ax.set_title('Zs std. Case '+str(i), fontsize=15, color='navy')
    Plot_Maps_XBeach(C_R_zs_var.isel(case=i), vmin=0, vmax=0.5, var_name = 'zs_var',p_coast =1, ax=ax)
    # ax.text(0.82, 0.95, textstr, transform=ax.transAxes, fontsize=14,verticalalignment='top', bbox=props)

Principal Component Analysis (PCA)

The Principal Component Analysis (PCA) is applied to the hourly fields of mean wave setup and its standard deviation to infer the infragravity wave component. This statistical technique is widely used in climatology to identify dominant variability patterns and reduce dimensionality (Camus et al. 2014). Hence, by PCA the main variability modes (patterns) of the target variables and the temporal coefficients of the identified patterns. PCA projects the original data on a new space, searching for the maximum variance of the sample data. The eigenvectors (empirical orthogonal functions, EOFs) of the data covariance matrix define the vectors of the new space. The transformed components of the original data over the new vectors are the principal components (PCs). Then, the target variables can be expressed as a linear combination of EOFs and PCs. As an example, following figure shows the first four empirical EOFs and PCs of the mean wave setup variable for one of the computational domains of Tongatapu, Nukualofa.

HyBeat

The following intermediate results are for Nukualofa domain in Tongatapu island

PCA of mean wave setup

Plot_EOFsMEAN(C_R_zs_mean, vmin=0, vmax=1,mask=1)

for n_pc in range(C_R_zs_mean.n_pcs.values):
    Plot_EOFs(C_R_zs_mean, n_pc=n_pc, vmin=-0.01, vmax=0.01,mask=1)

PCA of standard deviation of wave setup

Plot_EOFsMEAN(C_R_zs_var, vmin=0, vmax=0.05,mask=1)

for n_pc in range(C_R_zs_mean.n_pcs.values):
    Plot_EOFs(C_R_zs_var, n_pc=n_pc, vmin=-0.01, vmax=0.01,mask=1)

Reconstruction validation based pm RBFs

Finally, the spatial and temporal reconstruction of these variables for any given sea state are obtained using linear combination of EOFs and PCs combined with a non-linear interpolation technique based on radial basis functions (RBFs). This interpolation technique is very convenient for scattered and multivariate data (see details in Camus et al., 2011b). The numerical validation of the results (see next figrues) confirms the ability of the developed methodology to spatially reconstruct wave setup in coral reef areas at nearshore locations with a considerable reduction in the computational effort.

HyBeat

The following validation results are for mean wave setup variable in Nukualofa domain (Tongataopu island)

dataset = xr.open_dataset(os.path.join(p_results, 'LHS_mda_sel_XBeach_TX02.nc'))

var_rec = ['zs_mean']
vars_mda = ['Hs','Steepness','Dir','SWL']
subset=subset_mda.to_dataframe().loc[:,vars_mda].values
ix_scalar_subset = [0,1,3]
ix_directional_subset = [2]

target = pd.DataFrame(C_R_zs_mean.PCs[:,:C_R_zs_mean.n_pcs.values].values).values
ix_scalar_target = list(range(C_R_zs_mean.n_pcs.values))
ix_directional_target = []

dataset1=dataset.to_dataframe().loc[:,vars_mda].values

output = RBF_Reconstruction(subset, ix_scalar_subset, ix_directional_subset, target, ix_scalar_target, ix_directional_target, dataset1)

ix_scalar: 0,  optimization: 0.31 | interpolation: 0.00
ix_scalar: 1,  optimization: 0.21 | interpolation: 0.00
ix_scalar: 2,  optimization: 0.24 | interpolation: 0.00
ix_scalar: 3,  optimization: 0.24 | interpolation: 0.00

Z_ALL = RBF_var_reconstruct(C_R_zs_mean, output, 0, len(dataset.index), var=var_rec[0])
dataset=dataset.rename({'SWL':'zs0'})
dataset[var_rec[0]]=(('index', 'ny', 'nx'), Z_ALL[var_rec[0]].values)
dataset['globalx']=(('ny','nx'), C_R_zs_mean.globalx.values)
dataset['globaly']=(('ny','nx'), C_R_zs_mean.globaly.values)
dataset['zb']=(('ny','nx'), C_R_zs_mean.zb.values)
mask=np.where(dataset.zb<0,1,np.nan)
dataset['mask']=(('ny','nx'), mask)

for i in range(len(dataset.index)):
    
    fig = plt.figure(figsize=[30,5])
    gs = gridspec.GridSpec(1,3, hspace = 0.005, wspace = 0.005)
    if i==0 or i==5 or i==10 or i==21 or i==30 or i==34 or i==44 or i==49 or i==67 or i==90  or i==27:
        continue
    ax = fig.add_subplot(gs[0])
    ax.set_title('XBeach Simulation', fontsize=20, color='navy')
    
    
    Plot_Maps_XBeach(C_R_zs_mean.isel(case=i), vmin=-1.2, vmax=1.2, var_name = 'zs_mean',p_coast =1, ax=ax)

    cd = np.where((np.round(subset[i,0],5)==np.round(dataset.Hs.values,5)) & (np.round(subset[i,1],5)==np.round(dataset.Steepness.values,5)) 
         & (np.round(subset[i,2],5)==np.round(dataset.Dir.values,5)) & (np.round(subset[i,3],5)==np.round(dataset.zs0.values,5)))[0]
    
    ax = fig.add_subplot(gs[1])
    ax.set_title('RBF Reconstruction', fontsize=20, color='navy')
    Plot_Maps_XBeach(dataset.isel(index=cd).squeeze(),vmin=-1.2, vmax=1.2, var_name = 'zs_mean',p_coast =1, ax=ax, mask=1)
    
    
    ax = fig.add_subplot(gs[2])
    ax.set_title('XBeach Simulation - XBeach Reconstruction', fontsize=20, color='navy')
    
    Plot_Maps_XBeach(C_R_zs_mean.isel(case=i)-dataset.isel(index=cd).squeeze(),vmin=-0.5, vmax=0.5, var_name = 'zs_mean',p_coast =1, ax=ax)

    

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>

<Figure size 2160x360 with 0 Axes>