Coastal Flooding Impact#

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

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

# pip
import xarray as xr
import netCDF4 as nc
import matplotlib.pyplot as plt
import numpy as np
import numpy.ma as ma
import pandas as pd
import pickle as pkl

import geopandas as gpd
import json

from keplergl import KeplerGl

import rioxarray
from pyproj import CRS

# 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

Forecast Parameteres#

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

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

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

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

Study site: Savaii

# hazard type
hazard_type = 'TCs'

# Flooding models
models = ['Landcrops_coastal_dam', 'Build_coastal_dam', 'Especials_coastal_dam', 'People_coastal_dam', 'Roads_coastal_dam']

# event name
event = '06d_flooding'

if site == 'Tongatapu':

    site_g = 'Tongatapu'    
    code_utm = 32701

elif site == 'Upolu':
        
    site_g = 'Samoa'
    code_utm = 32702
        
elif site == 'Savaii':
    
    site_g = 'Samoa'   
    code_utm = 32702

Database#

# database
p_site = op.join(p_data, site_g)

# 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_kepler_config_dmg = op.join(p_riskscape, site, 'config_files', 'config_damage_{0}.json'.format(site.lower()))

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

Forecast Output Folder#

# forecast folder
p_forecast = op.join(p_site, 'forecast', '06_tc_inundation_forecast')

# last valid execution of 05a_swell_inundation_forecast_waves
dates_exec = sorted([x for x in os.listdir(p_forecast) if op.isfile(op.join(p_forecast, x, 'flooding_{0}'.format(site), 'Flooding_Maps.nc'))])

if len(dates_exec) == 0:
    raise ValueError('No flooding input available from 06c')

# choose last date
date = dates_exec[-1]

print('forecast execution: {0}'.format(date))

forecast execution: reforecast_evan

  • Select TC Mode, and Name

tc_mode = 'synthetic' # 'historical'/'synthetic'

# TC id is choosen at notebook 06b
name = 'tc'

# forecast date code
p_fore_date = op.join(p_forecast, date)
print('forecast date code: {0}'.format(date))

# prepare forecast subfolders for this date
p_output = op.join(p_fore_date, 'flooding_{0}'.format(site))  

# input file from 06c
p_flooding_maps = op.join(p_output, 'Flooding_Maps.nc')

# output event
event = 'Flooding_Max_{0}_{1}.tif'.format(hazard_type, name)
p_event = op.join(p_riskscape_data, 'coastal', event)

forecast date code: reforecast_evan

Flooding to raster#

# load flooding maps data
depth_max = xr.open_dataset(p_flooding_maps).max(dim='time').to_dataframe().reset_index()

# Resample x y coordinates
X, Y = np.meshgrid(
    np.arange(depth_max.x.min(), depth_max.x.max(), 5), 
    np.arange(depth_max.y.min(), depth_max.y.max(), 5),
)

df_aux = pd.DataFrame(
    {
        'x': np.reshape(X, -1),
        'y': np.reshape(Y, -1),
        'z': np.zeros((len(np.reshape(X, -1)))),
    }).set_index(['x', 'y'])

ds_rain = pd.concat([depth_max.set_index(['x','y']), df_aux]).reset_index().groupby(['x', 'y']).max().to_xarray()

raster = ds_rain.rio.write_crs('EPSG:{0}'.format(code_utm)).transpose('y', 'x')
    
raster = raster.where(raster < 5)

raster.rio.to_raster(p_event)

RiskScape#

os.chdir(p_riskscape_site)

start_time = time.time()

for model in models:

    path_output = op.join(p_riskscape_site, 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_coastal_dam

Build_coastal_dam

Especials_coastal_dam

People_coastal_dam

Roads_coastal_dam

--- 466.9858179092407 seconds ---

Output#

land_dmg = gpd.read_file(op.join(p_riskscape_site, event, 'Landcrops_coastal_dam', 'event-impact.shp'))
build_dmg = gpd.read_file(op.join(p_riskscape_site, event, 'Build_coastal_dam', 'event-impact.shp'))
spec_dmg = gpd.read_file(op.join(p_riskscape_site, event, 'Especials_coastal_dam', 'event-impact.shp'))
people_dmg = gpd.read_file(op.join(p_riskscape_site, event, 'People_coastal_dam', 'event-impact.shp'))
roads_dmg = gpd.read_file(op.join(p_riskscape_site, event, 'Roads_coastal_dam', 'event-impact.shp'))

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