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-03 00:57:52

# 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(''), '05e_flooding_impact')
site = nb_args['site']

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

Study site: Savaii

# hazard type
hazard_type = 'Swells'

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

# event name
event = '05e_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', '05_swell_inundation_system')

# 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, 'nearshore', 'flooding_{0}'.format(site), 'Flooding_Maps.nc'))])

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

# choose last date
date = dates_exec[-1]

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

forecast execution: 20220831

# 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_fore_near = op.join(p_fore_date, 'nearshore')
p_output = op.join(p_fore_near, 'flooding_{0}'.format(site))  

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

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

forecast date code: 20220831

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

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