TCs Winds
Contents
TCs Winds#
# TODO:
# - codigos CRS.from_epsg(4326) independientes de site?
This way, we have created a library of cases, with more than 7000 events, that allow for any new forecasted event, and based on the tracks from the Joint Typhon Warning Center, (JTWC), the winds that these tracks might produce are simulated with Holland’s empirical model at the ERA-5 resolution, then these forecasted winds are compared with the library of cases, choosing the most similar ones to reconstruct the high resolution winds based on k-nearest neighbors.
Forecast
In this section we will produce the probabilistic assessment of WINDS for the forecast information of the TC GITA issued some days previous to its arrival
from datetime import datetime
print("execution start: {0}".format(datetime.today().strftime('%Y-%m-%d %H:%M:%S')))
execution start: 2022-09-03 00:42:53
import os
import os.path as op
import sys
import time
import pandas as pd
import numpy as np
import xarray as xr
from pyproj import CRS
from sklearn.cluster import KMeans
from scipy.stats import qmc
from shapely.geometry import Point, Polygon
from matplotlib.colors import ListedColormap
from cartopy.feature import ShapelyFeature
from cartopy.io.shapereader import Reader
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# raster tools
from rasterio.crs import CRS
import rioxarray as rio
# bluemath
sys.path.insert(0, op.join(op.abspath(''), '..', '..'))
sys.path.insert(0, op.join(op.abspath(''), '..'))
from bluemath.tc_warnings import download_warning
from bluemath.tc_forecast import get_info_TCs_warning
from bluemath.tc_forecast import Plot_Tracks_Maxwinds, Plot_Tracks_Maxwinds_Density, Plot_Tracks_Maxwinds_zoom
from bluemath.winds.graffiti_rain import *
from bluemath.winds.plots.common import custom_cmap
from bluemath.winds.tc_forecast import Load_TCs_Warning, Plot_tracks, Plot_Density_Forecast_Tracks
from bluemath.winds.winds import windy
from bluemath.winds.plots.probs import plot_prob_assessment
# operational utils
from operational.util import clean_forecast_folder, read_config_args
Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.16.0
# hide warnings
import warnings
warnings.filterwarnings('ignore')
# HTML recipe to center notebook plots
from IPython.core.display import HTML
HTML("""
<style>
.output_png {
display: table-cell;
text-align: center;
vertical-align: middle;
}
</style>
""")
# matplotlib animatiom embed limit
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 2**128
# prepare listed colormap
cmap_probs = ListedColormap([
'forestgreen','limegreen','chartreuse','yellow','gold',
'goldenrod','darkorange','crimson','brown','darkmagenta',
])
JTWC Tropical Warnings#
https://www.metoc.navy.mil/jtwc/jtwc.html
Forecast Parameters#
# project database
p_data = r'/media/administrador/HD2/SamoaTonga/data'
p_db = r'/media/administrador/HD1/DATABASES'
#site = 'Savaii' # Savaii / Upolu / Tongatapu
# (optional) get site from config file
nb_args = read_config_args(op.abspath(''), '10a_wind_damage')
site = nb_args['site']
print('Study site: {0}'.format(site))
Study site: Upolu
# current date (for execution folder)
date = datetime.today().strftime('%Y%m%d%H%M')
print('forecast execution date: {0}'.format(date))
forecast execution date: 202209030042
# forecast parameteres
_clean_forecast_folder_days = 7 # days to keep after forecast folder cleaning
# choose agencies for TC warning
ags_sel = [
' UKMI', ' JMA', ' JGSI', ' AVNI', ' HWFI', ' NVGI',
' ECM2', ' JTWC', ' JNVI', ' NVG2', ' UEMI', ' UHM2', ' UKM2',
]
# site related parameteres
if site == 'Savaii':
site_center = (-172.40+360,-13.65)
# area of interest
area_interest = [[177, -165+360] , [-21.5, -7]]
elif site == 'Upolu':
site_center = (-171.78+360,-13.92)
# area of interest
area_interest = [[177, -165+360] , [-21.5, -7]]
elif site == 'Tongatapu':
site_center = (-175.19+360,-21.16)
# area of interest
area_interest = [[179, -169+360] , [-23.6, -18.1]]
area_interest_plot = [x for sl in area_interest for x in sl]
# point shape
radius = 0.15
site_shape = Point(*site_center).buffer(radius)
Database#
p_site = op.join(p_data, site)
# deliverable folder
p_deliv = op.join(p_site, 'd10_wind_damage')
# WRF database
p_db_winds = op.join(p_db, 'WRF_Wind')
# ERA5
p_db_era = op.join(p_db_winds, 'era5_data', 'era5_data_{0}.nc'.format(site))
# high resolution
p_db_hr = op.join(p_db_winds, 'hr_wind_data', 'HR_{0}_ALL.nc'.format(site))
# coast path
p_coast = op.join(p_data, 'gshhg-shp-2.3.7', 'GSHHS_shp', 'f', 'GSHHS_f_L1.shp')
# load the high-resolution coastline for the plots
coast_full = ShapelyFeature(
Reader(p_coast).geometries(),
ccrs.PlateCarree(), edgecolor='black',
)
Forecast Output Folder#
p_forecast = op.join(p_site, 'forecast', '10_wind_damage')
# forecast folder
p_fore_date = op.join(p_forecast, date)
print('forecast date code: {0}'.format(date))
# prepare forecast subfolders for this date
p_fore_warn = op.join(p_fore_date, 'warnings')
for p in [p_fore_date, p_fore_warn]:
if not op.isdir(p):
os.makedirs(p)
# clean forecast folder
clean_forecast_folder(p_forecast, days=_clean_forecast_folder_days)
forecast date code: 202209030042
deleted folder: /media/administrador/HD2/SamoaTonga/data/Upolu/forecast/10_wind_damage/202208140720
deleted folder: /media/administrador/HD2/SamoaTonga/data/Upolu/forecast/10_wind_damage/202208150720
Download TC warnings#
# download TC warnings
files_ash, files_bsh = download_warning(p_fore_warn, code_basin='sh')
At: 2022-09-02T22:43:19
Total number of active storms: 6
Total number of active storms in the sh: 0
# stop forecast if no warnings
if len(files_bsh) == 0:
raise ValueError('No warnings in the SH. Forecast execution cancelled.')
sys.exit()
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
/tmp/ipykernel_34885/2713157806.py in <module>
1 # stop forecast if no warnings
2 if len(files_bsh) == 0:
----> 3 raise ValueError('No warnings in the SH. Forecast execution cancelled.')
4 sys.exit()
5
ValueError: No warnings in the SH. Forecast execution cancelled.
# choose storm code
code = 0
file_active = files_ash[code]
file_best = files_bsh[code]
print('active TC: {0}'.format(op.basename(file_active)))
print('best track: {0}'.format(op.basename(file_best)))
# show available warning times
name_warning, available_times_warning = get_info_TCs_warning(file_active)
print('name warning: {0}'.format(name_warning))
print('last available warning time: {0}'.format(available_times_warning[-1]))
# use last available time warning
time_warning = str(available_times_warning[-1])
print('time warning: {0}'.format(time_warning))
TC Warning Output#
tc_name = str(name_warning)
p_fore_tc = op.join(p_fore_date, tc_name + '_forecast')
p_hr_tc = op.join(p_fore_tc, 'hr_{0}_tracks.nc'.format(tc_name))
p_vortex_tc = op.join(p_fore_tc, 'vortex_{0}_tracks.nc'.format(tc_name))
p_final_track_tif = op.join(p_fore_tc, 'hr_FINAL_track_{0}_lonlat.tif'.format(tc_name))
p_final_track_tif_utm = op.join(p_fore_tc, 'hr_FINAL_track_{0}_utm.tif'.format(tc_name))
#p_final_track_tif_utm_risk = op.join(p_data, 'riskscape_projects', site, 'data', 'hr_FINAL_track_{0}_utm_risk.tif'.format(tc_name))
# warnings
p_fore_save_track = op.join(p_fore_date, 'Forecast_{0}_time_{1}.nc'.format(name_warning, time_warning[:-3]))
p_fore_save_track_sel = op.join(p_fore_date, 'Forecast_{0}_time_{1}_selected.nc'.format(name_warning, time_warning[:-3]))
#p_aux_f = op.join(p_data, 'riskscape_projects', site, 'data')
for p in [p_fore_tc]: #, p_aux_f]:
if not op.isdir(p):
os.makedirs(p)
# Load TRACKS
TRACKS = Load_TCs_Warning(file_active, time = time_warning)
# save warning TRACKS
TRACKS.to_netcdf(p_fore_save_track)
TRACKS = xr.open_dataset(p_fore_save_track)
Area#
Figure Explanation
The figures below show TC warning tracks at the study site.
Need visual inspection on wether the tracks are inside the area of interest.
If there are no tracks in this area, don’t continue running the notebook.
Plot_Tracks_Maxwinds(TRACKS);
Plot_Tracks_Maxwinds(TRACKS, area_extent=area_interest_plot);
Plot_Tracks_Maxwinds_Density(TRACKS, area_extent=area_interest_plot);
TC Tracks#
Track warning predictions
For a given historic track, at any given day, the warning track predictions issued by different available agencies are collected!!
Agencies Selection#
print('Number of tracks: {0}'.format(len(TRACKS.track.values)))
print('Available Agencies:\n{0}\n------------'.format(TRACKS.Agency.values))
c, ia, ib = np.intersect1d(TRACKS.Agency.values, ags_sel, return_indices=True)
print('Total number of tracks: {0}'.format(len(ia)))
print('Agencies from list available : ' + str(np.array(ags_sel)[ib]))
TRACKS = TRACKS.isel(track=ia)
TRACKS['Index'] = (('track'), ia)
if TRACKS.track.size > 0:
Plot_Tracks_Maxwinds_Density(TRACKS, area_extent=area_interest_plot);
if TRACKS.track.size > 0:
Plot_Tracks_Maxwinds_zoom(TRACKS, area_interest_plot);
fig = plt.figure(figsize=[30,15]) #30,20
gs = gridspec.GridSpec(1,2,hspace=0.01,wspace=0.01)
ax1 = fig.add_subplot(gs[0], projection = ccrs.PlateCarree(central_longitude=190))
_,_,area = Plot_tracks(TRACKS, var='maxwinds', ax=ax1, fig=fig);
ax2 = fig.add_subplot(gs[1], projection = ccrs.PlateCarree(central_longitude=190))
Plot_Density_Forecast_Tracks(TRACKS, ax=ax2, fig=fig);
fig = plt.figure(figsize=[30,15])
gs = gridspec.GridSpec(1,2,hspace=0.01,wspace=0.01)
ax1 = fig.add_subplot(gs[0],projection = ccrs.PlateCarree(central_longitude=190))
_,_,area = Plot_tracks(TRACKS, var='maxwinds', ax=ax1, fig=fig, area_extent=area_interest_plot)
ax2 = fig.add_subplot(gs[1],projection = ccrs.PlateCarree(central_longitude=190))
Plot_Density_Forecast_Tracks(TRACKS, ax=ax2, fig=fig, area_extent=area_interest_plot, m_size=100)
Winds prediction#
Figure Explanation
The figures below show TC max. winds prediction for each tracks at the study site.
# Plot each track maxwinds
p = TRACKS.plot.scatter(
x='longitude',y='latitude',
hue='maxwinds',col='track',
col_wrap=3,
subplot_kws={'projection':ccrs.PlateCarree(central_longitude=180)},
transform=ccrs.PlateCarree(),
)
# add point and set extent
for ax in p.axes.flat:
ax.scatter(*site_center, s=100,c ='red',transform=ccrs.PlateCarree())
ax.coastlines('10m')
ax.stock_img()
ext_e = 20 # extra area
ax.set_extent([site_center[0]-ext_e, site_center[0]+ext_e, site_center[1]-ext_e, site_center[1]+ext_e])
hr_data_recon_tracks = []
vortex_models = []
for track in range(len(TRACKS.track.values)):
try:
wind1 = windy(
data_path = p_db_winds,
site_name = site,
site_center = site_center,
site_shape = site_shape,
tropical_cyclon = TRACKS.isel(track=track),
tropical_cyclon_name = f'TRACK {track}',
tc_type = 'forecast',
era5_data = xr.open_dataset(p_db_era),
hr_wind_data = xr.open_dataset(p_db_hr),
postprocess = True,
check_correlation = (False,0.8),
verbose = True,
plot = True,
)
vortex_models.append(wind1.vortex_model.expand_dims({'track':[track]}))
closest_era5_data, hr_data_recon, era5_data_used, scalers_data_used = wind1.calculate_nearest_era5(
back_cells = 2, clos_times = 7,
scale = 'custommean', min_M_era5 = 0,
minmax = (0,5),
)
hr_data_recon_tracks.append(hr_data_recon.expand_dims({'track':[track]}))
print(f'TRACK {track} analyzed!!')
except:
print(f'TRACK {track} NOT analyzed!!')
hr_data_recon_all_tracks = xr.concat(hr_data_recon_tracks, dim='track').M.max(dim='time')
hr_data_recon_all_tracks.to_netcdf(p_hr_tc)
vortex_data_recon_all_tracks = xr.concat(vortex_models, dim='track').max(dim='time')
vortex_data_recon_all_tracks.to_netcdf(p_vortex_tc)
hr_data_recon_all_tracks = xr.open_dataarray(p_hr_tc)
vortex_data_recon_all_tracks = xr.open_dataset(p_vortex_tc)
Winds pattern for each predicted Track#
Figure Explanation
TODO: FIGURE EXPLANATION.
# plot winds module
p = np.sqrt(vortex_data_recon_all_tracks.u**2 + vortex_data_recon_all_tracks.v**2).plot(
col = 'track',
col_wrap = 5,
vmin = 10,
vmax = 50,
cmap = custom_cmap(100, 'plasma_r', 0.05,0.9, 'viridis', 0.2, 1),
subplot_kws = {'projection': ccrs.PlateCarree(central_longitude=180)},
transform = ccrs.PlateCarree(),
)
# add coast, set extent
for ax in p.axes.flat:
ax.add_feature(coast_full, facecolor='None', alpha=0.8, edgecolor='black', linewidth=1.5, zorder=5)
ax.set_extent((180,190,-24,-18))
# plot hr data
p = hr_data_recon_all_tracks.plot(
col = 'track',
col_wrap = 5,
vmin = 10,
vmax = 50,
cmap = custom_cmap(100, 'plasma_r', 0.05, 0.9, 'viridis', 0.2, 1),
subplot_kws = {'projection':ccrs.PlateCarree(central_longitude=180)},
transform =ccrs.PlateCarree(),
)
# add coast, set extent
for ax in p.axes.flat:
ax.add_feature(coast_full, facecolor='None', alpha=0.8, edgecolor='black', linewidth=1.5, zorder=5)
ax.set_extent((
hr_data_recon_all_tracks.lon.min(),hr_data_recon_all_tracks.lon.max(),
hr_data_recon_all_tracks.lat.min(),hr_data_recon_all_tracks.lat.max(),
))
# TODO codigos from_epsg independientes de site?
# SAVE one final track for riskscape
raster = hr_data_recon_all_tracks.max(dim='track').to_dataset().rename({'M':'z','lon':'x','lat':'y'}).copy()
raster = raster.rio.write_crs(CRS.from_epsg(4326)).copy()
raster.rio.to_raster(p_final_track_tif)
# Reproject the data using the crs object
raster_rep = raster.rio.reproject(CRS.from_epsg(32701))
raster_rep.rio.to_raster(p_final_track_tif_utm)
#raster_rep.rio.to_raster(p_final_track_tif_utm_risk)
SWATH#
Figure Explanation
Maximun wind fields.
p = hr_data_recon_all_tracks.max(dim='track').plot(
vmin = 30,
vmax = 60,
cmap = custom_cmap(100, 'plasma_r', 0.05, 0.9, 'viridis', 0.2, 1),
subplot_kws = {'projection': ccrs.PlateCarree(central_longitude=180)},
transform = ccrs.PlateCarree(),
)
p.axes.add_feature(coast_full, facecolor = 'None', alpha=0.8, edgecolor='black', linewidth=1.5, zorder=5)
plt.title('SWATH [m/s] -- Maximum Wind Fields')
Probabilistic assessment#
Figure Explanation
Winds higher than threshold probabilities.
plot_prob_assessment(
vortex_data_recon_all_tracks,
hr_data_recon_all_tracks,
wind_ths = [15,20,25,30],
coast = coast_full,
)
print("execution end: {0}".format(datetime.today().strftime('%Y-%m-%d %H:%M:%S')))