Add TCs to synthetic data
Contents
13. Add TCs to synthetic data#
Results of simulated historical & synthetic TCs (SHyTCWaves) are added to the output of the waves and surges emulator
inputs required:
Synthetic DWTs
Synthetic storm surge
Synthetic spectra
Synthetic bulk parameters
Track and pressures of historical and synthetic TCs simulated with SHyTCWaves
Super-point of historical and synthetic TCs simulated with SHyTCWaves
in this notebook:
Add surge (inverse barometer) from TCs whenever the simulated TC enters the 7.5 degrees radio.
Add waves spectrum from TCs whenever the simulated TC enters the 7.5 degrees radio.
Correct the Hs (Hs_corrected = sqrt(offset^2 + alfa^2*Hs^2))
The synthetic generation of surge conditions and wave conditions on a TC climate (DWT37-42) is made by randomly selecting a TC that affects the site on a 14º radio and checking if the TC also affects the site on a 7.5º radio; if that is the case, the storm surge (inverse barometer) and the super-point spectrum generated from the random TC will be incorporated, replacing the surge for the day when the TC is closest to site, and adding the TC spectrum to the regular climate spectrum for the complete duration of the TC.
#!/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
import string
import pickle5 as pk
import wavespectra
# DEV: override installed teslakit
import sys
sys.path.insert(0, op.join(os.path.abspath(''), '..', '..', '..'))
# teslakit
from bluemath.teslakit.io.aux_nc import StoreBugXdset
from bluemath.teslakit.waves_emulator import add_TC_to_emulator
Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.20.0
/Users/albacid/Software/Bitbucket_repos/teslakit.2.0/venv/lib/python3.7/site-packages/wavespectra/specarray.py:29: UserWarning: Cannot import sympy, install "extra" dependencies for full functionality
'Cannot import sympy, install "extra" dependencies for full functionality'
13.1. Files and paths#
# project path
p_data = '/Users/albacid/Projects/TeslaKit.2.0_projects/SAMOA'
site = 'samoa'
p_in = op.join(p_data,'SIMULATION','climate_emulator','WAVES_noTCs' )
p_in_spec = op.join(p_in,'spec' )
p_out = op.join(p_data,'SIMULATION','climate_emulator')
p_out_bulk = op.join(p_out,'WAVES_TCs')
p_out_spec = op.join(p_out,'WAVES_TCs','spec')
if not os.path.isdir(p_out_bulk): os.makedirs(p_out_bulk)
if not os.path.isdir(p_out_spec): os.makedirs(p_out_spec)
# waves from historical and synthetic TCs in r2
p_TCs = '/Volumes/WD/SHyTCWaves_' + site
# input files
IB_noTCs_file = op.join(p_data, 'TIDE','IB_noTCs_daily_sim.nc')
dwt_sim_file = op.join(p_data,'ESTELA','DWT_sim.nc')
tcs_sim_r1_params_file = op.join(p_data,'TCs','TCs_r1_sim_params_' + site + '.nc') # Nakajo TCs parameters r1
tcs_hist_r1_params_file_sel = op.join(p_data,'TCs','TCs_hist_r1_params_' + site + '_raw.nc') # historical TCs parameters r1
tcs_r2_file = op.join(p_data, 'resources','TCs_ShyTCwaves','indices_histnak_samoa_r8.nc') # TCs in r2 (shyTCwaves)
# output files
IB_TCs_file = op.join(p_data, 'TIDE','IB_TCs_daily_sim.nc')
rnd_TCs_file = op.join(p_out,'rnd_TCs.nc')
rnd_TCs_textfile = op.join(p_out,'rnd_TCs_' + site + '.txt')
13.2. Parameters#
# Load sim DWTs
DWTs_sim = xr.open_dataset(dwt_sim_file)
# Load Surge
IB_noTCs = xr.open_dataset(IB_noTCs_file)
# Load historical and synthetic TCs from r1
TCs_his = xr.open_dataset(tcs_hist_r1_params_file_sel)
TCs_sim = xr.open_dataset(tcs_sim_r1_params_file)
# Join TCs from r1 in one file
TCs_his['id'] = (('storm'),np.arange(0, len(TCs_his.storm)))
TCs_his = TCs_his[['category','id']]
TCs_his = TCs_his.assign_coords(storm = np.arange(0, len(TCs_his.storm)))
TCs_his['type'] = (('storm'), np.zeros((len(TCs_his.storm))))
TCs_his['type'][:] = 1
TCs_sim['id'] = (('storm'), np.arange(0, len(TCs_sim.storm)))
TCs_sim = TCs_sim[['category','id']]
TCs_sim = TCs_sim.assign_coords(storm = np.arange(len(TCs_his.storm), len(TCs_his.storm)+len(TCs_sim.storm)))
TCs_sim['type'] = (('storm'), np.zeros((len(TCs_sim.storm))))
TCs_sim['type'][:] = 2
TCs = xr.concat([TCs_his, TCs_sim], dim='storm')
TCs['type'].attrs = {'0':'no TC', '1':'historic', '2':'synthetic'}
# Load TCs from r2 (simulated TCs, shytcwaves)
tcs_r2 = xr.open_dataset(tcs_r2_file)
13.3. Select random TCs#
# define output
rnd_TCs = IB_noTCs.copy(deep=True)
rnd_TCs = rnd_TCs.drop({'level_IB','evbmus_sims'})
rnd_TCs['TC_cat'] = (('time','n_sim'), np.full((len(IB_noTCs.time), len(IB_noTCs.n_sim)), np.nan))
rnd_TCs['TC_type'] = (('time','n_sim'), np.zeros((len(IB_noTCs.time), len(IB_noTCs.n_sim))))
rnd_TCs['TC_id'] = (('time','n_sim'), np.full((len(IB_noTCs.time), len(IB_noTCs.n_sim)), np.nan))
rnd_TCs['TC_type'].attrs = {'0':'no TC', '1':'historic', '2':'synthetic'}
rnd_TCs = rnd_TCs.drop({'level_IB','evbmus_sims'})
# Select random TC when bmus from TC (37-42)
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)
# data for each day
for ind_t, dwt in enumerate(DWTs_sim_n.evbmus_sims.values):
if dwt>36:
cat_TC = int(dwt -36 -1)
#------------------------------
# Select random TC for that category (between historical and simulated) in r1
TC_cat = TCs.where(TCs.category==cat_TC, drop=True)
rnd_TC = np.random.randint(0,len(TC_cat.storm))
TC_cat = TC_cat.isel(storm=rnd_TC)
# Check if r1 TC (TCs for obtaining DWTs) enters in r2 (simulated TCs)
if TC_cat.type==1:
tcs_r2_r1 = tcs_r2[['id_hist','cat_hist','cat_hist_bad']]
tcs_r2_r1 = tcs_r2_r1.where(tcs_r2_r1.id_hist==TC_cat.id,drop=True)
tc_in_r2 = len(tcs_r2_r1.hist)
elif TC_cat.type==2:
tcs_r2_r1 = tcs_r2[['id_syn','cat_syn','cat_syn_bad']]
tcs_r2_r1 = tcs_r2_r1.where(tcs_r2_r1.id_syn==TC_cat.id,drop=True)
tc_in_r2 = len(tcs_r2_r1.syn)
if tc_in_r2!=0 :
# Save selected TC indicators
rnd_TCs['TC_cat'][ind_t, n] = TC_cat.category.values
rnd_TCs['TC_type'][ind_t, n] = TC_cat.type.values
rnd_TCs['TC_id'][ind_t, n] = TC_cat.id.values
# save
###rnd_TCs.to_netcdf(rnd_TCs_file)
# select synthetic TCs to simulate with Green Surge
rnd_TCs = rnd_TCs.where(rnd_TCs.TC_type==2, drop=True)
TC_id = rnd_TCs.TC_id.values.flatten()
TC_id = TC_id[~np.isnan(TC_id)]
TC_id = np.unique(TC_id)
np.savetxt(rnd_TCs_textfile, TC_id, '%d',delimiter='\n')
# Load random selected TCs
rnd_TCs = xr.open_dataset(rnd_TCs_file)
13.4. Add TCs Surge#
# define output dataset
IB_TCs = IB_noTCs.copy(deep=True)
IB_TCs['TC_cat'] = (('time','n_sim'), np.full((len(IB_noTCs.time), len(IB_noTCs.n_sim)), np.nan))
IB_TCs['TC_type'] = (('time','n_sim'), np.zeros((len(IB_noTCs.time), len(IB_noTCs.n_sim))))
IB_TCs['TC_id'] = (('time','n_sim'), np.full((len(IB_noTCs.time), len(IB_noTCs.n_sim)), np.nan))
IB_TCs['TC_type'].attrs = {'0':'no TC', '1':'historic', '2':'synthetic'}
# Use randomly selected TCs from above
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)
# select random selected TCs for that simulation
rnd_TCs_n = rnd_TCs.isel(n_sim=n)
rnd_TCs_n = rnd_TCs_n.where(~np.isnan(rnd_TCs_n.TC_cat), drop=True)
# data for each day
for ind_t, (t, dwt) in enumerate(zip(DWTs_sim_n.time.values, DWTs_sim_n.evbmus_sims.values)):
if dwt>36:
cat_TC = int(dwt -36 -1)
# There was a TC in r1, check if entered r2
TC_cat = rnd_TCs_n.where(rnd_TCs_n.time==t, drop=True)
if len(TC_cat.time)!=0:
# Load selected random TC track
path_rnd_TC = op.join(p_TCs, 'shy_tesla_' + site + '_' + TC_cat.TC_type.attrs[str(int(TC_cat.TC_type.values))],)
file = glob.glob(op.join(path_rnd_TC,'track_samoa_cat*_id' + str(int(TC_cat.TC_id.values)) + '_' + TC_cat.TC_type.attrs[str(int(TC_cat.TC_type.values))] + '.pkl'))
trackTC = pk.load(open(file[0],'rb'))
# Add IB from when TC passes closest to site
min_dist = np.argmin(trackTC.dist2land.values)
press_min_dist = trackTC.p0[min_dist]
IB_TC_mindist = press_min_dist - 1013
IB_TCs['level_IB'][ind_t, n] = (-1*IB_TC_mindist)/100.0
# Save selected TC indicators
IB_TCs['TC_cat'][ind_t, n] = TC_cat.TC_cat.values[0]
IB_TCs['TC_type'][ind_t, n] = TC_cat.TC_type.values[0]
IB_TCs['TC_id'][ind_t, n] = TC_cat.TC_id.values[0]
# save
IB_TCs.to_netcdf(IB_TCs_file)
# Plot Inverse barometer of regular climate and TCs
surge = xr.open_dataset(IB_noTCs_file)
surge_TCs = xr.open_dataset(IB_TCs_file)
surge_n = surge.isel(n_sim=0)
surge_TCs_n = surge_TCs.isel(n_sim=0)
plt.figure(figsize=(20,10))
plt.plot(surge_TCs_n.time, surge_TCs_n['level_IB'], 'r',label = 'IB with TCs')
plt.plot(surge_n.time, surge_n['level_IB'], label = 'IB')
plt.legend()
plt.grid()
#plt.xlim(datetime(2032,2,15), datetime(2032,3,15))
13.5. Add TCs wave spectrum#
for n in DWTs_sim.n_sim:
print()
print('- Sim: {0} -'.format(int(n)+1))
# Select simulation of DWTs
DWTs_sim_n = DWTs_sim.isel(n_sim=n)
# select random selected TCs for that simulation
rnd_TCs_n = rnd_TCs.isel(n_sim=n)
rnd_TCs_n = rnd_TCs_n.where(~np.isnan(rnd_TCs_n.TC_cat), drop=True)
# data for each year
years = DWTs_sim_n.time.dt.year.values
for yy in np.arange(years[0],years[-1]):
y_ini = str(yy) + '-01-01'
y_end = str(yy) + '-12-31'
DWTs_sim_n_yy = DWTs_sim_n.sel(time=slice(y_ini,y_end))
# load simulated spectrum from regular climate
spec = xr.open_dataset(op.join(p_in_spec, 'spectrum_hourly_{0:08d}_{1}.nc'.format(int(n),yy)))
# variables to save when DWT from TCs
spec['TC_cat'] = (('time'), np.full((len(spec.time)), np.nan))
spec['TC_type'] = (('time'), np.zeros((len(spec.time))))
spec['TC_id'] = (('time'), np.full((len(spec.time)), np.nan))
spec['TC_type'].attrs = {'0':'no TC', '1':'historic', '2':'synthetic'}
# data for each day
for t, dwt in zip(DWTs_sim_n_yy.time.values, DWTs_sim_n_yy.evbmus_sims.values):
if dwt>36:
# There was a TC in r1, check if entered r2
TC_cat = rnd_TCs_n.where(rnd_TCs_n.time==t, drop=True)
if len(TC_cat.time)!=0:
path_rnd_TC = op.join(p_TCs, 'shy_tesla_' + site + '_' + TC_cat.TC_type.attrs[str(int(TC_cat.TC_type.values))],)
file_TC = glob.glob(op.join(path_rnd_TC,'superpoint_samoa_cat*_id' + str(int(TC_cat.TC_id.values)) + '_' + TC_cat.TC_type.attrs[str(int(TC_cat.TC_type.values))] + '.nc'))
# Add TC spectrum to regular climate spectrum
spec_TC = xr.open_dataset(file_TC[0])
spec = add_TC_to_emulator(t, spec_TC, spec, TC_cat)
else:
# Add spectrum for last day
spec_t = spec.efth.sel(time=t)
spec_t = spec.efth.where(spec.time==spec_t.time)
spec_t = spec_t.isel(dir=0,freq=0)
ind_t = np.where(~np.isnan(spec_t))[0][0]
# add spectrum for last day
if ind_t-24<0: # inicio de año
print('dwt tc al inicio')
elif ind_t+24>len(spec.time): # fin año
print('dwt tc al final')
else:
# add spectrum from last day
spec['efth'][:,:, ind_t:ind_t+24] = spec.efth[:,:, ind_t:ind_t+24].values + spec.efth[:,:, ind_t-24:ind_t].values
#------------------------------
# save spectrum
StoreBugXdset(spec, op.join(p_out_spec, 'spectrum_hourly_{0:08d}_{1}.nc'.format(int(n),yy)))
#------------------------------
# compute bulk params
stats = spec.spec.stats(['hs','tp','tm02','dpm','dspr'])
# add TC parameters
stats['TC_cat'] = (('time'), spec['TC_cat'].values)
stats['TC_id'] = (('time'), spec['TC_id'].values)
stats['TC_type'] = (('time'), spec['TC_type'].values)
stats['TC_type'].attrs = {'0':'no TC', '1':'historic', '2':'synthetic'}
if yy == years[0]:
bulk_TCs = stats.copy(deep=True)
else:
bulk_TCs = xr.concat([bulk_TCs, stats], dim='time')
#------------------------------
# save bulk
StoreBugXdset(bulk_TCs, op.join(p_out_bulk, 'bulk_params_hourly_{0:08d}.nc'.format(int(n))))
13.5.1. Correct Hs#
# plot Hs from regular climate and Hs with TCs added
n = 6 # choose sim
bulk = xr.open_dataset(op.join(p_data,'SIMULATION','climate_emulator','WAVES_noTCs', 'bulk_params_hourly_{0:08d}.nc'.format(int(n))))
bulk_TCs = xr.open_dataset(op.join(p_out_bulk, 'bulk_params_hourly_{0:08d}.nc'.format(int(n))))
plt.figure(figsize=(20,10))
plt.plot(bulk_TCs.time, bulk_TCs['hs'], '.-',label = 'hs with TCs')
plt.plot(bulk.time, bulk['hs'], label = 'hs')
plt.legend()
plt.grid()
# Correct Hs of synthetic TCs
alfa_1 = 1.6
alfa_2 = 1.3
alfa_3 = 1.1
offset = 0 # no offset since TCs are added to regular waves not replaced as in the hindcast
for n in DWTs_sim.n_sim.values:
print(n)
bulk_TCs = xr.open_dataset(op.join(p_out_bulk, 'bulk_params_hourly_{0:08d}.nc'.format(int(n))))
bulk_TCs['hs_corrected'] = (('time'), bulk_TCs.hs.copy(deep=True))
ind = np.where(bulk_TCs.TC_cat<=3)
bulk_TCs['hs_corrected'][ind] = np.sqrt(offset**2 + alfa_1**2*bulk_TCs['hs'].values[ind]**2)
ind = np.where(bulk_TCs.TC_cat==4)
bulk_TCs['hs_corrected'][ind] = np.sqrt(offset**2 + alfa_2**2*bulk_TCs['hs'].values[ind]**2)
ind = np.where(bulk_TCs.TC_cat==5)
bulk_TCs['hs_corrected'][ind] = np.sqrt(offset**2 + alfa_3**2*bulk_TCs['hs'].values[ind]**2)
StoreBugXdset(bulk_TCs, op.join(p_out_bulk, 'bulk_params_hourly_{0:08d}.nc'.format(int(n))))
# plot Hs from regular climate and Hs with TCs added (corrected)
n = 6 # choose sim
bulk = xr.open_dataset(op.join(p_data,'SIMULATION','climate_emulator','WAVES_noTCs', 'bulk_params_hourly_{0:08d}.nc'.format(int(n))))
bulk_TCs = xr.open_dataset(op.join(p_out_bulk, 'bulk_params_hourly_{0:08d}.nc'.format(int(n))))
plt.figure(figsize=(20,10))
plt.plot(bulk_TCs.time, bulk_TCs['hs_corrected'], '.-',label = 'hs with TCs corrected')
plt.plot(bulk.time, bulk['hs'], label = 'hs')
plt.legend()
plt.grid()
13.5.2. Correct spectrum#
# Correct spectrum of synthetic TCs
alfa_1 = 1.6
alfa_2 = 1.3
alfa_3 = 1.1
offset = 0 # no offset since TCs are added to regular waves not replaced as in the hindcast
years = DWTs_sim.isel(n_sim=0).time.dt.year.values
for n in DWTs_sim.n_sim[6:7]:
print(n.values)
for yy in np.arange(years[0],years[-1]):
spec_TCs = xr.open_dataset(op.join(p_out_spec, 'spectrum_hourly_{0:08d}_{1}.nc'.format(int(n),yy)))
spec_TCs['efth_corrected'] = (('dir','freq','time'), spec_TCs.efth.copy(deep=True))
ind = np.where(spec_TCs.TC_cat<=3)[0]
spec_TCs['efth_corrected'][:,:,ind] = offset + alfa_1**2*spec_TCs['efth'].values[:,:,ind]
ind = np.where(spec_TCs.TC_cat==4)[0]
spec_TCs['efth_corrected'][:,:,ind] = offset + alfa_2**2*spec_TCs['efth'].values[:,:,ind]
ind = np.where(spec_TCs.TC_cat==5)[0]
spec_TCs['efth_corrected'][:,:,ind] = offset + alfa_3**2*spec_TCs['efth'].values[:,:,ind]
StoreBugXdset(spec_TCs, op.join(p_out_spec, 'spectrum_hourly_{0:08d}_{1}.nc'.format(int(n),yy)))