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

Obtain wave metrics#

if wave_metrics:  
    
    # calculate wave metrics
    hindcast_stats, hindcast_monthly = obtain_wave_metrics_grid(
        p_reco_hnd,
        out_sim, 
        perc = wave_metrics_perc,
    )
    
    # store metrics and monthly mean
    hindcast_stats.to_netcdf(p_hindcast_stats)
    hindcast_monthly.to_netcdf(p_hindcast_monthly)

else:
    
    # load metrics and monthly mean
    hindcast_stats = xr.open_dataset(p_hindcast_stats)
    hindcast_monthly = xr.open_dataset(p_hindcast_monthly)
Preprocessing efth and stats...: 100%|██████████████████████████████| 36846/36846 [19:36:46<00:00,  1.92s/it]                                                                                         
hindcast_monthly
<xarray.Dataset>
Dimensions:  (lat: 138, lon: 267, time: 502)
Coordinates:
  * lat      (lat) float64 -14.14 -14.13 -14.13 -14.12 ... -13.33 -13.32 -13.32
  * lon      (lon) float64 187.1 187.1 187.1 187.1 ... 188.7 188.7 188.7 188.7
  * time     (time) datetime64[ns] 1979-01-31 1979-02-28 ... 2020-10-31
Data variables:
    hs       (lon, lat, time) float64 1.623 1.677 1.651 ... 2.11 1.912 1.481
    tp       (lon, lat, time) float64 9.679 13.95 14.32 ... 9.45 8.775 8.461
    tm       (lon, lat, time) float64 6.437 7.345 8.069 ... 6.3 6.054 5.752
    dm       (lon, lat, time) float64 159.4 184.9 191.2 ... 108.6 108.9 87.7
    dp       (lon, lat, time) float64 205.1 263.8 248.9 ... 123.4 139.2 82.5
    dspr     (lon, lat, time) float64 68.23 61.66 63.1 ... 41.84 44.72 38.43
# Plot hindcast at selected point
point = [187.8, -13.88]

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

# monthly hindcast at point
ilon = np.argmin(np.abs(hindcast_monthly.lon.values - point[0]))
ilat = np.argmin(np.abs(hindcast_monthly.lat.values - point[1]))
month_data = hindcast_monthly.isel(lon=ilon, lat=ilat)

# plot hs, tp and dm time series
Plot_hindcast_series_point(recon_data, month_data, 'hs');
Plot_hindcast_series_point(recon_data, month_data, 'tp');
Plot_hindcast_series_point(recon_data, month_data, 'dm');
_images/05_Hindcast_Reconstruction_16_0.png _images/05_Hindcast_Reconstruction_16_1.png _images/05_Hindcast_Reconstruction_16_2.png

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