Risk Assessment
Contents
Risk Assessment#
import os
import os.path as op
import sys
import glob
import natsort
import xarray as xr
import matplotlib.pyplot as plt
from matplotlib import colors
from matplotlib.colors import ListedColormap
from matplotlib.gridspec import GridSpec
import cmocean
import matplotlib.colors as mcolors
import matplotlib
import numpy as np
import numpy.ma as ma
import pandas as pd
from natsort import natsorted
import contextily as ctx
# geospatial
import rioxarray
sys.path.insert(0, op.join(os.getcwd(), '..', 'lib'))
from plot import Map
from config import *
from data_process import agg_damage_syn, calculate_risk_syn
island = 'samoa'
mode = 'syn'
code_utm = code_utm_samoa
# main project folder
p_main = op.join(os.getcwd(), '..')
p_data = op.join(p_main ,'data')
p_output = op.join(p_main, 'output')
p_tcs = op.join(p_main, 'data', 'TCs_')
# load bathymetry and dem
bathy = rioxarray.open_rasterio(op.join(p_data, '{0}_smooth_utm.tif'.format(island))).drop('band').squeeze()
dem_savaii = xr.open_dataset(op.join(p_data, 'savaii_5m.nc')).isel(x=slice(None, None, 20), y=slice(None, None, 20))
dem_upolu = xr.open_dataset(op.join(p_data, 'upolu_5m.nc')).isel(x=slice(None, None, 20), y=slice(None, None, 20))
tcs = xr.open_dataset(op.join(p_tcs, 'TCs_r8_{0}_{1}.nc'.format(mode, island)))
tcs_ids = xr.open_dataset(op.join(p_tcs, 'Tesla_sim_newTCIds_Samoa.nc'))
# select only first simulation
#tcs = tcs.sel
# load dem
dem_savaii = xr.open_dataset(op.join(p_data, 'savaii_5m.nc')).isel(x=slice(None, None, 20), y=slice(None, None, 20))
dem_upolu = xr.open_dataset(op.join(p_data, 'upolu_5m.nc')).isel(x=slice(None, None, 20), y=slice(None, None, 20))
dem_samoa = xr.merge([dem_savaii.sel(x=slice(None, dem_upolu.x.min()-10)), dem_upolu])
dem_samoa = dem_samoa.where(dem_samoa > 0)
mapRk = Map()
mapRk.dem = dem_samoa
# load damages
pr_savaii_rainTC = glob.glob(op.join(p_main, 'riskscape', 'io', 'io_savaii', 'results', mode, 'rainTC*', 'damage_100_TC*.nc'))
pr_savaii_windTC = glob.glob(op.join(p_main, 'riskscape', 'io', 'io_savaii', 'results', mode, 'windTC*', 'damage_100_TC*.nc'))
pr_savaii_coastTC = glob.glob(op.join(p_main, 'riskscape', 'io', 'io_savaii', 'results', mode, 'coastTC*', 'damage_100_TC*.nc'))
pr_upolu_rainTC = glob.glob(op.join(p_main, 'riskscape', 'io', 'io_upolu', 'results', mode, 'rainTC*', 'damage_100_TC*.nc'))
pr_upolu_windTC = glob.glob(op.join(p_main, 'riskscape', 'io', 'io_upolu', 'results', mode, 'windTC*', 'damage_100_TC*.nc'))
pr_upolu_coastTC = glob.glob(op.join(p_main, 'riskscape', 'io', 'io_upolu', 'results', mode, 'coastTC*', 'damage_100_TC*.nc'))
ids = tcs_ids.sel(n_sim=0).storm.to_dataframe().dropna().storm.values.astype(int)
pr_savaii_rainTC_new = []
for file in pr_savaii_rainTC:
if int(file.split('damage_100_TC')[1].split('.nc')[0]) in ids:
pr_savaii_rainTC_new.append(file)
pr_upolu_rainTC_new = []
for file in pr_upolu_rainTC:
if int(file.split('damage_100_TC')[1].split('.nc')[0]) in ids:
pr_upolu_rainTC_new.append(file)
EAD, VaR 95 % & Worst Scenario#
Calculate Risk Statistics
n_years = 80
print('Rain Savaii & Upolu')
EAD_rain_sv, VaR95_rain_sv, Ws_rain_sv = calculate_risk_syn(pr_savaii_rainTC_new, n_years)
EAD_rain_up, VaR95_rain_up, Ws_rain_up = calculate_risk_syn(pr_upolu_rainTC_new, n_years)
print('\nWind Savaii & Upolu')
EAD_wind_sv, VaR95_wind_sv, Ws_wind_sv = calculate_risk_syn(pr_savaii_windTC, n_years)
EAD_wind_up, VaR95_wind_up, Ws_wind_up = calculate_risk_syn(pr_upolu_windTC, n_years)
print('\nCoast Savaii & Upolu')
#EAD_coast_sv, VaR95_coast_sv, Ws_coast_sv = calculate_risk(pr_savaii_coastTC, tcs)
#EAD_coast_up, VaR95_coast_up, Ws_coast_up = calculate_risk(pr_upolu_coastTC, tcs)
Rain Savaii & Upolu
File 40 / 41
Wind Savaii & Upolu
File 168 / 169
Coast Savaii & Upolu
EAD_rain = xr.concat([EAD_rain_sv, EAD_rain_up], dim='x').sortby('x')
EAD_wind = xr.concat([EAD_wind_sv, EAD_wind_up], dim='x').sortby('x')
VaR95_rain = xr.concat([VaR95_rain_sv, VaR95_rain_up], dim='x').sortby('x')
VaR95_wind = xr.concat([VaR95_wind_sv, VaR95_wind_up], dim='x').sortby('x')
Ws_rain = xr.concat([Ws_rain_sv, Ws_rain_up], dim='x').sortby('x')
Ws_wind = xr.concat([Ws_wind_sv, Ws_wind_up], dim='x').sortby('x')
#EAD_coast = xr.concat([EAD_coast_sv, EAD_coast_up], dim='x').sortby('x')
EAD_coast, VaR95_coast, Ws_coast = [], [], []
EAD_coast, VaR95_coast, Ws_coast = xr.Dataset(), xr.Dataset(), xr.Dataset()
Plot
mapRk.gridDamageResults(EAD_rain, EAD_wind, EAD_coast, VaR95_rain, VaR95_wind, VaR95_coast, Ws_rain, Ws_wind, Ws_coast, samoa_extent, code_utm, plot_coast=False)
TC-Wind, TC-Rain & TC-Coast damage integration#
EAD_samoa = EAD_rain + EAD_wind #+ EAD_coast
fig = plt.figure(constrained_layout=False, figsize=(25,8))
gs = GridSpec(1, 1, figure=fig)
ax0 = fig.add_subplot(gs[0,0])
ax0 = mapRk.BackGround_1cmaps(fig, ax0, samoa_extent, cmocean.cm.turbid, code_utm, cbar=False, vmin=-2000, vmax=1000, alpha=0.9)
im = EAD_samoa.EAD.where(EAD_samoa.EAD > 0).plot(ax=ax0, x='x', cmap='hot_r', add_colorbar=False, norm=mcolors.LogNorm(), vmin=1, vmax=170000, zorder=5)
cb = fig.colorbar(im, ax=ax0, shrink=0.99, label='Economic Loss [2010 USD $]', format='%d')
ax0.set_title('Aggregated EAD $(TC-Rain + TC-Wind + TC-Coast)$', fontweight='bold')
plt.show()
Qualitative Damage#
ds_damage_1km = EAD_samoa.coarsen(x=10, y=10, boundary='pad').sum()
ds_damage_1km['EAD'] = ds_damage_1km.EAD / 1000
col_dict = {1:"#159A36", 2:"#79C533", 3:"#FFFF0C", 4:"#FEAB0B", 5:"#F80007"}
cm = ListedColormap([col_dict[x] for x in col_dict.keys()])
labels = np.array( ['Minimum', 'Low', 'Moderate', 'High', 'Extreme'])
len_lab = len(labels)
bounds = [0, 10, 50, 500, 1000, 2000]
norm = matplotlib.colors.BoundaryNorm(bounds, len_lab, clip=True)
fmt = matplotlib.ticker.FuncFormatter(lambda x, pos: labels[norm(x)])
fig = plt.figure(constrained_layout=False, figsize=(25,16))
gs = GridSpec(2, 1, figure=fig)
ax0 = fig.add_subplot(gs[0])
ax0 = mapRk.BackGround_1cmaps(fig, ax0, samoa_extent, cmocean.cm.turbid, code_utm, cbar=False, vmin=-2000, vmax=1000, alpha=0.9, source=ctx.providers.OpenStreetMap.Mapnik)
im = ds_damage_1km.EAD.where(ds_damage_1km.EAD > 0).plot(ax=ax0, x='x', cmap='hot_r', add_colorbar=False, norm=mcolors.LogNorm(), vmin=1, vmax=2000, zorder=5, alpha=0.9)
cb = fig.colorbar(im, ax=ax0, shrink=0.9, label='Economic Loss [2010 USD millions $]', format='%d')
ax0.set_title('Expected Annual Damage', fontweight='bold')
ax1 = fig.add_subplot(gs[1])
ax1 = mapRk.BackGround_1cmaps(fig, ax1, samoa_extent, cmocean.cm.turbid, code_utm, cbar=False, vmin=-2000, vmax=1000, alpha=0.9, source=ctx.providers.OpenStreetMap.Mapnik)
im = ds_damage_1km.EAD.where(ds_damage_1km.EAD > 0).plot(x='x', cmap=cm, add_colorbar=False, norm=norm, alpha=0.7, zorder=5)
tickz = np.diff(bounds) / 2 + bounds[:-1]
cb = fig.colorbar(im, ax=ax1, shrink=0.9, format=fmt, ticks=tickz, label='Coastal Risk Level')
ax1.set_title('')
plt.show()
