KMA Reconstruction#

# common
import warnings
warnings.filterwarnings('ignore')
import os
import os.path as op
import sys
import pickle as pkl

# pip
import xarray as xr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

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

# bluemath modules 
from bluemath.binwaves.reconstruction import generate_kp_coefficients, reconstruct_kmeans_clusters_hindcast

from bluemath.binwaves.plotting.swan_case import Plot_swan_case
from bluemath.binwaves.plotting.buoy import Plot_comparison_buoy, Plot_buoy_location
from bluemath.binwaves.plotting.spectra import Plot_example_spectra, Plot_spectra_propagation

Warning: cannot import cf-json, install "metocean" dependencies for full functionality

Database and site parameters#

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

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

p_kma = op.join(p_deliv, 'kma')
p_swan = op.join(p_deliv, 'swan')

# SWAN simulation parameters
name = 'binwaves'
p_swan_subset = op.join(p_swan, name, 'subset.nc')
p_swan_output = op.join(p_swan, name, 'out_main_binwaves.nc')
p_swan_input_spec = op.join(p_swan, name, 'input_spec.nc')

# SuperPoint KMeans
num_clusters = 2000
p_superpoint_kma = op.join(p_kma, 'Spec_KMA_{0}_NS.nc'.format(num_clusters))

# buoy data
point_buoy = [187.8, -13.88]
p_buoy = op.join(p_site, 'buoy', 'samoaapolima_lon{0}_lat{1}_year1989_1990.pkl'.format(point_buoy[0], point_buoy[1]))
                 

# output files
p_site_out = op.join(p_deliv, 'reconstruction')

p_kp_coeffs = op.join(p_site_out, 'kp_grid')          # SWAN propagation coefficients
p_efth_kma = op.join(p_site_out, 'efth_kma')          # KMA reconstruction
p_reco_hnd = op.join(p_site_out, 'hindcast')          # hindcast reconstruction folder

# aux functions

def load_point_kp(point, out_sim, p_kp_coeffs):
    '''
    Load kp propagation coefficients for a choosen point
    
    point   - point to load [lon, lat] (find nearest)
    out_sim - SWAN simulation output
    '''
    
    # find point 
    ilon = np.argmin(np.abs(out_sim.lon.values - point[0]))
    ilat = np.argmin(np.abs(out_sim.lat.values - point[1]))

    # load coefficient
    fn = 'kp_lon_{0}_lat_{1}.nc'.format(out_sim.lon.values[ilon], out_sim.lat.values[ilat])
    p_cf = op.join(p_kp_coeffs, fn)
    
    return xr.open_dataset(p_cf)

def load_point_efth_kma(point, out_sim, p_efth):
    '''
    Load spectra efth_kma for a choosen point
    
    point - point to load [lon, lat] (find nearest)
    out_sim - SWAN simulation output
    '''
    
    # find point 
    ilon = np.argmin(np.abs(out_sim.lon.values - point[0]))
    ilat = np.argmin(np.abs(out_sim.lat.values - point[1]))

    # load coefficient
    fn = 'efth_lon_{0}_lat_{1}.nc'.format(out_sim.lon.values[ilon], out_sim.lat.values[ilat])
    p_cf = op.join(p_efth, fn)
    
    return xr.open_dataset(p_cf)

def load_hindcast_recon(point, out_sim, p_reco):
    '''
    Load hindcast reconstruction for a choosen point
    
    point - point to load [lon, lat] (find nearest)
    out_sim - SWAN simulation output
    '''
    
    # find point 
    ilon = np.argmin(np.abs(out_sim.lon.values - point[0]))
    ilat = np.argmin(np.abs(out_sim.lat.values - point[1]))

    # load coefficient
    fn = 'reconst_hind_lon_{0}_lat_{1}.nc'.format(out_sim.lon.values[ilon], out_sim.lat.values[ilat])
    p_cf = op.join(p_reco, fn)
    
    return xr.open_dataset(p_cf)

KMA clusters#

We reconstruct the spec in the grid points, for the 2000 spectral centroids of the cluster

sp_kma = xr.open_dataset(p_superpoint_kma)

One cluster example OFFSHORE

# example cluster
cluster = 4

# Plot spectra
Plot_example_spectra(sp_kma, cluster=cluster);

_images/04_KMeans_Reconstruction_14_0.png

Simulations preprocess#

This is where we postprocess all the simulations, to obtain the propagation coefficients (kp) for all the 696 cases and for each point of the grid
This process takes a long time but it only needs to be run once

# generate propagation coefficients for all points and cases
generate_kp_coefficients(p_kp_coeffs, out_sim, sp_kma, override=False)

Preprocessing kp coefficients...: 100%|██████████████████████████████| 36846/368

Explanation

Here are a couple examples of the process made here for one of the propagation cases.
In the plots below, the propagation for a case out of the 696 has been plotted with the propagations for two different points of the grid (The ones marked with a star)

The first example is for the buoy location, which receives all the energy from the 142.5º direction. The Kp coefficient is 1 with the same direction as the input offshore spectra, meaning that any energy is lost.
The second example is for a location in the northern part of Tognatapu, where waves are shadowed by the many small islands in the area. For this reason the same input spectra is here transformed into 4 different systems with different directions, kist a kp coefficient of around 0.1/0.2

# Select one propagation case
case = 250
point_example = [187.93, -13.76]

# Load spectra kp
output_spec = load_point_kp(point_example, out_sim, p_kp_coeffs)

# Plot SWAN case 
Plot_swan_case(out_sim, subset, case, point_example[0], point_example[1]);

# Plot SWAN Hs propagation spectra
Plot_spectra_propagation(input_spec, output_spec, sp_kma, case, ylim=0.5, cmap='gnuplot2_r');

_images/04_KMeans_Reconstruction_19_0.png

_images/04_KMeans_Reconstruction_19_1.png

# Select one propagation case
case = 250
point_example = [187.93, -14.00]

# Load spectra kp
output_spec = load_point_kp(point_example, out_sim, p_kp_coeffs)

# Plot SWAN case 
Plot_swan_case(out_sim, subset, case, point_example[0], point_example[1]);

# Plot SWAN Hs propagation spectra
Plot_spectra_propagation(input_spec, output_spec, sp_kma, case, ylim=0.5, cmap='gnuplot2_r');

_images/04_KMeans_Reconstruction_20_0.png

_images/04_KMeans_Reconstruction_20_1.png

Reconstruction based on KMA#

Validation#

Below is the Buoy Location available for validation at our study site

# Plot buoy location
Plot_buoy_location(point_buoy, area_extraction);

_images/04_KMeans_Reconstruction_24_0.png

Buoy Location 1#

# reconstruct model data for that specific point
reconstruct_kmeans_clusters_hindcast(
    p_site_out, out_sim, sp_kma, 
    point = point_buoy, 
    reconst_hindcast = True,
    override = False,
)

Preprocessing efth and stats...: 100%|██████████████████████████████| 1/1 [00:23

Same cluster as in the example above, propagated to buoy at Location 1

# Load propagated spectra
efth_point = load_point_efth_kma(point_buoy, out_sim, p_efth_kma)

# Plot spectra
Plot_example_spectra(efth_point, cluster=cluster);

_images/04_KMeans_Reconstruction_28_0.png

# load buoy data and parse it to xarray.Dataset
with open(p_buoy, 'rb') as fr:
    buoy = pkl.load(fr)

buoy_data = buoy.to_xarray().isel(index = np.where((buoy.hs > 0) & (buoy.tm10 > 0))[0])

# Load hindcast reconstruction at nearest point to buoy
recon_data = load_hindcast_recon(point_buoy, out_sim, p_reco_hnd)

# compare hindcast reconstruction with buoy data
Plot_comparison_buoy(recon_data, buoy_data, ss=3);

_images/04_KMeans_Reconstruction_31_0.png

KMA Reconstruction

Contents