Coastal Flooding Impact
Contents
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-18 01:09:41
# 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 = 'Tongatapu' # 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: Tongatapu
# 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_2018-05-26_2018-05-28
# 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_2018-05-26_2018-05-28
Flooding to raster#
# load flooding maps data
ds = xr.open_dataset(p_flooding_maps).max(dim='time').to_dataframe().reset_index()
ds['y'] = np.round(ds['y'],4)
ds = ds.set_index(['x', 'y']).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 | |
|---|---|---|---|---|
| 2671246 | -175.329523 | -21.107935 | 0.752804 | 0 |
| 2671258 | -175.329523 | -21.108501 | 1.027955 | 0 |
| 2671259 | -175.329523 | -21.108548 | 1.032495 | 0 |
| 2671260 | -175.329523 | -21.108595 | 1.029771 | 0 |
| 2671261 | -175.329523 | -21.108642 | 1.045208 | 0 |
| ... | ... | ... | ... | ... |
| 15353199 | -175.199195 | -21.141932 | 0.301485 | 0 |
| 15353200 | -175.199195 | -21.141979 | 0.295128 | 0 |
| 15353201 | -175.199195 | -21.142026 | 0.281507 | 0 |
| 15353202 | -175.199195 | -21.142073 | 0.287864 | 0 |
| 15353203 | -175.199195 | -21.142120 | 0.273334 | 0 |
503974 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
--- 118.3973298072815 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