Storm Surge emulator
Contents
11. Storm Surge emulator#
inputs required:
Historical DWTs
Historical SLP
Synthetic DWTs
in this notebook:
Obtain historical storm surge as the inverse barometer
Simulation of storm surge for each synthetic DWT (daily resolution)
Workflow:
The offshore storm surge can be approximated as the inverse barometer effect: A decrease in atmospheric pressure of 1 mb will produce an increase in sea level of around 1 cm.
The storm surge generated by the wind set-up will only be considered at nearshore locations.
The synthetic generation of surge conditions on a regular climate (DWT1-36) is made by randomly selecting one of the historical surge conditions for each DWT.
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# common
import os
import os.path as op
import glob
# pip
import numpy as np
import xarray as xr
import pandas as pd
from datetime import datetime
import matplotlib.pyplot as plt
from matplotlib import gridspec
# DEV: override installed teslakit
import sys
sys.path.insert(0, op.join(os.path.abspath(''), '..', '..', '..'))
# teslakit
from bluemath.teslakit.io.aux_nc import StoreBugXdset
11.1. Files and paths#
# project path
p_data = '/Users/albacid/Projects/TeslaKit.2.0_projects/SAMOA'
p_out = op.join(p_data,'TIDE')
# input files
dwt_his_file = op.join(p_data,'ESTELA','pred_slp_grd','kma.nc')
dwt_sim_file = op.join(p_data,'ESTELA','DWT_sim.nc')
slp_his_file = op.join(p_data,'resources','SLP_hind_daily_wgrad.nc')
# output files
IB_his_file = op.join(p_out, 'IB_noTCs_daily_his.nc')
IB_sim_file = op.join(p_out, 'IB_noTCs_daily_sim.nc')
11.2. Parameters#
site = 'samoa'
#pnt_lon = -172.07+360 (187.93)
#pnt_lat = -13.76
lonp = 187 # closest coordinates in SLP dataset
latp = -14
SLP_his = xr.open_dataset(slp_his_file)
DWTs_his = xr.open_dataset(dwt_his_file)
DWTs_sim = xr.open_dataset(dwt_sim_file)
DWTs_his = xr.Dataset({'bmus':(('time'), DWTs_his.sorted_bmus_storms.values.astype(int))},
coords = {'time':(('time'), DWTs_his.time.values)})
n_clusters = 42
11.3. Obtain historical IB#
SLP_his = SLP_his.sel(time=slice(DWTs_his.time.values[0], DWTs_his.time.values[-1]))
# --------------------------------------
# Inverse Barometer (IB)
# select closest grid point to Site
IB_his = SLP_his.sel(longitude = lonp, latitude = latp)
# Calculate anomalies and change units
IB_his['slp'] = IB_his.slp*0.01 # (Pa to mb)
IB_his['slp_anomaly'] = IB_his.slp - np.mean(IB_his.slp.values)
IB_his['level_IB'] = -1*IB_his.slp_anomaly # (Inverse Barometer: mb to cm)
IB_his['level_IB'] = IB_his.level_IB/100.0 # (cm to m)
# Add DWT to dataset
IB_his['DWT'] = DWTs_his.bmus
IB_his = IB_his.drop({'slp','slp_gradient','slp_anomaly' })
# save
IB_his.to_netcdf(IB_his_file)
IB_his = xr.open_dataset(IB_his_file)
11.4. IB emulator#
# Define output
IB_sim = DWTs_sim.copy(deep=True)
IB_sim['level_IB']=(('time','n_sim'), np.zeros((len(DWTs_sim.time), len(DWTs_sim.n_sim) ))*np.nan)
# Group historical IB data by DWT
IB_his_dwts = np.zeros((n_clusters,len(IB_his.time)))*np.nan
for dd in range(n_clusters):
IB_his_dd = IB_his.where(IB_his.DWT==dd, drop=True)
IB_his_dwts[dd,:len(IB_his_dd.time)] = IB_his_dd.level_IB.values
# Simulate
for n in DWTs_sim.n_sim:
print('- Sim: {0} -'.format(int(n)+1))
# Select simulation of DWTs
DWTs_sim_n = DWTs_sim.isel(n_sim=n)
for ind_t, dwt in enumerate(DWTs_sim_n.evbmus_sims.values):
# select randomly one historical value for that DWT (including TCs in case it doesn't enter r2)
IB_his_dwt = IB_his_dwts[int(dwt-1),:]
IB_his_dwt = IB_his_dwt[~np.isnan(IB_his_dwt)]
rnd_t = np.random.randint(0, len(IB_his_dwt))
IB_sim['level_IB'][ind_t,n] = IB_his_dwt[rnd_t]
# save
StoreBugXdset(IB_sim, IB_sim_file)
11.5. Validation#
IB_sim = xr.open_dataset(IB_sim_file)
# choose simulation to plot
n = 0
IB_sim_n = IB_sim.isel(n_sim=n)
# All DWTs
plt.figure(figsize=(5,5))
n_bins = np.linspace(-0.1,.3,20)
plt.hist(IB_his.level_IB, n_bins,density=True,label='Historical',color='salmon',alpha=0.5,edgecolor='black');
plt.hist(IB_sim_n.level_IB, n_bins,density=True,label='Emulator',color='skyblue',alpha=0.5,edgecolor='black');
plt.legend()
plt.title('IB (cm) from all DWTs');
![_images/11_IB_Emulator_14_0.png](_images/11_IB_Emulator_14_0.png)
import warnings
warnings.filterwarnings('ignore')
# For each DWT
n_bins = np.linspace(-0.1,.3,20)
fig = plt.figure(figsize=[15.6,11.4])
gs = gridspec.GridSpec(7, 6, hspace=0.5) #, wspace=0.0, hspace=0.0)
gr, gc = 0, 0
for dwt in range(0,36):
IB_his_dwt = IB_his.level_IB.where(IB_his.DWT==dwt, drop=True)
IB_sim_n_dwt = IB_sim_n.level_IB.where(IB_sim_n.evbmus_sims==dwt+1, drop=True)
ax = plt.subplot(gs[gr, gc])
ax.hist(IB_his_dwt, n_bins,density=True,label='Historical',color='salmon',alpha=0.5,edgecolor='black');
ax.hist(IB_sim_n_dwt, n_bins,density=True,label='Emulator',color='skyblue',alpha=0.5,edgecolor='black');
ax.legend(prop={'size':6})
ax.set_title('IB (cm) from DWT ' + str(dwt+1) + '. ' + str(len(IB_his_dwt)) + ' data', fontsize=8)
ax.set_xticks([-.1,0,.1,.2,.3])
plt.xticks(fontsize=7)
# counter
gc += 1
if gc >= 6:
gc = 0
gr += 1
![_images/11_IB_Emulator_15_0.png](_images/11_IB_Emulator_15_0.png)