hurricane_ike_water_levels.ipynb

{
 "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 using [HoloViz](https://holoviz.org/) TriMesh rendering with Datashader. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import fsspec\n",
    "from dask.distributed import Client, progress, LocalCluster"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Start a dask cluster to crunch the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#cluster_type = 'local'\n",
    "cluster_type = 'coiled'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if cluster_type == 'coiled':\n",
    "    import coiled\n",
    "    cluster = coiled.Cluster(\n",
    "        region=\"us-west-2\",\n",
    "        arm=True,\n",
    "        worker_vm_types=[\"t4g.small\"],  # cheap, small ARM instances, 2cpus, 2GB RAM\n",
    "        worker_options={'nthreads':2},\n",
    "        n_workers=30,\n",
    "        wait_for_workers=True,\n",
    "        compute_purchase_option=\"spot_with_fallback\",\n",
    "        software='esip-pangeo-arm',\n",
    "        workspace='esip-lab',\n",
    "        timeout=300   # leave cluster running for 5 min in case we want to use it again\n",
    "    )\n",
    "\n",
    "    client = cluster.get_client()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if cluster_type=='local':\n",
    "    cluster = LocalCluster()\n",
    "    client = Client(cluster)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if cluster_type=='gateway':\n",
    "    from dask_gateway import Gateway\n",
    "    gateway = Gateway()\n",
    "    cluster = gateway.new_cluster()\n",
    "    cluster.adapt(minimum=4, maximum=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if cluster_type=='kubernetes':\n",
    "    from dask_kubernetes import KubeCluster\n",
    "    cluster = KubeCluster()\n",
    "    cluster.adapt(minimum=4, maximum=20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read the data using the cloud-friendly zarr data format"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = fsspec.filesystem('s3', anon=True, endpoint_url='https://mghp.osn.xsede.org')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_dataset(fs.get_mapper('s3://rsignellbucket1/esip/adcirc/ike'), engine='zarr', chunks={'time':40})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds['zeta']"
   ]
  },
  {
   "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": [
    "Take the maximum over the time dimension.  This is the computationally intensive step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "max_var = ds['zeta'].max(dim='time').compute()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize data on mesh using HoloViz.org tools"
   ]
  },
  {
   "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",
    "import hvplot.xarray\n",
    "import holoviews.operation.datashader as dshade\n",
    "\n",
    "dshade.datashade.precompute = True\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "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.tile_sources.OSM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "value = 'max water level'\n",
    "label = '{} (m)'.format(value)\n",
    "trimesh = gv.TriMesh((tris, points), label=label)\n",
    "mesh = dshade.rasterize(trimesh).opts(\n",
    "              cmap='rainbow', colorbar=True, width=600, height=400)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tiles * mesh"
   ]
  },
  {
   "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": [
    "ds['zeta'][:,ind].hvplot(x='time', grid=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.9"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}