<jupyter_start><jupyter_text>Hurricane Ike Maximum Water LevelsCompute the maximum water level during Hurricane Ike on a 9 million node triangular mesh storm surge model.  Plot the results with Datashader.<jupyter_code>%matplotlib inline

import warnings
warnings.filterwarnings("ignore")

from dask.distributed import Client, progress
from dask_kubernetes import KubeCluster
import xarray as xr
import numpy as np<jupyter_output><empty_output><jupyter_text>Start a dask cluster to crunch the data<jupyter_code>cluster = KubeCluster()
cluster.adapt(maximum=30);
cluster
client = Client(cluster)<jupyter_output><empty_output><jupyter_text>Read the data using the cloud-friendly zarr data format<jupyter_code>pangeo = 'ESIP-AWS-S3'
if pangeo=='PANGEO-GCS-S3':
    import gcsfs
    fs = gcsfs.GCSFileSystem(project='pangeo-181919')
    map = gcsfs.mapping.GCSMap('pangeo-data/rsignell/adcirc_test01', gcs=fs, check=False, create=False)
if pangeo=='ESIP-AWS-S3':
    import s3fs
    fs = s3fs.S3FileSystem(anon=True)
#    map = s3fs.S3Map('esip-pangeo/pangeo/adcirc/ike', s3=fs)
#    map = s3fs.S3Map('esip-pangeo/pangeo/adcirc/adcirc_01', s3=fs)
    map = s3fs.S3Map('esip-data-oregon/adcirc/adcirc_01', s3=fs)
%%time
ds = xr.open_zarr(map)
ds<jupyter_output><empty_output><jupyter_text>The zarr file was chunked with only 5 mb chunks.  So one might think that rechunking to 10 times as big chunks in the time dimension might be a good idea.  This turned out to be a terrible idea as dask did 4813 `zarr` operations, and then spent several minutes on `rechunk-merge-nanmax` operations before throwing an exception and dying.  Without rechunking, 60 workers do the job in 25 seconds.<jupyter_code># ds['zeta'] = ds['zeta'].chunk(chunks={'time':100, 'node':9228245})
ds['zeta']<jupyter_output><empty_output><jupyter_text>Note that unfortunately Xarray does not yet understand that `x` and `y` are coordinate variables, becuase we have a triangular mesh. Hopefully Xarray will gain some capabilities in this area in the future.  How many GB of sea surface height data do we have?<jupyter_code>ds['zeta'].nbytes/1.e9
max_var = ds['zeta'].max(dim='time').persist()
progress(max_var)<jupyter_output><empty_output><jupyter_text>Visualize data on mesh using Datashader<jupyter_code>import numpy as np
import datashader as dshade
import holoviews as hv
import geoviews as gv
import cartopy.crs as ccrs

import holoviews.operation.datashader as dshade

dshade.datashade.precompute = True

hv.extension('bokeh')
import warnings
warnings.filterwarnings("ignore")
import pandas as pd
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.WMTS('https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png')
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).opts(active_tools=['wheel_zoom'])<jupyter_output><empty_output><jupyter_text>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.<jupyter_code># 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])
%%time
ind = 7454940
ds['zeta'][:,ind].plot()
import matplotlib.pyplot as plt
plt.plot([1, 10, 20, 30, 60 ],[1/12,1/1.5,1/1.0, 1/0.75, 1/0.5],'o')
#cluster.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1325449618132545, 'fraction_numerical': 0.028578467602857847, 'mean_word_length': 5.7218543046357615, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 29, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '554124', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1536, 'n_tokens_neox': 1461, 'n_words': 395, 'path': 'hurricane_ike_water_level.ipynb', 'repo_name': 'cdi-workshops/pangeo_ocean_examples', 'size': '4059'}