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 = 'Samoa'
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 = [-172.92+360, -171.3+360, -14.14, -13.3]
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 = [12, 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,
    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 = 3.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#

def Plot_hindcast_elnino(data, quant=None, var_s = 'hs', vmin=0, vmax=2.5,
                         quiv_color = 'teal', figsize = [26, 20]):
    '''
    Plot hindcast El Nino events. Needs "n34" index  variable inside data with values
    "1,2,3" for "ElNino, Neutral, LaNina"

    data         - xarray.Dataset with "n34" index values
    p_coastfile  - path to .shp file from GSHHS
    quant        - quantile, None for mean, float to calculate quantile
    var_s        - variable to plot
    vmin         - minimun value for colormap
    vmax         - maximun value for colormap
    '''

    # generate figure and gridsizer
    fig = plt.figure(figsize=figsize)
    gs = gridspec.GridSpec(1, 3, hspace=0.01, wspace=0.01)

    # years type
    # TODO enhe en el label
    years = ['El Niño', 'Neutral', 'La Niña']
    for s in range(len(years)):

        # add subplot and select year type
        ax = fig.add_subplot(
            gs[s],
            projection = ccrs.PlateCarree(central_longitude=180),
        )

        # mean or quantile
        if quant == None:
            s_p = data.isel(time = np.where(data.n34.values==s)[0]).mean(dim = 'time')

        else:
            s_p = data.isel(time = np.where(data.n34.values==s)[0]).quantile(0.99, dim = 'time')

        # plot season hindcast
        Plot_hindcast_maps(
            s_p, var_s,
            vmin = vmin, vmax = vmax,
            quiv_fact = 10, quiv_color = quiv_color,
            ax = ax,
        )

        ax.set_title(years[s], fontsize=16, color='firebrick', fontweight='bold')

    return fig

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

_images/05_Hindcast_Reconstruction_35_0.png

Extreme waves: 99% percentile#

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

_images/05_Hindcast_Reconstruction_37_0.png