Historical TC : Evan (2012)#

# common
import os
import os.path as op
import sys

# pip
import numpy as np
import xarray as xr
import wavespectra

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

# bluemath modules
from bluemath.greensurgerepo.storms import historic_track_preprocessing, historic_track_interpolation, get_category
from bluemath.greensurgerepo.vortex import vortex_model_general
from bluemath.greensurgerepo.greensurge import GS_wind_partition, print_GFD_properties, GS_windsetup_reconstruction, presure_to_IB

from bluemath.binwaves.forecast import cut_resample_area, reconstruct_spec

from bluemath.greensurgerepo.plots.plots_greensurge import plot_case_input, plot_presure_vs_IB

from bluemath.astronomical_tide import calculate_AT_TPXO9v4

from bluemath.wswan.plots.common import custom_cmap, bathy_cmap

from bluemath.tc_forecast import Plot_Swath_ShyTCWaves, Plot_Swath_ShyTCWaves_coords, \
Plot_SuperPoint_TC_Forecast, Plot_Grid_HsTmDp, Plot_Grid_TWL_max, Plot_Grid_SWATH, \
Plot_TWL_timeseries, Plot_Grid_HsPd_series
from bluemath.tc_forecast import animate_tc_hs_ss_tide_twl_forecast

# hide warnings
import warnings
warnings.filterwarnings('ignore')

## Methodology

BinWaves

In this section we will make use of the combination of SHyTCWaves and Greensurge at a regional scale, and BinWaves at a local scale, to obtain the evolution of waves and water levels caused by TC Evan in 2012

TITLE

Database and site parameters#

p_data = r'/media/administrador/HD2/SamoaTonga/data'
site = 'Samoa'
p_site = op.join(p_data)

#path databases
p_db = r'/media/administrador/HD1/DATABASES'

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

# historical TCs (ibtracs)
p_ibtracs = op.join(p_deliv, 'IBTrACS.ALL.v04r00.nc')

# tc forecast
p_tcforecast = op.join(p_deliv, 'TC_forecast')

# ShyTCWaves
p_shy_grid = os.path.join(p_tcforecast, 'Super_Points', 'xds_shy_bulk_Samoa_EVAN_2012.nc')
p_superpoint_shytcwaves = os.path.join(p_tcforecast, 'Super_Points', 'xds_shy_sp_Samoa_EVAN_2012.nc')

# BinWaves: superpoint and reconstruction coefficients
deg_sup = 15
area_extraction =  [-172.85+360, -171.38+360, -14.09, -13.41]
resample_factor = 10

p_deliv_d05 = op.join(p_site, 'd05_swell_inundation_forecast')
p_superpoint_binwaves = op.join(p_deliv_d05, 'spec', 'super_point_superposition_{0}_.nc'.format(deg_sup))
p_swan_output = op.join(p_deliv_d05, 'swan', 'binwaves', 'out_main_binwaves.nc')
p_kp_coeffs = op.join(p_deliv_d05, 'reconstruction', 'kp_grid')

# GreenSurge
domain_study='Samoa'
p_greensurge = op.join(p_db, 'greensurge')
p_GFD_libdir = op.join(p_greensurge, 'GFD_lib', site, 'GFD_Samoa_T84_W40_D15_24h_CDdefault')
p_GFD_info = op.join(p_GFD_libdir, 'GFD_lib_info.nc')

# astronomical tide
model_directory = op.join(p_db, 'TPX09_atlas_v4')


# model execution parameters
run_binwaves = False
run_green_WS = True

TC Selection#

With the name and the year we can select here any historical TC from IBTrACS

ibtracs = xr.open_dataset(p_ibtracs)

name = b'EVAN'
year = 2012
centerID = 'WMO'

storm = ibtracs.isel(
    storm = np.where(
        (ibtracs.name == name) & \
        (ibtracs.time[:,0].dt.year.values == year)
    )[0]
).isel(storm=0)

# preprocess selected historic TC: obtain strom variables
df_TC_hist = historic_track_preprocessing(storm, center = centerID)

# computational time step mandatory for GS methodology: 
dt_comp = 20  # minutes

# generate interpolated storm track and mandatory data to apply vortex model 
storm_track, time_input = historic_track_interpolation(df_TC_hist, dt_comp)

plot_case_input(storm_track, name);

ShyTCWaves#

SHyTCWaves

Using SHyTCWaves we are able of obtaining the SWATH for the wave parameters and the evolution of the SuperPoint directional spectra over time.

SWATH#

TC_grid = xr.open_dataset(p_shy_grid)

Plot_Swath_ShyTCWaves(
    TC_grid, storm_track, 
    name, 
    var = 'hswath', 
    cmap = custom_cmap(50, 'YlOrRd', 0.15, 0.9, 'YlGnBu_r', 0, 0.85),
);

SuperPoint#

Locations od the SuperPoint surrounding our study site

sp = xr.open_dataset(p_superpoint_shytcwaves)

Plot_Swath_ShyTCWaves_coords(
    TC_grid, storm_track, 
    name, 
    sp,
    mksize = 20,
    var = 'hswath',
    cmap = custom_cmap(15, 'YlOrRd', 0.15, 0.9, 'YlGnBu_r', 0, 0.85),
);

time_s = sp.time[np.nanargmax(sp.efth.spec.hs().values)]

Plot_SuperPoint_TC_Forecast(sp, time_s, ymax=5, figsize=[20,6]);

BinWaves#

BinWaves

Using the SuperPoint information from SHyTCWaves at a regional scale, we downscale the waves locally my means of BinWaves1.

out_sim_swan = xr.open_dataset(p_swan_output)

out_sim = cut_resample_area(
    out_sim_swan, area_extraction, 
    resample_factor, 
    case = 15,
    plot = False,
)

We interpolate the spectra to the same coordinates used to run BinWaves

coords_binwaves = xr.open_dataset(p_superpoint_binwaves).isel(time=0)

sp = sp.interp(
    dir = coords_binwaves.dir.values, 
    freq = coords_binwaves.freq.values,
)

# binwaves spec reconstruction
name_recon = 'binwaves_{0}_{1}.nc'.format(name.decode("utf-8"), year)
p_recon = op.join(p_tcforecast, name_recon)

if run_binwaves:
    sp_p = reconstruct_spec(p_kp_coeffs, sp, out_sim)
    sp_p.to_netcdf(p_recon)

else:
    sp_p = xr.open_dataset(p_recon)
    

Generating nearshore forecast...: 100%|██████████████████████████████| 1127/1127

# use wavespectra to calculate Hs, Tp, Tm and Dp
sp_p['Hs'] = sp_p.efth.spec.hs()
sp_p['Tp'] = sp_p.efth.spec.tp()
sp_p['Tm'] = sp_p.efth.spec.tm02()
sp_p['Dp'] = sp_p.efth.spec.dp()

Below are the maps for the time in which the SuperPoint energy was the highest

Plot_Grid_HsTmDp(sp_p, time_s);

SWATH#

Here we define a list of interest points to see the evolution of the wave height and direction over time

points_interest = np.array([
    [-172.23,-171.699, -172.675,-172.58],
    [-13.95, -13.794 , -13.47, -13.785],
])

points_names = ['P1', 'P2', 'P3', 'P4' ]

Plot_Grid_HsPd_series(sp, sp_p, points_interest, p_kp_coeffs, out_sim, points_names);

Generating nearshore forecast...: 100%|██████████████████████████████| 1/1 [00:0
Generating nearshore forecast...: 100%|██████████████████████████████| 1/1 [00:0
Generating nearshore forecast...: 100%|██████████████████████████████| 1/1 [00:0
Generating nearshore forecast...: 100%|██████████████████████████████| 1/1 [00:0

Astronomical Tide#

# domain coordinates
lon_tide_domain = [TWL.lon.min().values-0.5, TWL.lon.max().values+0.5]
lat_tide_domain = [TWL.lat.min().values-0.5, TWL.lat.max().values+0.5]

ds_tide = calculate_AT_TPXO9v4(
    model_directory, lon_tide_domain, lat_tide_domain,
    TWL.time.values[0], TWL.time.values[-1] + np.timedelta64(1, 'h'),
)

ds_tide['lon'] = ds_tide['lon'] + 360

tide = (
    (ds_tide.interpolate_na(dim=('lon'), method='nearest').tide.values) + \
    (ds_tide.interpolate_na(dim=('lat'), method='nearest').tide.values)
)/2

ds_tide['tide'] = (('lat','lon','time'), tide)

TWL['AT'] = (('time','lon', 'lat'), ds_tide.interp(
    lon = TWL.lon.values,
    lat = TWL.lat.values,
    method = 'linear',
).tide.transpose('time','lon','lat').values)

TWL['TWL'] = 0.3 * TWL['Hs'] + TWL['SS'] + TWL['AT']
TWL['TWL_Swath'] = (('lon', 'lat'), TWL.TWL.max(dim='time'))

TWL = TWL.sel(lon=slice(sp_p.lon.min(), sp_p.lon.max()), lat=slice(sp_p.lat.min(), sp_p.lat.max()))

TWL#

TWL

For calculating the total water level (TWL) proxy, we use the following formula:

\[ \eta _{TWL} = \eta _{AT} + \eta _{SS} + 0.3* H_s \]

where the astronomical tide (\(\eta _{AT}\)) is extracted from TOPEX, the storm surge (\(\eta _{SS}\)) from GreenSurge model and \(0.3* H_s\) is a proxy of the setup produced by waves, which are obtained from a combination of SHyTCWaves at a regional scale and downscaled using BinWaves

Below are the SWATH maps associated to the Storm Surge (Wind SetUp + Inverse Barometer), Hs, and TWL

Plot_Grid_SWATH(TWL);

Point Analysis

Below, any point from the ones defined above can be selected to see the evolution of the waves and different level parameters over time.

# point selection
p = 3
p_int = points_interest[0][p] + 360, points_interest[1][p]

twl = TWL.isel(
    lon = np.argmin(np.abs(TWL.lon.values-p_int[0])), 
    lat = np.argmin(np.abs(TWL.lat.values-p_int[1])),
)

mp = np.argmax(twl.TWL.values)

print('Point selected: {0}'.format(points_names[p]))
print('Coordinates: {0}'.format(p_int))

Point selected: P4
Coordinates: (187.42, -13.785)

point_code = points_names[p] + '   ' + str(p_int)

Plot_TWL_timeseries(twl, mp, point_code);

Tip

Below there is an animation showing the storm surge (\(\eta _{SS}\)), the astronomical tide (\(\eta _{AT}\)), the \(H_s\) and the combination of the different variables in the TWL.

Historical TC : Evan (2012)

Contents

Historical TC : Evan (2012)#

Database and site parameters#

TC Selection#

ShyTCWaves#

SWATH#

SuperPoint#

BinWaves#

SWATH#

GreenSurge#

Astronomical Tide#

TWL#