Damage

# TODO:
# - donde coge el input para riskscape generado en 10a?

from datetime import datetime

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

execution start: 2022-09-13 08:30: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 = 'Savaii'  # 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()

    # riskscape
    districts_fn = 'savaii_district.shp'     
    site_main = 'Samoa'
    name_dist = 'District'

    
elif site == 'Upolu':

    # site
    site_ = site.lower()

    # riskscape
    districts_fn = 'upolu_district.shp'     
    site_main = 'Samoa'

    name_dist = 'District'

elif site == 'Tongatapu':

    # site
    site_ = site
    
    # riskscape
    districts_fn = 'tongatapu_divisions.shp'    
    site_main = 'Tongatapu'
    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_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'))

Forecast Output Folder

p_forecast = op.join(p_site, 'forecast', '10_wind_damage')

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

# choose last date
date = dates_exec[-1]

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

last available date: 202208310800

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

tk = 'max'
event = 'Flooding_Metamodel_' + tc_name + '_tk_' + tk + '_' + date 
p_event = op.join(p_riskscape, site, 'data', 'wind', 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)

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

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

NameError: name 'p_riskscape_site_wind' is not defined

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

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