Vaisigano Flooding Event

import os
import os.path as op
import sys
import warnings
warnings.filterwarnings('ignore')
from IPython.display import Image
import numpy as np
import earthpy.spatial as es
import cmocean

sys.path.insert(0, op.join(op.abspath(''), '..','..',  '..'))

from bluemath.rainfall_forecast.rainfall_forecast import *

Warning: cannot import cf-json, install "metocean" dependencies for full functionality

# database
p_data = r'/media/administrador/HD2/SamoaTonga/data'
site, nm = 'upolu', 'up'
p_site = op.join(p_data, site)

p_deliv = op.join(p_site, 'd09_rainfall_forecast')

#Basins path
p_basins = op.join(p_deliv, 'basins', nm + '_db' +  '_fixed.shp')

#Dem
p_dem = os.path.join(p_deliv, site + '_5m.nc')

#Vaisigano case
p_apia2001 = op.join(p_deliv, 'apia2001')
n_case = 'up20_out_2001'
n_case1 = 'up1_out_2001'
d_o1 = os.path.join(p_apia2001, 'up1_dem.asc')
n_case2 = 'up20_out_2001'
d_o2 = os.path.join(p_apia2001, 'up20_dem.asc')

#Area extent
extent=[414360, 419410, 8468395, 8471551 ]

bas_raw = gpd.read_file(p_basins)
basins = Obtain_Basins(bas_raw, site)

dem = xr.open_dataset(p_dem)
dem=dem.sel(x=slice(extent[0],extent[1]), y=slice(extent[2],extent[3]))
dem['z']=(('y','x'),np.where(dem.z.values<0, np.nan, dem.z.values))
dem=dem.rename({'z':'elevation'})

dem2 = xr.open_dataset(p_dem)
dem2=dem2.sel(x=slice(extent[0],extent[1]), y=slice(extent[2],extent[3]))
dem2=dem2.rename({'z':'elevation'})
elevations2=dem2.elevation.values

Hydrograph

Information regarding the hydrograph of the Vaisigano flooding Event ocurred in Apia, Samoa in April 2001 has been extracted from the SOPAC Technical Report 338: A REVIEW OF FLOODING IN APIA, SAMOA, APRIL 2001, Developed by the Department of Physical Geography of the Macquarie University

Image('./resources/Vaisigano_hydrograph.png', width=1000)    

p_out = os.path.join(p_apia2001, n_case1, n_case1)
coefs = np.loadtxt(d_o1, max_rows=5, usecols=1)
ncols, nrows, xllcorner, yllcorner, cellsize = coefs[0], coefs[1], coefs[2], coefs[3], coefs[4]
out = xds_res_out(p_out, xllcorner, yllcorner, cellsize, ncols,nrows).squeeze()
out['max']=(('Y','X'),np.flipud(out['max'].values))

p_out = os.path.join(p_apia2001, n_case2, n_case2)
coefs = np.loadtxt(d_o2, max_rows=5, usecols=1)
ncols, nrows, xllcorner, yllcorner, cellsize = coefs[0], coefs[1], coefs[2], coefs[3], coefs[4]
out2 = xds_res_out(p_out, xllcorner, yllcorner, cellsize, ncols,nrows).squeeze()
out2['max']=(('Y','X'),np.flipud(out2['max'].values))

out=out.sel(X=slice(extent[0],extent[1]), Y=slice(extent[2],extent[3]))
out['max']=(('y','x'),np.where(out['max'].values<0.15, np.nan, out['max'].values))

out2=out2.sel(X=slice(extent[0],extent[1]), Y=slice(extent[2],extent[3]))
out2['max']=(('y','x'),np.where(out2['max'].values<0.15, np.nan, out2['max'].values))

Report Flooding

Image('./resources/Vaisigano_flood.png', width=1000)    

LisFlood-FP Flooding

ccmap = custom_cmap(100, 'gist_rainbow', 0.7,1, 'bwr_r', 0.5, 0.85)
hillshade = es.hillshade(elevations2,azimuth=140,altitude=40)

fig = plt.figure(figsize=[20,10])
gs = gridspec.GridSpec(1, 1)
ax = fig.add_subplot(gs[0, 0])
cmap = cmocean.cm.delta
cmap.set_under(color=cmap(0))

im1 = ax.pcolormesh(dem2.x.values, dem2.y.values, elevations2, cmap=cmap, 
                norm = colors.SymLogNorm(linthresh=0.1,vmin=-35,vmax=90),alpha=1)

ax.pcolormesh(dem2.x.values, dem2.y.values, hillshade, cmap="Greys", alpha=0.3)  # Shades
ax.pcolor(out.X.values, out.Y.values, out['max'].squeeze(), cmap=ccmap,vmin=0, vmax=2, alpha=.5)  # Shades
ax.pcolor(out2.X.values, out2.Y.values, out2['max'].squeeze(), cmap=ccmap,vmin=0, vmax=2, alpha=.5)  # Shades
ax.set_aspect('equal')
ax.set_xlim(extent[0], extent[1])
ax.set_ylim(extent[2], extent[3])
ax.set_xlabel('X_utm (m)',fontsize=16)
ax.set_ylabel('Y_utm (m)',fontsize=16)

Text(0, 0.5, 'Y_utm (m)')

Comparison

Image('./resources/Vaisigano_Superposition.png', width=1000)