surge_updates.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "<table>\n",
    " <tr align=left><td><img align=left src=\"https://i.creativecommons.org/l/by/4.0/88x31.png\">\n",
    " <td>Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved MIT license. (c) Kyle T. Mandli</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# GeoClaw Storm Surge Modeling:  Updates and New Features"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Thank You to Contributors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Storm Object\n",
    "\n",
    "> This object contains a time series of time data that describe a particular\n",
    "    storm.  This includes the attributes below and the ability to read from\n",
    "    multiple sources for data such as the U.S. National Hurricane Center (NHC),\n",
    "    the Japanese Meterological Agency (JMA), and the Indian Meteorlogical\n",
    "    Department (IMD).  This class can then write out in any of these formats,\n",
    "    construct the wind and pressure fields using a supported parameterized\n",
    "    model, or output the GeoClaw supported storm format used for running storm\n",
    "    surge simulations.\n",
    "\n",
    "**Important Point** - Python data container for time series description of a storm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Core Capabilities\n",
    "\n",
    " - `read` in lots of storm formats (e.g. *ATCF*, *HURDAT*, *JMA*)\n",
    "\n",
    " - `write` out in some storm formats (e.g. *GeoClaw*, *ATCF*, *HURDAT*)\n",
    "\n",
    " - `plot` storm track and intensity\n",
    "\n",
    " - `category` of storm given a categorization scheme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Utilities\n",
    "\n",
    " - `construct_fields` of a storm given a parameterization (e.g. Holland 1980, CLE 2015)\n",
    " - Others"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "### Example Use\n",
    "```python\n",
    "# Storm parameters - Parameterized storm (Holland 1980)\n",
    "data.storm_specification_type = 'holland80'  # (type 1)\n",
    "data.storm_file = os.path.expandvars(os.path.join(os.getcwd(), 'ike.storm'))\n",
    "```\n",
    "\n",
    "```python\n",
    "# Convert ATCF data to GeoClaw format\n",
    "clawutil.data.get_remote_file(\n",
    "               \"http://ftp.nhc.noaa.gov/atcf/archive/2008/bal092008.dat.gz\")\n",
    "atcf_path = os.path.join(scratch_dir, \"bal092008.dat\")\n",
    "# Note that the get_remote_file function does not support gzip files which\n",
    "# are not also tar files.  The following code handles this\n",
    "with gzip.open(\".\".join((atcf_path, 'gz')), 'rb') as atcf_file,    \\\n",
    "        open(atcf_path, 'w') as atcf_unzipped_file:\n",
    "    atcf_unzipped_file.write(atcf_file.read().decode('ascii'))\n",
    "```\n",
    "\n",
    "```python\n",
    "# Read in unzipped file\n",
    "ike = Storm(path=atcf_path, file_format=\"ATCF\")\n",
    "```\n",
    "\n",
    "```python\n",
    "# Calculate landfall time - Need to specify as the file above does not\n",
    "# include this info (9/13/2008 ~ 7 UTC)\n",
    "ike.time_offset = datetime.datetime(2008, 9, 13, 7)\n",
    "```\n",
    "\n",
    "```python\n",
    "# Write storm data to GeoClaw format\n",
    "ike.write(data.storm_file, file_format='geoclaw')\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "**GeoClaw Storm Format**\n",
    "```bash\n",
    "53\n",
    "2008-09-13T07:00:00\n",
    "\n",
    "    -1.04040000e+06     -3.70000000e+01      1.72000000e+01      1.54333332e+01      1.66680000e+05      1.00600000e+05      4.63000000e+05\n",
    "    -1.01880000e+06     -3.84000000e+01      1.73000000e+01      1.80055554e+01      1.66680000e+05      1.00500000e+05      4.63000000e+05\n",
    "    -9.97200000e+05     -3.99000000e+01      1.75000000e+01      2.31499998e+01      1.66680000e+05      1.00300000e+05      4.63000000e+05\n",
    "    -9.75600000e+05     -4.13000000e+01      1.78000000e+01      2.31499998e+01      3.70400000e+04      1.00200000e+05      4.63000000e+05\n",
    "    -9.54000000e+05     -4.28000000e+01      1.82000000e+01      2.31499998e+01      3.70400000e+04      1.00000000e+05      5.09300000e+05\n",
    "    -9.32400000e+05     -4.43000000e+01      1.87000000e+01      2.57222220e+01      1.66680000e+05      9.99000000e+04      4.63000000e+05\n",
    "...\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import gzip\n",
    "import datetime\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from clawpack.geoclaw.surge.storm import Storm\n",
    "import clawpack.clawutil as clawutil\n",
    "\n",
    "# Scratch directory for storing topo and dtopo files:\n",
    "scratch_dir = os.path.join(os.environ[\"CLAW\"], 'geoclaw', 'scratch')\n",
    "\n",
    "# Convert ATCF data to GeoClaw format\n",
    "clawutil.data.get_remote_file(\"http://ftp.nhc.noaa.gov/atcf/archive/2008/bal092008.dat.gz\")\n",
    "atcf_path = os.path.join(scratch_dir, \"bal092008.dat\")\n",
    "# Note that the get_remote_file function does not support gzip files which\n",
    "# are not also tar files.  The following code handles this\n",
    "with gzip.open(\".\".join((atcf_path, 'gz')), 'rb') as atcf_file,    \\\n",
    "        open(atcf_path, 'w') as atcf_unzipped_file:\n",
    "    atcf_unzipped_file.write(atcf_file.read().decode('ascii'))\n",
    "\n",
    "# Uncomment/comment out to use the old version of the Ike storm file\n",
    "# ike = Storm(path=\"old_ike.storm\", file_format=\"ATCF\")\n",
    "ike = Storm(path=atcf_path, file_format=\"ATCF\")\n",
    "\n",
    "# Calculate landfall time - Need to specify as the file above does not\n",
    "# include this info (9/13/2008 ~ 7 UTC)\n",
    "ike.time_offset = datetime.datetime(2008, 9, 13, 7)\n",
    "\n",
    "# Plot - Incorporated into the head branch\n",
    "fig = plt.figure()\n",
    "axes = fig.add_subplot(1, 1, 1)\n",
    "category_color = {5: 'red',\n",
    "                  4: 'yellow',\n",
    "                  3: 'orange',\n",
    "                  2: 'green',\n",
    "                  1: 'blue',\n",
    "                  0: 'gray'}\n",
    "category = ike.category(categorization=\"NHC\")\n",
    "longitude = ike.eye_location[:, 0]\n",
    "latitude = ike.eye_location[:, 1]\n",
    "for i in range(len(longitude)):\n",
    "    color = category_color[category[i]]\n",
    "    axes.plot(longitude[i:i + 2], latitude[i:i + 2], color=color)\n",
    "\n",
    "axes.set_title(\"Hurricane Ike Track\")\n",
    "axes.set_xlabel(\"Longitude\")\n",
    "axes.set_ylabel(\"Latitude\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Paramaterized Wind Models\n",
    "\n",
    "**Parameterizations:** Holland 1980, Holland 2010, CLE, SLOSH, Rankine, Modified Rankine, deMaria, Willoughby\n",
    "\n",
    "Added to both the Python and Fortran code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import numpy\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Parameters\n",
    "N = 1000\n",
    "radius = 100e3\n",
    "Pn = 1005\n",
    "Pc = 950\n",
    "A = 23.0\n",
    "B = 1.5\n",
    "rho_air = 1.15\n",
    "OMEGA = 7.2722052166430395e-5\n",
    "THETA_0 = 0.52359877559829882 \n",
    "f = 2.0 * OMEGA * numpy.sin(THETA_0)\n",
    "f = 0.0\n",
    "\n",
    "# Evaluate profiles\n",
    "x = numpy.concatenate((numpy.linspace(-radius, -0.01, N), \n",
    "                       numpy.linspace(0.01, radius, N)), axis=0)\n",
    "r = numpy.abs(x) * 1e-3\n",
    "p = Pc + (Pn - Pc) * numpy.exp(-A/(r)**B)\n",
    "C = 1e1**2 *A * B * (Pn - Pc) / (rho_air)\n",
    "v = numpy.sqrt(C * numpy.exp(-A / r**B) / r**B + r**2 * f**2 / 4.0) - r * f / 2.0\n",
    "\n",
    "fig = plt.figure()\n",
    "fig.set_figwidth(fig.get_figwidth() * 2)\n",
    "\n",
    "# Convert to kms\n",
    "x /= 1e3\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 1)\n",
    "axes.plot(x, v)\n",
    "axes.set_title(\"Wind Velocity Profile\")\n",
    "axes.set_xlabel('km')\n",
    "axes.set_ylabel('m/s')\n",
    "axes.set_xlim([numpy.min(x), numpy.max(x)])\n",
    "axes.set_ylim([0.0, numpy.max(v) + 5])\n",
    "\n",
    "axes = fig.add_subplot(1, 2, 2)\n",
    "axes.plot(x, p)\n",
    "axes.set_title(\"Pressure Profile\")\n",
    "axes.set_xlabel('km')\n",
    "axes.set_ylabel('mb')\n",
    "axes.set_xlim([numpy.min(x), numpy.max(x)])\n",
    "axes.set_ylim([Pc - 5, Pn + 5])\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Examples Added\n",
    "\n",
    " - Synthetic basin example\n",
    " - [Storm Surge Tutorial](https://github.com/clawpack/apps/tree/master/notebooks/geoclaw/katrina) - https://github.com/clawpack/apps/tree/master/notebooks/geoclaw/katrina\n",
    " - [Full Examples](https://github.com/mandli/surge-examples) - https://github.com/mandli/surge-examples - Please Contribute!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Ensembles of Storms\n",
    "\n",
    "Create and manage an [ensemble of storms](https://github.com/mandli/batch/).\n",
    "\n",
    "Tutorial is incoming.\n",
    "\n",
    "https://github.com/mandli/batch/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bugs Squashed\n",
    "\n",
    " - Recording storm fields at gauges"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Thanks!\n",
    "\n",
    " - Ask questions via [google groups](groups.google.com), [slack](clawpack.slack.com) or [github](www.github.com/clawpack).\n",
    " - Contribute!"
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}