Risk Assessment

# TODO: Insert statistical framework for risk calculation

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
import matplotlib.colors as mcolors
from matplotlib.gridspec import GridSpec
import cmocean
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 calculate_risk

island = 'samoa'              
mode = 'hist'          
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)))

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

EAD, VaR 95 % & Worst Scenario

Calculate Risk Statistics

print('Rain Savaii & Upolu')
EAD_rain_sv, VaR95_rain_sv, Ws_rain_sv = calculate_risk(pr_savaii_rainTC)
EAD_rain_up, VaR95_rain_up, Ws_rain_up = calculate_risk(pr_upolu_rainTC)
print('\nWind Savaii & Upolu')
EAD_wind_sv, VaR95_wind_sv, Ws_wind_sv = calculate_risk(pr_savaii_windTC)
EAD_wind_up, VaR95_wind_up, Ws_wind_up = calculate_risk(pr_upolu_windTC)
print('\nCoast Savaii & Upolu')
EAD_coast_sv, VaR95_coast_sv, Ws_coast_sv = calculate_risk(pr_savaii_coastTC)
EAD_coast_up, VaR95_coast_up, Ws_coast_up = calculate_risk(pr_upolu_coastTC)

Rain Savaii & Upolu
 File 10 / 11
Wind Savaii & Upolu
 File 48 / 49
Coast Savaii & Upolu
 File 67 / 68

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')
EAD_coast = xr.concat([EAD_coast_sv, EAD_coast_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')
VaR95_coast = xr.concat([VaR95_coast_sv, VaR95_coast_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')
Ws_coast = xr.concat([Ws_coast_sv, Ws_coast_up], dim='x').sortby('x')

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)

Expected Annual Damage (TC-Wind, TC-Rain & TC-Coast)

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

Value at Risk 95% (TC-Wind, TC-Rain & TC-Coast)

VaR95_samoa = VaR95_rain + VaR95_wind + VaR95_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 = VaR95_samoa.VaR95.where(VaR95_samoa.VaR95 > 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()

Worst Scenario (TC-Wind, TC-Rain & TC-Coast)

Ws_samoa = Ws_rain + Ws_wind + Ws_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 = Ws_samoa.Ws.where(Ws_samoa.Ws > 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()

EAD at coarser spatial resolution (1km2)

ds_damage_1km = EAD_samoa.coarsen(x=10, y=10, boundary='pad').sum()
ds_damage_1km['EAD'] = ds_damage_1km.EAD / 1000

fig, ax0 = plt.subplots(1, constrained_layout=False, figsize=(25,16))
                 
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=None, zorder=5, alpha=0.9)
cb = fig.colorbar(im, ax=ax0, shrink=0.6, label='Economic Loss [2010 USD millions $]', format='%d')
ax0.set_title('Expected Annual Damage', fontweight='bold')

plt.show()