Damage
Contents
Damage#
from datetime import datetime
print("execution start: {0}".format(datetime.today().strftime('%Y-%m-%d %H:%M:%S')))
execution start: 2022-09-12 20:11:47
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
# 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))
# 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)
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))
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-12 20:12:58