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:45:30

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(''), '10d_damage')
site = nb_args['site']

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

Study site: Savaii

# Rainfall Flooding
if site != 'Savaii':
    models = ['Landcrops_wind_dam', 'Build_wind_dam', 'Especials_wind_dam', 'People_wind_dam']
else:
    models = ['Landcrops_wind_dam', 'Build_wind_dam', 'People_wind_dam'] 

# event name
# event = '10d_flooding'

# 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'
    event = 'hr_FINAL_track_AMOS_utm_risk'
    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'
    event = 'hr_FINAL_track_AMOS_utm_risk'
    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'
    event = 'hr_FINAL_track_GITA_utm_risk'
    name_dist = 'name'

Database

p_site = op.join(p_data, site)

# riskscape folder
p_riskscape = op.join(p_data, 'riskscape_projects')
p_riskscape_site = op.join(p_riskscape, site)
p_riskscape_site_wind = op.join(p_riskscape, site, 'data', 'wind')

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

## Reforecast

tk = 'max'
event = 'Flooding_Metamodel_' + tc_name + '_tk_' + tk 
p_event = op.join(p_riskscape, site, 'data', 'wind', event + '.tif')

Forecast Output Folder

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)

#a = rioxarray.open_rasterio(op.join('/media/administrador/HD2/SamoaTonga/data/riskscape_projects/Upolu/data/wind','hr_FINAL_track_amos_utm_risk.tif'))

start_time = time.time()

for model in models:

    path_output = op.join(p_riskscape_site_wind, 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))

Landcrops_wind_dam

Build_wind_dam

WARNING: No output has been produced from pipeline step(s) 'sample_area_layer:[select], sampled:[select], select_5:[select], raw_results:[select], event_impact_table:[select], report_event-impact:[select], filter:[filter]'. Possible causes may be: the pipeline uses a filter condition that does not match any data; the model uses a join condition based on geometry, but the datasets do not overlap geographically; or for CSV-based data, the incorrect CRS (or `crs-longitude-first`) bookmark setting is being used.

People_wind_dam

--- 126.05045056343079 seconds ---

Impact over assets - Output RiskScape

land_dmg = gpd.read_file(op.join(p_riskscape_site_wind, event, 'Landcrops_wind_dam', 'event-impact.shp'))
build_dmg = gpd.read_file(op.join(p_riskscape_site_wind, event, 'Build_wind_dam', 'event-impact.shp'))
people_dmg = gpd.read_file(op.join(p_riskscape_site_wind, event, 'People_wind_dam', 'event-impact.shp'))
if site != 'Savaii':
    spec_dmg = gpd.read_file(op.join(p_riskscape_site_wind, event, 'Especials_wind_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 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_wind, event, 'Landcrops_wind_dam', 'event-impact.shp'))
build_sim_dmg = gpd.read_file(op.join(p_riskscape_site_wind, event, 'Build_wind_dam', 'event-impact.shp'))
if site != 'Savaii':
    spec_sim_dmg = gpd.read_file(op.join(p_riskscape_site_wind, event, 'Especials_wind_dam', 'event-impact.shp'))
    spec_sim_dmg['Loss'] = spec_sim_dmg['Loss'].replace(-9999.0, 0)
people_sim_dmg = gpd.read_file(op.join(p_riskscape_site_wind, event, 'People_wind_dam', 'event-impact.shp'))
                                                
build_sim_dmg['Loss'] = build_sim_dmg['Loss'].replace(-9999.0, 0)
land_sim_dmg['Loss'] = land_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)

if not build_sim_dmg.empty:

    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 land_sim_dmg.empty:
    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 = spec_sim_dmg[['TYPE', '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'}))

if not people_sim_dmg.empty:
    people_sim_dmg = people_sim_dmg[['Safety_Thr', 'Residents_', 'geometry', name_dist]]
    people_sim_dmg = people_sim_dmg.to_crs(epsg=4326)
    
A=[]
for x in range(len(districts.columns.values)):
    if 'Loss' in districts.columns.values[x]:
        A.append(x)
districts['TDamage_sim'] = np.sum(districts.iloc[:,A], axis=1).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

with open(p_kepler_config_dist, “w”) as outfile: json.dump(map_1.config, outfile)

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: 555911 2010 US dollars (USD).

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

execution end: 2022-09-08 22:47:58