Hazard Winds

from datetime import datetime

print("execution start: {0}".format(datetime.today().strftime('%Y-%m-%d %H:%M:%S')))

execution start: 2022-09-13 08:10:02

import os
import os.path as op
import sys
import glob

import pandas as pd
import numpy as np
import xarray as xr

# raster tools
from rasterio.crs import CRS
import json
import rioxarray

# kepler 
from keplergl import KeplerGl

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

# operational utils
from operational.util import read_config_args

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.16.0

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(''), '10b_hazard')
site = nb_args['site']

print('Study site: {0}'.format(site))

Study site: Savaii

# site related parameteres
if site == 'Savaii':
    
    # reforecast track
    n, tc_name = 166, 'amos' # AMOS 2016
    tif_fn = 'hr_FINAL_track_AMOS_lonlat.tif'
    code_utm = 32702
    
elif site == 'Upolu':
    
    # reforecast track
    n, tc_name = 166, 'amos' # AMOS 2016
    tif_fn = 'hr_FINAL_track_AMOS_lonlat.tif'
    code_utm = 32702

elif site == 'Tongatapu':
    
    # reforecast track
    n, tc_name = 180, 'gita'
    tif_fn = 'hr_FINAL_track_GITA_lonlat.tif'
    code_utm = 32701

# site related parameteres
if site == 'Savaii': 
    code_utm = 32702    
elif site == 'Upolu': 
    code_utm = 32702
elif site == 'Tongatapu':    
    code_utm = 32702

Database

p_site = op.join(p_data, site)

# riskscape folder
p_riskscape = op.join(p_data, 'riskscape_projects')
p_riskscape_site = op.join(p_riskscape, site)

# kepler config files
p_kepler_config = op.join(p_riskscape_site, 'config_files', 'config_hazard_rain_{0}.json'.format(site.lower()))

Forecast Output Folder

p_forecast = op.join(p_site, 'forecast', '10_wind_damage')

# last valid execution of 10a
dates_exec = sorted([x for x in os.listdir(p_forecast) if len(glob.glob(op.join(p_forecast, x, 'Forecast_*.nc'))) > 0])
dates_exec = [x for x in dates_exec if 'reforecast' not in x]

if len(dates_exec) == 0:
    raise ValueError('No solved TCs available from 10a')

# choose last date
date = dates_exec[-1]

print('last available date: {0}'.format(date))

last available date: 202208310800

# available TC names
fs = glob.glob(op.join(p_forecast, date, 'Forecast_*.nc'))
names = list(set([op.basename(x).split('_')[1] for x in fs]))

print('available TCs: {0}'.format(names))

available TCs: ['amos']

# select tc name
tc_name = names[0]

# forecast folder
p_fore_date = op.join(p_forecast, date)
print('forecast date code: {0}'.format(date))

forecast date code: 202208310800

# tc forecast
p_fore_tc = op.join(p_fore_date, tc_name + '_forecast')

# wind hazard files
p_hr_tc = op.join(p_fore_tc, 'hr_{0}_tracks.nc'.format(tc_name))

if not op.isfile(p_hr_tc):
    raise ValueError('{0} not found'.format(p_hr_tc))
    

tk = 'max' #Track number or max

if tk == 'max':
    # file_name = op.join(p_out_tracks, 'Flooding_Metamodel_Analogues_envelope.nc')
    ds_wind = xr.open_dataset(p_hr_tc).max(dim = 'track')
else:
    ds_wind = xr.open_dataset(p_hr_tc).isel(track = tk)

ds_wind

<xarray.Dataset>
Dimensions:  (lat: 423, lon: 701)
Coordinates:
  * lat      (lat) float32 -13.84 -13.84 -13.84 -13.84 ... -13.42 -13.42 -13.42
  * lon      (lon) float32 -172.8 -172.8 -172.8 -172.8 ... -172.1 -172.1 -172.1
Data variables:
    M        (lat, lon) float64 nan nan nan nan nan ... 9.13 9.132 9.135 9.13

tif_riskscape = op.join(p_riskscape, site, 'data', 'wind', 'Flooding_Metamodel_' + tc_name + '_tk_' + tk + '_' + date + '.tif') #Always save in the same folder

print('Saving tif at: ' + tif_riskscape)

Saving tif at: /media/administrador/HD2/SamoaTonga/data/riskscape_projects/Savaii/data/wind/Flooding_Metamodel_amos_tk_max_202208310800.tif

Hazard Winds

# TODO codigos from_epsg independientes de site?

raster = ds_wind.rename({'M':'z','lon':'x','lat':'y'}).rio.write_crs(CRS.from_epsg(4326)).copy()
raster_rep = raster.rio.reproject(CRS.from_epsg(code_utm))
raster.rio.to_raster(tif_riskscape)

data_sim_84 = raster.copy()
data_sim_84 = data_sim_84.where((data_sim_84>=0.2) & (data_sim_84<=100))

df_sim_84 = data_sim_84.to_dataframe().reset_index().dropna()\
    .rename(columns={'x':'Lon', 'y':'Lat'})

Plot Hazard - SWATH winds

a_file = open(p_kepler_config, "rb")
config_winds = json.load(a_file)

# Hazard - flooding depth
map_1 = KeplerGl(height=700, data={"Winds - SWATH": df_sim_84}, config=config_winds, show_docs=False)
map_1

print("execution end: {0}".format(datetime.today().strftime('%Y-%m-%d %H:%M:%S')))

execution end: 2022-09-13 08:10:05