Contents

Waves Hindcast

Contents

Waves Hindcast#

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

# 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 reconstruct_kmeans_clusters_hindcast, obtain_wave_metrics_grid
from bluemath.binwaves.reconstruction import get_hindcast_snapshot
from bluemath.binwaves.nino34 import n34_format

from bluemath.binwaves.plotting.hindcast import Plot_hindcast_series_point, Plot_hindcast_snapshot, Plot_hindcast_maps
from bluemath.binwaves.plotting.hindcast import Plot_hindcast_seasonality, Plot_hindcast_elnino
from bluemath.binwaves.plotting.nino34 import Plot_n34

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

Methodology#

Hindcast

Once we have reconstructed the K-Means clusters for the KMA clusters, and validated the results against instrumental data, we can efficiently reconstruct the hindcast following the steps below

TITLE

Database and site parameters#

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

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

# nino34 dataset
p_n34 = op.join(p_data, 'ninho34.txt')

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


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

p_reco_hnd = op.join(p_site_out, 'hindcast')                     # hindcast reconstruction folder
p_hindcast_stats = op.join(p_reco_hnd, 'hindcast_stats.nc')      # hindcast reconstruction metrics
p_hindcast_monthly = op.join(p_reco_hnd, 'hindcast_monthly.nc')  # hindcast monthly mean

# hindcast reconstruction
recon_hindcast = False
wave_metrics = True
wave_metrics_perc = [25, 50, 75, 90, 95, 99]

# hindcast snapshot
hindcast_snapshot = True
time_snapshot = '1993-05-12-00:00:00'

# aux functions

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)

# Load SWAN simulation output and kmeans clusters
out_sim_swan = xr.open_dataset(p_swan_output)
sp_kma = xr.open_dataset(p_superpoint_kma)

Here we need to define the area for the extraction and the resample factor
The plot is an example of the area that has been cut

# this will reduce memory usage
area_extraction=[184.5, 185.3, -21.62, -20.97]
resample_factor = 2

# cut area
out_sim = out_sim_swan.sel(
    lon = slice(area_extraction[0], area_extraction[1]), 
    lat = slice(area_extraction[2], area_extraction[3]),
)

# resample extraction points
out_sim = out_sim.isel(
    lon = np.arange(0, len(out_sim.lon), resample_factor), 
    lat = np.arange(0, len(out_sim.lat), resample_factor),
)

# check that the area is OK
plt.figure(figsize = [10, 5])
out_sim.isel(case = 10).Hsig.plot(cmap='RdBu', vmin=0, vmax=2);

_images/05_Hindcast_Reconstruction_10_0.png

Reconstruct Hindcast#

# reconstruct model data for that area
if recon_hindcast:
    
    reconstruct_kmeans_clusters_hindcast(
        p_site_out, out_sim, sp_kma, 
        reconst_hindcast = True,
        override = False,
    )
    

Wave climate#

Snapshot analysis#

if hindcast_snapshot:
    data_snapshot = get_hindcast_snapshot(p_reco_hnd, out_sim, time_snapshot)

    # plot hindcast snapshot
    Plot_hindcast_snapshot(data_snapshot) 

_images/05_Hindcast_Reconstruction_19_0.png

Average climate#

Plot_hindcast_maps(
    hindcast_stats, 'Dir_m', 
    vmin = 0, vmax = 360,
    quiv_color = 'palevioletred', 
    figsize = [15, 15],
);

_images/05_Hindcast_Reconstruction_21_0.png

Plot_hindcast_maps(
    hindcast_stats, 'Hs', 
    percentile = 50,
    vmin = 0, vmax = 2.5,
    quiv_fact = 8, figsize = [15, 15],
);

_images/05_Hindcast_Reconstruction_22_0.png

Plot_hindcast_maps(
    hindcast_stats, 'Hs', 
    percentile = 90,
    vmin = 0, vmax = 3.5,
    quiv_fact = 8, figsize = [15, 15],
);

_images/05_Hindcast_Reconstruction_23_0.png

Plot_hindcast_maps(
    hindcast_stats, 'Hs', 
    percentile = 99,
    vmin = 0, vmax = 5,
    quiv_fact = 8, quiv_color = 'teal',
    figsize = [15, 15],
);

_images/05_Hindcast_Reconstruction_24_0.png

Plot_hindcast_maps(
    hindcast_stats, 'Tm', 
    percentile = 50,
    vmin = 5, vmax = 10,
    quiv_fact = 8, quiv_color = 'cornsilk',
    figsize = [15, 15],
);

_images/05_Hindcast_Reconstruction_25_0.png

Seasonality#

seasons = hindcast_monthly.groupby('time.season').mean()

Plot_hindcast_seasonality(seasons, var_s = 'hs', vmin = 0, vmax = 2.7);
    

_images/05_Hindcast_Reconstruction_28_0.png

Plot_hindcast_seasonality(seasons, var_s = 'tm', vmin = 4, vmax = 10, quiv_color = 'cornsilk');
    

_images/05_Hindcast_Reconstruction_29_0.png

Plot_hindcast_seasonality(seasons, var_s = 'tp', vmin = 6.5, vmax = 16, quiv_color = 'cornsilk');

_images/05_Hindcast_Reconstruction_30_0.png

Interannual variability#

# load ninho34 data
n34 = np.loadtxt(p_n34,skiprows=1, max_rows=74)
n34 = n34_format(n34, rolling_mean=True)

n34 = n34.sel(time = hindcast_monthly.time)

# plot ninho34 and thresholds
thr1 = 1   # over this value, years are classified as El Nino (1)
thr2 = -1  # under this value, years are classified as La Nina (3)
Plot_n34(n34, l1 = thr1, l2 = thr2, figsize=[19, 6]);

# add ninho34 classification to monthly hindcast data
hindcast_monthly['n34'] = (('time'), n34.classification.values) # 1:El Nino, 2:Neutral, 3:La Nina

_images/05_Hindcast_Reconstruction_32_0.png

Average waves 50% percentile#

Plot_hindcast_elnino(hindcast_monthly, var_s = 'hs', vmin = 0, vmax = 2.5, figsize=[25,12]);

_images/05_Hindcast_Reconstruction_34_0.png

Extreme waves: 99% percentile#

Plot_hindcast_elnino(hindcast_monthly, quant = 0.99, var_s = 'hs', vmin = 0, vmax = 3, figsize=[25,12]);

_images/05_Hindcast_Reconstruction_36_0.png