Hurricane Ike Maximum Water Levels

Compute maximum water levels during Hurricane Ike on a 9 million node triangular mesh ADCIRC storm surge model. Visualize the results using HoloViz TriMesh rendering with Datashader.

[ ]:
import xarray as xr
import numpy as np
import pandas as pd
import fsspec

Start a dask cluster to crunch the data

[ ]:
from dask.distributed import Client
[ ]:
from dask_gateway import Gateway
gateway = Gateway()
cluster = gateway.new_cluster()
[ ]:
#from dask_kubernetes import KubeCluster
#cluster = KubeCluster()
[ ]:
cluster
[ ]:
cluster.adapt(minimum=4, maximum=20);
[ ]:
client = Client(cluster)

Read the data using the cloud-friendly zarr data format

[ ]:
ds = xr.open_zarr(fsspec.get_mapper('s3://pangeo-data-uswest2/esip/adcirc/ike', anon=False, requester_pays=True))
[ ]:
#ds = xr.open_zarr(fsspec.get_mapper('gcs://pangeo-data/rsignell/adcirc_test01'))
[ ]:
ds['zeta']

How many GB of sea surface height data do we have?

[ ]:
ds['zeta'].nbytes/1.e9

Take the maximum over the time dimension and persist the data on the workers to use later. This is the computationally intensive step.

[ ]:
max_var = ds['zeta'].max(dim='time').persist()

Visualize data on mesh using HoloViz.org tools

[ ]:
import numpy as np
import datashader as dshade
import holoviews as hv
import geoviews as gv
import cartopy.crs as ccrs
import hvplot.xarray
import holoviews.operation.datashader as dshade

dshade.datashade.precompute = True
hv.extension('bokeh')
[ ]:
v = np.vstack((ds['x'], ds['y'], max_var)).T
verts = pd.DataFrame(v, columns=['x','y','vmax'])
[ ]:
points = gv.operation.project_points(gv.Points(verts, vdims=['vmax']))
[ ]:
tris = pd.DataFrame(ds['element'].values.astype('int')-1, columns=['v0','v1','v2'])
[ ]:
tiles = gv.tile_sources.OSM
[ ]:
value = 'max water level'
label = '{} (m)'.format(value)
trimesh = gv.TriMesh((tris, points), label=label)
mesh = dshade.rasterize(trimesh).opts(
              cmap='rainbow', colorbar=True, width=600, height=400)
[ ]:
tiles * mesh

Extract a time series at a specified lon, lat location

Because Xarray does not yet understand that x and y are coordinate variables on this triangular mesh, we create our own simple function to find the closest point. If we had a lot of these, we could use a more fancy tree algorithm.

[ ]:
# find the indices of the points in (x,y) closest to the points in (xi,yi)
def nearxy(x,y,xi,yi):
    ind = np.ones(len(xi),dtype=int)
    for i in range(len(xi)):
        dist = np.sqrt((x-xi[i])**2+(y-yi[i])**2)
        ind[i] = dist.argmin()
    return ind
[ ]:
#just offshore of Galveston
lat = 29.2329856
lon = -95.1535041
[ ]:
ind = nearxy(ds['x'].values,ds['y'].values,[lon], [lat])
[ ]:
ds['zeta'][:,ind].hvplot(x='time', grid=True)