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-10 17:22:27

# 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 = 'Savaii'  # 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'    
    nm_reforecast = 'harold'
    code_utm = 32701

elif site == 'Upolu':
        
    site_g = 'Samoa'
    nm_reforecast = 'evan'
    code_utm = 32702
        
elif site == 'Savaii':
    
    site_g = 'Samoa'   
    nm_reforecast = 'evan'
    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()))
p_kepler_config_hazard = op.join(p_riskscape_site, 'config_files', 'config_hazard_flooding_{0}.json'.format(site))

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

Reforecast

date = 'reforecast_{0}'.format(nm_reforecast) 

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

forecast execution: reforecast_evan

• Select TC Mode, and Name

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

name = nm_reforecast.upper()

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

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
p_event = op.join(p_riskscape_data, 'coastal', 'Flooding_Max_{0}_{1}.tif'.format(hazard_type, name))

forecast date code: reforecast_evan

Flooding to raster

# load flooding maps data
depth_max = xr.open_dataset(p_flooding_maps).max(dim='time').rename({'depth':'z'}).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)

data_sim_84 = raster.rio.reproject(CRS.from_epsg(4326)).copy()
data_sim_84 = data_sim_84.where((data_sim_84>=0.2) & (data_sim_84<=100))
df_sim_84 = data_sim_84.to_dataframe().reset_index().dropna().rename(columns={'x':'Lon', 'y':'Lat'})#, 'depth':'z'})

Plot Hazard - SWATH

a_file = open(p_kepler_config_hazard, "rb")
config_hazard= json.load(a_file)

# Hazard - flooding depth
map_1 = KeplerGl(height=700, data={"Depth ": df_sim_84}, config=config_hazard, show_docs=False)

map_1

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

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