Damage

from datetime import datetime

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

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

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

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.inset_locator import inset_axes


# raster tools
from rasterio.crs import CRS
import geopandas as gpd
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

# hide warnings
import warnings
warnings.filterwarnings('ignore')

Forecast Parameters

# project database
p_data = r'/media/administrador/HD2/SamoaTonga/data'
p_db = r'/media/administrador/HD1/DATABASES'

#site = 'Upolu'  # Savaii / Upolu / Tongatapu

# (optional) get site from config file
nb_args = read_config_args(op.abspath(''), '09d_damage')
site = nb_args['site']

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

Study site: Upolu

# Rainfall Flooding

if site != 'Savaii':
    models = ['Especials_rain_dam', 'Build_rain_dam', 'People_rain_dam', 'Roads_rain_dam', 'Landcrops_rain_dam']
else:
    models = ['Build_rain_dam', 'People_rain_dam', 'Roads_rain_dam', 'Landcrops_rain_dam']

# site related parameteres
if site == 'Savaii':

    # site
    site_ = site.lower()
    
    code_utm = 32702
    
    # riskscape
    districts_fn = 'savaii_district.shp'     
    site_main = 'Samoa'
    name_dist = 'District'

    
elif site == 'Upolu':

    # site
    site_ = site.lower()
    
    code_utm = 32702
    
    # riskscape
    districts_fn = 'upolu_district.shp'     
    site_main = 'Samoa'
    name_dist = 'District'


elif site == 'Tongatapu':

    # site
    site_ = site

    code_utm = 32701
    
    # riskscape
    districts_fn = 'tongatapu_divisions.shp'    
    site_main = 'Tongatapu'
    name_dist = 'name'

Database

tk = 'max' #Track number or max

p_site = op.join(p_data, site)

# deliverable folder
p_deliv = op.join(p_site, 'd09_rainfall_forecast')

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

p_riskscape_site = op.join(p_riskscape, site)
p_riskscape_data = op.join(p_riskscape, site, 'data')
p_riskscape_data_exp = op.join(p_riskscape, site_main, 'data')

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

p_kepler_config_dist = os.path.join(p_riskscape_site, 'config_files', 'config_damageDist_' + site.lower() + '.json')

# riskscape executable
p_bluemath = op.join(op.abspath(''), '..', '..', 'bluemath')
p_riskscape_bin = op.abspath(op.join(p_bluemath, 'riskscape','bin','riskscape'))

Forecast Output Folder

p_forecast = op.join(p_site, 'forecast', '09_rainfall_tc_inundation')

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

# choose last date
date = dates_exec[-1]

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

last available date: 202209010245

# 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: 202209010245

event = 'Flooding_Metamodel_' + tc_name + '_tk_' + tk + '_' + date + '_.tif'
p_event = op.join(p_riskscape, site, 'data', 'rainfall', event + '.tif')

Impact caused by TC

Buildings: Each building has associated a replacement cost expressed in 2010 US dollars (USD).

Roads and Special Infrastructures: None of these layers have a replacement cost associated as an attribute, therefore it has been obtained from the PCRAFI report (The World Bank 2013). The roads have also been classified based on the hazard value in transitable (water depth<0.3) or not transitable (water depth>0.3).

Crops: The typology has been divided into three main classes, tree crops, root crops, and intercrops, since the effects of the wind and flooding are different for each category.

Population: No economic damage is assessed to the population, however according to Cox et al. (2010) different safety threads have been defined based on water depth (similar with wind) what can be more useful on a forecasted event.

run riskscape model

os.chdir(p_riskscape_site)

# TODO como indicamos el input de riskscape? debe usar lo generado en 09a

start_time = time.time()
for model in models:
    
    path_output = op.join(p_riskscape_site_rain, event, model)
    list = [p_riskscape_bin, ' model run ', model, ' --parameter ' '"input-hazards.layer=' + p_event + '" ', ' --output ', path_output, ' -r']

    to_run = ''.join(list)
    stream = os.popen(to_run)
    print(model)
    output = stream.read()
    output
    
print("--- %s seconds ---" % (time.time() - start_time))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
/tmp/ipykernel_30748/58202454.py in <module>
      4 for model in models:
      5 
----> 6     path_output = op.join(p_riskscape_site_rain, event, model)
      7     list = [p_riskscape_bin, ' model run ', model, ' --parameter ' '"input-hazards.layer=' + p_event + '" ', ' --output ', path_output, ' -r']
      8 

NameError: name 'p_riskscape_site_rain' is not defined

Impact over assets - Output RiskScape

land_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Landcrops_rain_dam', 'event-impact.shp'))
build_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Build_rain_dam', 'event-impact.shp'))
if site != 'Savaii':
    spec_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Especials_rain_dam', 'event-impact.shp'))
people_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'People_rain_dam', 'event-impact.shp'))
roads_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Roads_rain_dam', 'event-impact.shp'))

Convert X, Y coordinates to Lon, Lat

Plot Impact

# Load Hazard - winds SWATH
a_file = open(p_kepler_config_dmg, "rb")
config_damage = json.load(a_file)
map_1 = KeplerGl(height=700, config=config_damage, show_docs=False)


if not build_dmg.empty:
    build_dmg = build_dmg.to_crs(epsg=4326)
    build_dmg = build_dmg[['Occ', 'Exposed', 'Loss', 'geometry']]
    map_1.add_data(data=build_dmg, name='Buildings Damage')

if not roads_dmg.empty:
    roads_dmg = roads_dmg.to_crs(epsg=4326)
    roads_dmg = roads_dmg[['CLASS', 'Impassab_1', 'Loss', 'geometry']]
    map_1.add_data(data=roads_dmg, name='Roads Damage')

if not people_dmg.empty:
    people_dmg = people_dmg.to_crs(epsg=4326)
    people_dmg = people_dmg[['Safety_Thr', 'Residents_', 'geometry']]
    map_1.add_data(data=people_dmg, name='People Damage')


if not land_dmg.empty:
    land_dmg = land_dmg.to_crs(epsg=4326)
    land_dmg = land_dmg[['CLASS_NAME', 'Loss', 'geometry']]
    map_1.add_data(data=land_dmg, name='Land Damage')


if site != 'Savaii':
    if not spec_dmg.empty:
        spec_dmg = spec_dmg.to_crs(epsg=4326)
        spec_dmg = spec_dmg[['TYPE', 'Loss', 'geometry']]
        map_1.add_data(data=spec_dmg, name='Specials Damage')

map_1

map_1

Damage by district

a_file = open(p_kepler_config_dist, "rb")
config_damage_dist = json.load(a_file)

Join damage by district

# p_riskscape_event = op.join(p_riskscape_site, event)

land_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event,  'Landcrops_rain_dam', 'event-impact.shp'))
build_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Build_rain_dam', 'event-impact.shp'))
if site != 'Savaii':
    spec_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Especials_rain_dam', 'event-impact.shp'))
people_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'People_rain_dam', 'event-impact.shp'))
roads_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Roads_rain_dam', 'event-impact.shp'))

districts = gpd.read_file(op.join(p_riskscape_site, 'data', districts_fn)).to_crs(epsg=4326)
districts = districts.set_index(name_dist)

if not build_dmg.empty:
    build_sim_dmg['Loss'] = build_sim_dmg['Loss'].replace(-9999.0, 0)
    build_sim_dmg = build_sim_dmg[['Occ', 'Exposed', 'Loss', 'geometry', name_dist]]
    build_sim_dmg = build_sim_dmg.to_crs(epsg=4326)
    districts = districts.join(build_sim_dmg.groupby(name_dist).sum()[['Loss']]\
                               .rename(columns={'Loss':'Loss Sim Buildings'}))

if not roads_dmg.empty:
    roads_sim_dmg['Loss'] = roads_sim_dmg['Loss'].replace(-9999.0, 0)
    roads_sim_dmg = roads_sim_dmg[['CLASS', 'Impassab_1', 'Loss', 'geometry',name_dist]]
    roads_sim_dmg = roads_sim_dmg.to_crs(epsg=4326)
    districts = districts.join(roads_sim_dmg.groupby(name_dist).sum()[['Loss']]\
                               .rename(columns={'Loss':'Loss Sim Roads'}))

if not people_dmg.empty:
    people_sim_dmg = people_sim_dmg[['Safety_Thr', 'Residents_', 'geometry',name_dist]]
    people_sim_dmg = people_sim_dmg.to_crs(epsg=4326)

if not land_sim_dmg.empty:
    land_sim_dmg['Loss'] = land_sim_dmg['Loss'].replace(-9999.0, 0)
    land_sim_dmg = land_sim_dmg[['CLASS_NAME', 'Loss', 'geometry',name_dist]]
    land_sim_dmg = land_sim_dmg.to_crs(epsg=4326)
    districts = districts.join(land_sim_dmg.groupby(name_dist).sum()[['Loss']]\
                               .rename(columns={'Loss':'Loss Sim Land'}))
    
if site != 'Savaii':
    if not spec_sim_dmg.empty:
        spec_sim_dmg['Loss'] = spec_sim_dmg['Loss'].replace(-9999.0, 0)
        spec_sim_dmg = spec_sim_dmg[['Loss', 'geometry']]#,name_dist]]
        spec_sim_dmg = spec_sim_dmg.to_crs(epsg=4326)
        # districts = districts.join(spec_sim_dmg.groupby(name_dist).sum()[['Loss']]\
        #                        .rename(columns={'Loss':'Loss Sim Specials'}))

districts = districts.fillna(0)
if site != 'Savaii':
    
    districts['TDamage_sim'] = districts['Loss Sim Buildings'] + districts['Loss Sim Roads']\
        + districts['Loss Sim Land'] #+ districts['Loss Sim Specials']
else:
    districts['TDamage_sim'] = districts['Loss Sim Buildings'] + districts['Loss Sim Roads']\
        + districts['Loss Sim Land'] # + districts['Loss Sim Specials']
    
districts['TDamage_sim'] = districts['TDamage_sim'].astype(int) 

Plot damage by district

map_1 = KeplerGl(height=700, data={"Damage": districts}, config=config_damage_dist, show_docs=False)
map_1

Damage components

land_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Landcrops_rain_dam', 'event-impact.shp')).to_crs(epsg=4326)
build_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Build_rain_dam', 'event-impact.shp')).to_crs(epsg=4326)
#spec_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Especials_rain_dam', 'event-impact.shp')).to_crs(epsg=4326)
people_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'People_rain_dam', 'event-impact.shp')).to_crs(epsg=4326)
roads_sim_dmg = gpd.read_file(op.join(p_riskscape_site_rain, event, 'Roads_rain_dam', 'event-impact.shp')).to_crs(epsg=4326)

build_sim_dmg['Loss'] = build_sim_dmg['Loss'].replace(-9999.0, 0)
#spec_sim_dmg['Loss'] = spec_sim_dmg['Loss'].replace(-9999.0, 0)
land_sim_dmg['Loss'] = land_sim_dmg['Loss'].replace(-9999.0, 0)
roads_sim_dmg['Loss'] = roads_sim_dmg['Loss'].replace(-9999.0, 0)

districts = gpd.read_file(op.join(p_riskscape_site, 'data', districts_fn)).to_crs(epsg=4326)
districts = districts.set_index(name_dist)

build_sim_dmg = build_sim_dmg[['Occ', 'Exposed', 'Loss', 'geometry', name_dist]]
roads_sim_dmg = roads_sim_dmg[['CLASS', 'Impassab_1', 'Loss', 'geometry',name_dist]]
land_sim_dmg = land_sim_dmg[['CLASS_NAME', 'Loss', 'geometry',name_dist]]
#spec_sim_dmg = spec_sim_dmg[['TYPE', 'Loss', 'geometry',name_dist]]
people_sim_dmg = people_sim_dmg[['Safety_Thr', 'Residents_', 'geometry',name_dist]]

districts = districts.join(build_sim_dmg.groupby(name_dist).sum()[['Loss']]\
    .rename(columns={'Loss':'Loss Sim Buildings'}))
districts = districts.join(land_sim_dmg.groupby(name_dist).sum()[['Loss']]\
    .rename(columns={'Loss':'Loss Sim Land'}))
#districts = districts.join(spec_sim_dmg.groupby(name_dist).sum()[['Loss']]\
#    .rename(columns={'Loss':'Loss Sim Specials'}))
districts = districts.join(roads_sim_dmg.groupby(name_dist).sum()[['Loss']]\
    .rename(columns={'Loss':'Loss Sim Roads'}))

districts = districts.fillna(0)



districts['TDamage_sim'] = districts['Loss Sim Buildings'] + districts['Loss Sim Roads']\
    + districts['Loss Sim Land'] # + districts['Loss Sim Specials']
districts['TDamage_sim'] = districts['TDamage_sim'].astype(int) 

centroids = districts.centroid
ranges = [0, 250000, 500000, 1000000, 1250000, 1500000, 1750000, 2000000]
cd = ['#FFFFFF', '#FFFFA6', '#FDFF09', '#EE8208', '#FB2908', '#680002', '#000000']

districts['color_sim'] = cd[0]
districts['color_sim'].where(districts['TDamage_sim'] < ranges[0], cd[0], inplace=True)
districts['color_sim'].where(districts['TDamage_sim'] < ranges[1], cd[1], inplace=True)
districts['color_sim'].where(districts['TDamage_sim'] < ranges[2], cd[2], inplace=True)
districts['color_sim'].where(districts['TDamage_sim'] < ranges[3], cd[3], inplace=True)
districts['color_sim'].where(districts['TDamage_sim'] < ranges[4], cd[4], inplace=True)
districts['color_sim'].where(districts['TDamage_sim'] < ranges[5], cd[5], inplace=True)
districts['color_sim'].where(districts['TDamage_sim'] < ranges[6], cd[6], inplace=True)

colors = ['#5052F8', '#1AC484', '#9943F9', '#FC4A80']

fig = plt.figure(figsize=(15,12))
ax_s = fig.add_subplot(1,1,1)

districts.plot(ax=ax_s, color=districts['color_sim'],  alpha=1, edgecolor='grey')

for dst in districts.index[:-1]:
    sizes_s = districts.loc[dst][
        #['Loss Sim Buildings','Loss Sim Roads', 'Loss Sim Land', 'Loss Sim Specials']].values
        ['Loss Sim Buildings', 'Loss Sim Land']].values
    lon, lat = districts.centroid[dst].x, districts.centroid[dst].y
    
    if not np.all((sizes_s == 0)):
        ax_sub_r = inset_axes(
            ax_s, width=0.6, height=0.6, loc=10, 
            bbox_to_anchor=(lon, lat), bbox_transform=ax_s.transData,
        )
        wedges, tests = ax_sub_r.pie(sizes_s, colors=colors, wedgeprops={'edgecolor':None})
     
ax_s.legend(wedges, ['Buildings', 'Roads', 'Land', 'Specials'])
ax_s.set_title('Damage by district', fontsize=15, fontweight='bold')

plt.show()

Aggregated damage

print('Total damage caused by TC ' + tc_name + ' Rainfall Flooding in '+ site +': {0} 2010 US dollars (USD).'\
    .format(np.nansum(districts['TDamage_sim'].values)))

Total damage caused by TC gita Rainfall Flooding in Tongatapu: 1571446 2010 US dollars (USD).

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