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-02 18:12:36

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: Savaii

# current date (for execution folder)
#date = datetime.today().strftime('%Y%m%d%H%M')
#print('forecast execution date: {0}'.format(date))

# 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',
]

# TODO reforecast
ags_sel = [' AVNI', ' NVGI', ' JTWC', ' NVGM',' JAVI' ]

# site related parameteres
if site == 'Savaii':
    site_center = (-172.40+360,-13.65)
    
    # area of interest
    area_interest = [[177, -165+360] , [-21.5, -7]]
    
    # reforecast track
    n, tc_name = 166, 'amos' # AMOS 2016
    name_f_file = 'ash202016_amos_aids.dat'

elif site == 'Upolu':
    site_center = (-171.78+360,-13.92)
    
    # area of interest
    area_interest = [[177, -165+360] , [-21.5, -7]]
    
    # reforecast track
    n, tc_name = 166, 'amos' # AMOS 2016
    name_f_file = 'ash202016_amos_aids.dat'

elif site == 'Tongatapu':
    site_center = (-175.19+360,-21.16)
    
    # area of interest
    area_interest = [[179, -169+360] , [-23.6, -18.1]]
    
    # reforecast track
    n, tc_name = 180, 'gita'
    name_f_file = 'ash092018_gita_aids.dat'
    
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')

# TODO reforecast tracks
path_tc_rf = os.path.join(p_deliv, name_f_file)

# load the high-resolution coastline for the plots
coast_full = ShapelyFeature(
    Reader(p_coast).geometries(),
    ccrs.PlateCarree(), edgecolor='black',
)

## Reforecast

date = 'reforecast_{0}'.format(tc_name) 

print('forecast execution: {0}'.format(date))

forecast execution: reforecast_amos

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: reforecast_amos

Download TC warnings

# download TC warnings
#files_ash, files_bsh = download_warning(p_fore_warn, code_basin='sh')

# TODO REFORECAST
files_ash = [path_tc_rf]
files_bsh = ['test']

# stop forecast if no warnings
if len(files_bsh) == 0:
    raise ValueError('No warnings in the SH. Forecast execution cancelled.')
    sys.exit()
    

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

active TC: ash202016_amos_aids.dat
best track: test

# 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]))

name warning: 20
last available warning time: 2016-04-25T12:00:00.000000000

# use last available time warning 
time_warning = str(available_times_warning[-1])
print('time warning: {0}'.format(time_warning))

time warning: 2016-04-25T12:00:00.000000000

# TODO reforecast
name_warning = tc_name

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!!

# TODO reforecast
#p_track = os.path.join(p_deliv,'sp_{0}.nc'.format(n))

#tco = xr.open_dataset(p_track)
#tc = process_tc_file(tco)
#tc = tc.isel(date_time=slice(150,385))

#Plot_track(tc, figsize=[20,10], cmap='jet_r', m_size=10);

Agencies Selection

print('Number of tracks: {0}'.format(len(TRACKS.track.values)))
print('Available Agencies:\n{0}\n------------'.format(TRACKS.Agency.values))

Number of tracks: 44
Available Agencies:
[' CARQ' ' AHNI' ' AVNI' ' CLIP' ' CLSW' ' CMES' ' CMSD' ' CONW' ' COT2'
 ' CTC2' ' DSHA' ' ECS5' ' EEM2' ' GFN2' ' JAVI' ' JJG2' ' JNVI' ' LGEA'
 ' LGEN' ' N012' ' N032' ' N042' ' N052' ' N062' ' N072' ' N082' ' N092'
 ' N102' ' N112' ' N122' ' N132' ' N142' ' N152' ' N162' ' N172' ' N202'
 ' NVGI' ' NVS1' ' NVS5' ' RVCN' ' S5XX' ' S5YY' ' SHIA' ' XTRP']
------------

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)

Total number of tracks: 3
Agencies from list available : [' AVNI' ' JAVI' ' NVGI']

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)

(<Figure size 2160x1080 with 4 Axes>,
 <cartopy.mpl.geoaxes.GeoAxesSubplot at 0x7f7f7c5afd90>,
 [177, 195, -21.5, -7])

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!!')
        

 Saving class attributes... 


 All data will be saved in /media/administrador/HD1/DATABASES/WRF_Wind    to study the Savaii domain,    which center location is (187.6, -13.65)    and forecast cyclon TRACK 0 will be analyzed!! 


 Calculating/saving vortex model... 

Done!!

 Plotting vortex model figures... 

Done!!

 Calculating/saving ERA5 data... 

Done!!

 Calculating/saving high resolution winds data... 

Done!!

 Post-processing the ERA5 and HR datasets... 

Done!!

 Getting the lons/lats to reconstruct the HR winds 

Done!!
TRACK 0 analyzed!!

 Saving class attributes... 


 All data will be saved in /media/administrador/HD1/DATABASES/WRF_Wind    to study the Savaii domain,    which center location is (187.6, -13.65)    and forecast cyclon TRACK 1 will be analyzed!! 


 Calculating/saving vortex model... 

Done!!

 Plotting vortex model figures... 

Done!!

 Calculating/saving ERA5 data... 

Done!!

 Calculating/saving high resolution winds data... 

Done!!

 Post-processing the ERA5 and HR datasets... 

Done!!

 Getting the lons/lats to reconstruct the HR winds 

Done!!
TRACK 1 analyzed!!

 Saving class attributes... 


 All data will be saved in /media/administrador/HD1/DATABASES/WRF_Wind    to study the Savaii domain,    which center location is (187.6, -13.65)    and forecast cyclon TRACK 2 will be analyzed!! 


 Calculating/saving vortex model... 

Done!!

 Plotting vortex model figures... 

Done!!

 Calculating/saving ERA5 data... 

Done!!

 Calculating/saving high resolution winds data... 

Done!!

 Post-processing the ERA5 and HR datasets... 

Done!!

 Getting the lons/lats to reconstruct the HR winds 

Done!!
TRACK 2 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')

Text(0.5, 1.0, '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')))

execution end: 2022-09-02 18:17:30