{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Hurricane Ike Maximum Water Levels\n", "Compute the maximum water level during Hurricane Ike on a 9 million node triangular mesh storm surge model. Plot the results with Datashader. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "from dask.distributed import Client, progress\n", "from dask_kubernetes import KubeCluster\n", "import xarray as xr\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Start a dask cluster to crunch the data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cluster = KubeCluster()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cluster.scale(25);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cluster" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "client = Client(cluster)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Read the data using the cloud-friendly zarr data format" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "pangeo = 'PANGEO-GCS-S3'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if pangeo=='ESIP-AWS-S3':\n", " import s3fs\n", " fs = s3fs.S3FileSystem(anon=True)\n", " map = s3fs.S3Map('rsignell/adcirc_test01', s3=fs)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "if pangeo=='PANGEO-GCS-S3':\n", " import gcsfs\n", " fs = gcsfs.GCSFileSystem(project='pangeo-181919')\n", " map = gcsfs.mapping.GCSMap('pangeo-data/rsignell/adcirc_test01', gcs=fs, check=False, create=False)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds = xr.open_zarr(map)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "ds" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds['zeta']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "How many GB of sea surface height data do we have?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds['zeta'].nbytes/1.e9" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want to calculate the maximum water level at every cell over all time steps" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ds['zeta'].max(dim='time')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Dask will read data only when necessary, so nothing has been loaded yet. In the next cell, however, we use `persist` to persist the data to the Dask worker nodes, and watch the progress" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "max_var = ds['zeta'].max(dim='time').persist()\n", "progress(max_var)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Visualize data on mesh using Datashader" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import datashader as dshade\n", "import holoviews as hv\n", "import geoviews as gv\n", "import cartopy.crs as ccrs\n", "\n", "from holoviews.operation.datashader import datashade, rasterize\n", "\n", "datashade.precompute = True\n", "\n", "hv.extension('bokeh')\n", "%opts Image RGB VectorField [width=800 height=600]" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "warnings.filterwarnings(\"ignore\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import pandas as pd" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "v = np.vstack((ds['x'], ds['y'], max_var)).T\n", "verts = pd.DataFrame(v, columns=['x','y','vmax'])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "points = gv.operation.project_points(gv.Points(verts, vdims=['vmax']))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tris = pd.DataFrame(ds['element'].values.astype('int')-1, columns=['v0','v1','v2'])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "tiles = gv.WMTS('https://maps.wikimedia.org/osm-intl/{Z}/{X}/{Y}@2x.png')\n", "value = 'max water level'\n", "label = '{} (m)'.format(value)\n", "trimesh = gv.TriMesh((tris, points), label=label)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "%%opts Image [colorbar=True] (cmap='rainbow')\n", "\n", "meshes = rasterize(trimesh,aggregator=dshade.mean('vmax'))\n", "tiles * meshes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Extract a time series at a specified lon, lat location" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# find the indices of the points in (x,y) closest to the points in (xi,yi)\n", "def nearxy(x,y,xi,yi):\n", " ind = np.ones(len(xi),dtype=int)\n", " for i in range(len(xi)):\n", " dist = np.sqrt((x-xi[i])**2+(y-yi[i])**2)\n", " ind[i] = dist.argmin()\n", " return ind" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#just offshore of Galveston\n", "lat = 29.2329856\n", "lon = -95.1535041" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ind = nearxy(ds['x'].values,ds['y'].values,[lon], [lat])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%time\n", "ds['zeta'][:,ind].plot()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.6" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": { "01fcfe3453574be4ac3d37e28425d0b0": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.2.0", "model_name": "HBoxModel", "state": { "children": [ "IPY_MODEL_fad4f62c9399419993153fb961830488", "IPY_MODEL_bdfc56e2940c4c739ff1d5faeec40d7f", "IPY_MODEL_14478fbdfc0a4fa9a41fb87209cfb76b" ], "layout": "IPY_MODEL_dc161608c8bc4a019de49bad6cc0cb2f" } }, "061e3b19ad0944bb98c6b44bd10a989c": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.2.0", "model_name": "DescriptionStyleModel", "state": { "description_width": "" } }, "07d213ec9c774483937f447bf3dd86ae": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.2.0", "model_name": "FloatProgressModel", "state": { "bar_style": "success", "layout": "IPY_MODEL_f860565f03c9460a9fcd197828d4a043", "max": 1, "style": "IPY_MODEL_86cda597474343f990730a66d71dce98", "value": 1 } }, "13c18d134d3e495b9dff310d81b0c2cc": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.2.0", "model_name": "HBoxModel", "state": { "children": [ "IPY_MODEL_fbf901b76b014589bcd19fb800b66b6a", "IPY_MODEL_5b148f98c7da47678f5e38de89b5c215", "IPY_MODEL_e667ec2530bb4d9bb801700b164a1ccf" ], "layout": "IPY_MODEL_6d01a84902be460aa27a72b0da24383a" } }, "14478fbdfc0a4fa9a41fb87209cfb76b": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.2.0", "model_name": "HTMLModel", "state": { "layout": "IPY_MODEL_6b3abb3e1d9f4a9685960ac30e86037b", "style": "IPY_MODEL_4f1f297a1c524a489d284d1f5ebf9f18", "value": "