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-11 17:26:36

# 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'
p_lf = op.join(p_db, 'Lisflood_Coastal')

# site = 'Savaii'  # 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: Upolu

# 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
    ev_date = ['2018-05-26', '2018-05-28']

elif site == 'Upolu':
        
    site_g = 'Samoa'
    code_utm = 32702
    ev_date = ['2015-01-22', '2015-01-24']
        
elif site == 'Savaii':
    
    site_g = 'Samoa'   
    code_utm = 32702
    ev_date = ['2015-01-22', '2015-01-24']
    

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_hazard = op.join(p_riskscape_site, 'config_files', 'config_hazard_rain_{0}.json'.format(site.lower()))
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'))

Forecast Output Folder

# forecast folder
p_forecast = op.join(p_site, 'forecast', '05_swell_inundation_system')

date = 'reforecast_{0}_{1}'.format(ev_date[0], ev_date[1])

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

date: reforecast_2015-01-22_2015-01-24

# 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: reforecast_2015-01-22_2015-01-24

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),
        'depth': np.zeros((len(np.reshape(X, -1)))),
    }).set_index(['x', 'y'])

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

raster = ds.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'})
df_sim_84

	Lon	Lat	z	spatial_ref
113067	-172.073475	-13.866101	1.918868	0
129784	-172.073337	-13.866009	0.508727	0
135356	-172.073291	-13.865962	0.592278	0
135357	-172.073291	-13.866009	0.570926	0
135358	-172.073291	-13.866055	0.396399	0
...	...	...	...	...
78601066	-171.423755	-14.024806	1.276123	0
78601067	-171.423755	-14.024852	1.307839	0
78601068	-171.423755	-14.024898	1.538250	0
78601069	-171.423755	-14.024944	1.005600	0
78601070	-171.423755	-14.024990	1.001868	0

347271 rows × 4 columns

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

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