Damage

from datetime import datetime

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

execution start: 2022-09-08 22:16:50

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

# 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()
    
    # reforecast track
    n, tc_name = 166, 'amos' # AMOS 2016
    
    # riskscape
    districts_fn = 'savaii_district.shp'     
    site_main = 'Samoa'
    code_utm = 32702
    name_dist = 'District'
    
elif site == 'Upolu':

    # site
    site_ = site.lower()
    
    # reforecast track
    n, tc_name = 166, 'amos' # AMOS 2016
    
    # riskscape
    districts_fn = 'upolu_district.shp'     
    site_main = 'Samoa'
    code_utm = 32702
    name_dist = 'District'

elif site == 'Tongatapu':

    # site
    site_ = site
    
    # reforecast track
    n, tc_name = 180, 'gita'
    
    # riskscape
    districts_fn = 'tongatapu_divisions.shp'    
    site_main = 'Tongatapu'
    code_utm = 32701
    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')
event = 'Flooding_Metamodel_' + tc_name + '_tk_' + tk 
p_event = op.join(p_riskscape, site, 'data', 'rainfall', event + '.tif')

p_riskscape_site = op.join(p_riskscape, site)
p_riskscape_site_rain = op.join(p_riskscape, site, 'data', 'rainfall')

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

p_kepler_config_dmg

'/media/administrador/HD2/SamoaTonga/data/riskscape_projects/Savaii/config_files/config_damage_savaii.json'

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: es necesario? podemos darle otro nombre  al evento?
# event folder
#if tc_name == 'amos': aux = 'Amos'
#if tc_name == 'gita': aux = 'Gita'

#event = 'Rainfall_Flooding_{0}_Reforecast'.format(aux)

print('Running event: ' + p_event)

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

    # list = [p_riskscape_bin, ' model run ', model, ' --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))

Running event: /media/administrador/HD2/SamoaTonga/data/riskscape_projects/Savaii/data/rainfall/Flooding_Metamodel_amos_tk_max.tif
Build_rain_dam

People_rain_dam

Roads_rain_dam

Landcrops_rain_dam

--- 194.23740649223328 seconds ---

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

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) 

if site == 'Tongatapu':
    districts = districts.loc[
        (districts.index != 'Tonga') & (districts.index != 'Tongatapu') & \
        (districts.index!= "'Eua") & (districts.index != "Vahe 'Eua")]
    

Plot damage by district

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

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 amos Rainfall Flooding in Savaii: 11294607 2010 US dollars (USD).

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

execution end: 2022-09-08 22:20:30