Fitting a radial profile to tropical cyclone observations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting a raidal profile to tropical cyclone observations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook enables users to load an analysis track file and interactively fit a radial profile to the range of wind radii provided in the track file. The aim is to improve the quality of the simulated wind field by improving the fit of parameters in the profile (namely the Holland $\\beta$ value and the environmental pressure).\n",
    "\n",
    "At a minimum, it requires estimates of the central pressure and the radius to maximum winds. In the absence of pressure at the outermost closed isobar, the daily long-term-mean sea level pressure at the location of the storm is used. \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preliminaries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the notebook with the required modules for plotting, analysis and the fiddly things like file loaders. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%matplotlib inline\n",
    "\n",
    "from __future__ import print_function # py 2.7 compat.\n",
    "\n",
    "import sys\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use [``matplotlib``](http://matplotlib.org/index.html) for plotting the data up. We use [``basemap``](http://matplotlib.org/basemap/) to plot the track on a simple map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.collections import LineCollection\n",
    "from matplotlib.colors import Normalize, BoundaryNorm, ListedColormap\n",
    "from matplotlib.cm import get_cmap\n",
    "from matplotlib.dates import HourLocator, DateFormatter\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "from mpl_toolkits.basemap import cm as basemapcm\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import seaborn; seaborn.set()\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "# Import widgets for interactive notebook\n",
    "from IPython.html.widgets import interact, fixed\n",
    "from IPython.html import widgets\n",
    "\n",
    "from Utilities.metutils import convert\n",
    "from Utilities.track import Track\n",
    "from PlotInterface.tracks import TrackMapFigure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ``NotebookFinder`` function allows me to import other IPython notebooks as Python modules. I use another notebook to set up the interactive widget for loading a file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from NoteBookFinder import NotebookFinder\n",
    "sys.meta_path.append(NotebookFinder())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We build a widget to allow selection of a track file. This is from [Interactive Widgets](http://nbviewer.ipython.org/github/adrn/ipython/blob/2.x/examples/Interactive%20Widgets/File%20Upload%20Widget.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import FileWidgets\n",
    "\n",
    "FileLoader = FileWidgets.FileLoaderWidget()\n",
    "FileLoader"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the existing wind modules from TCRM to generate the profile. These are part of the ``wind.windmodels`` module. For comparison, we use the Powell and Willoughby profiles as well. The former calculates $\\beta$ as a function of maximum wind speed and latitude, while the latter uses maximum wind speed only."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The track data needs to be loaded from a csv file. Here's the challenge - different track files have different formats, with different field names for the same data. Some files have fields that others don't, while they all have different date formats and potentially different units for data like pressure and wind speed. \n",
    "\n",
    "For initial prototyping, I'm using the format of the BoM's Analysis Track data for TC Marcia, provided by David Grant (Qld Regional Office Severe Weather Section). The intensity and track data is preliminary, based on operational estimates and subject to change following post analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "COLNAMES = ['Datetime', 'Latitude', 'Longitude', 'Symbol', 'Category','CentralPressure', 'PressureOCI', 'RadiusOCI', 'Radius1000hPa',\n",
    "            'RadiusMaxWinds', 'MeanWind', 'WindSpeed', 'VerticalExtent', 'Uncertainty', 'FinalT', 'CurrentIntensity',\n",
    "            'R34', 'R41', 'R48', 'R64', 'NER34', 'SER34', 'SWR34', 'NWR34', 'NER41', 'SER41', 'SWR41', \n",
    "            'NWR41', 'NER48', 'SER48', 'SWR48', 'NWR48', 'NER64', 'SER64', 'SWR64', 'NWR64']\n",
    "COLTYPES = [datetime, 'f8', 'f8', 'i', 'i', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'i', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', \n",
    "            'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f', 'f']\n",
    "COLUNITS = ('', 'degree', 'degree', '', '', 'hPa', 'hPa', 'nm', 'nm',\n",
    "            'nm', 'kts', 'kts', '', 'nm', '', '', \n",
    "            'nm', 'nm', 'nm', 'nm', 'nm', 'nm', 'nm', 'nm', \n",
    "            'nm', 'nm', 'nm', 'nm', 'nm', 'nm', 'nm', 'nm', \n",
    "            'nm', 'nm', 'nm', 'nm')\n",
    "DATEFORMAT = \"%Y-%m-%dT%H:%M:%SZ\"\n",
    "dtype = np.dtype({'names':COLNAMES, 'formats':COLTYPES})\n",
    "converters = {\n",
    "    0: lambda s: datetime.strptime(s.strip(), DATEFORMAT),\n",
    "    5: lambda s: convert(float(s.strip() or 0), COLUNITS[5], 'Pa'),\n",
    "    6: lambda s: convert(float(s.strip() or 0), COLUNITS[6], 'Pa'),\n",
    "    7: lambda s: convert(float(s.strip() or 0), COLUNITS[7], 'km'),\n",
    "    8: lambda s: convert(float(s.strip() or 0), COLUNITS[8], 'km'),\n",
    "    9: lambda s: convert(float(s.strip() or 0), COLUNITS[9], 'km'),\n",
    "    10: lambda s: convert(float(s.strip()), COLUNITS[10], 'mps'),\n",
    "    11: lambda s: convert(float(s.strip()), COLUNITS[11], 'mps'),\n",
    "    13: lambda s: convert(float(s.strip()), COLUNITS[13], 'km'),\n",
    "    16: lambda s: convert(float(s.strip()), COLUNITS[16], 'km'),\n",
    "    16: lambda s: convert(float(s.strip()), COLUNITS[16], 'km'),\n",
    "    17: lambda s: convert(float(s.strip()), COLUNITS[17], 'km'),\n",
    "    18: lambda s: convert(float(s.strip()), COLUNITS[18], 'km'),\n",
    "    19: lambda s: convert(float(s.strip()), COLUNITS[19], 'km'),\n",
    "    20: lambda s: convert(float(s.strip()), COLUNITS[20], 'km'),\n",
    "    21: lambda s: convert(float(s.strip()), COLUNITS[21], 'km'),\n",
    "    22: lambda s: convert(float(s.strip()), COLUNITS[22], 'km'),\n",
    "    23: lambda s: convert(float(s.strip()), COLUNITS[23], 'km'),\n",
    "    24: lambda s: convert(float(s.strip()), COLUNITS[24], 'km'),\n",
    "    25: lambda s: convert(float(s.strip()), COLUNITS[25], 'km'),\n",
    "    26: lambda s: convert(float(s.strip()), COLUNITS[26], 'km'),\n",
    "    27: lambda s: convert(float(s.strip()), COLUNITS[27], 'km'),\n",
    "    28: lambda s: convert(float(s.strip()), COLUNITS[28], 'km'),\n",
    "    29: lambda s: convert(float(s.strip()), COLUNITS[29], 'km'),\n",
    "    30: lambda s: convert(float(s.strip()), COLUNITS[30], 'km'),\n",
    "    31: lambda s: convert(float(s.strip()), COLUNITS[31], 'km'),\n",
    "    32: lambda s: convert(float(s.strip()), COLUNITS[32], 'km'),\n",
    "    33: lambda s: convert(float(s.strip()), COLUNITS[33], 'km'),\n",
    "    34: lambda s: convert(float(s.strip()), COLUNITS[34], 'km')\n",
    "}\n",
    "delimiter = ','\n",
    "skip_header = 1\n",
    "usecols = tuple(range(36))\n",
    "missing_value = \"NaN\"\n",
    "filling_values = 0\n",
    "filename = FileLoader.filename\n",
    "\n",
    "data = np.genfromtxt(filename, dtype, delimiter=delimiter, skip_header=skip_header, converters=converters,\n",
    "                     missing_values=missing_value, filling_values=filling_values, usecols=usecols, autostrip=True)\n",
    "\n",
    "times = data['Datetime']\n",
    "timelist = [d.strftime('%Y-%m-%d %H:%M') for d in data['Datetime']]\n",
    "\n",
    "diffs = [d1 - d0 for d0, d1 in zip(times[:-1], times[1:])]\n",
    "dt =  np.array([0.0] + [round(d.days * 24. + d.seconds / 3600.)\n",
    "                             for d in diffs], 'f')\n",
    "index=np.zeros(len(times), 'i')\n",
    "from Utilities.loadData import getSpeedBearing\n",
    "\n",
    "speed, bearing = getSpeedBearing(index, data['Longitude'], data['Latitude'], dt)\n",
    "speed = convert(speed, 'kmh', 'mps')\n",
    "age = np.cumsum(dt)\n",
    "lon = data['Longitude']\n",
    "lat = data['Latitude']\n",
    "penv = data['PressureOCI']\n",
    "pcentre = data['CentralPressure']\n",
    "rmax = data['RadiusMaxWinds']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data is now stored as a ``numpy.recarray`` object, which makes it easier to work with. \n",
    "\n",
    "Check the date range:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print(\"First time: {0}\".format(timelist[0]))\n",
    "print(\"Final time: {0}\".format(timelist[-1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the track on a map to give some context. These two functions enable you to colourize the line segments (e.g. based on intensity). The classification is based on the estimated maximum gust wind speed and uses the Australian TC intnsity scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def make_segments(x, y):\n",
    "    points = np.array([x, y]).T.reshape(-1, 1, 2)\n",
    "    segments = np.concatenate([points[:-1], points[1:]], axis=1)\n",
    "\n",
    "    return segments\n",
    "\n",
    "def colorline(x, y, z=None, linewidth=1.0, alpha=1.0):\n",
    "    if z is None:\n",
    "        z = np.linspace(0.0, 1.0, len(x))\n",
    "\n",
    "    if not hasattr(z, '__iter__'):\n",
    "        z = np.array([z])\n",
    "\n",
    "    z = np.asarray(z)\n",
    "\n",
    "    segments = make_segments(x, y)\n",
    "    cmap = ListedColormap(['0.75', '#0FABF6', '#0000FF',\n",
    "                            '#00FF00', '#FF8100', '#ff0000'])\n",
    "    norm = BoundaryNorm([0, 17.5, 24.5, 32.5, 44.2, 55.5, 1000.], cmap.N)\n",
    "    lc = LineCollection(segments, array=z, cmap=cmap, \n",
    "                        norm=norm, linewidth=linewidth, alpha=alpha)\n",
    "    \n",
    "    ax = plt.gca()\n",
    "    ax.add_collection(lc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The map is set up to cover the extent of the track, with a 2$\\^{circ}$ margin around it. I haven't added parallels and meridians at this time (but could be a good exercise to do). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "minLon = np.floor(lon.min()) - 2\n",
    "maxLon = np.ceil(lon.max()) + 2\n",
    "minLat = np.floor(lat.min()) - 2\n",
    "maxLat = np.ceil(lat.max()) + 2\n",
    "\n",
    "xx = np.arange(minLon, maxLon + 0.1, 0.1)\n",
    "yy = np.arange(minLat, maxLat + 0.1, 0.1)\n",
    "xgrid, ygrid = np.meshgrid(xx, yy)\n",
    "\n",
    "m = Basemap(projection='mill', llcrnrlon=minLon, llcrnrlat=minLat, \n",
    "            urcrnrlon=maxLon, urcrnrlat=maxLat, resolution='h')\n",
    "m.drawcoastlines()\n",
    "m.drawstates()\n",
    "m.drawcountries()\n",
    "m.fillcontinents('coral')\n",
    "x,y = m(lon, lat)\n",
    "colorline(x,y,z=data['WindSpeed'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a simple function to calculate $\\beta$ using the method of Powell *et al.* (2005)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def powell_beta(lat, rmax):\n",
    "    rm = ma.array(rmax, mask=np.isnan(rmax))\n",
    "    beta = 1.881093 - 0.010917 * np.abs(lat) - 0.005567 * rm\n",
    "    return beta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the main parameters over the lifetime of the event. This can help identify the times where sufficient parameters are available to do some sort of profile fitting. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plottimes(data, beta=None):\n",
    "    dayLocator = DayLocator()\n",
    "    hourLocator = HourLocator(interval=6)\n",
    "    dateFormat = DateFormatter('%H:%MZ\\n%Y-%m-%d')\n",
    "    times = data['Datetime']\n",
    "    fig = plt.figure(figsize=(16,10))\n",
    "    ax1 = plt.subplot(2,1,1)\n",
    "    ax1.plot(times, convert(data['CentralPressure'], \"Pa\", \"hPa\"), linestyle='-', \n",
    "             color='k', linewidth=2, label=\"Central pressure (hPa)\")\n",
    "    ax1.scatter(times, convert(data['CentralPressure'], \"Pa\", \"hPa\"), marker='+', color='k' )\n",
    "\n",
    "    ax1.plot(times, convert(data['PressureOCI'], \"Pa\", \"hPa\"), linestyle='-', \n",
    "             color='0.75', linewidth=2, label=\"Environmental pressure (hPa)\")\n",
    "    ax1.scatter(times, convert(data['PressureOCI'], \"Pa\", \"hPa\"), marker='+', color='0.75' )\n",
    "\n",
    "    ax1.set_ylabel(\"Pressure (hPa)\")\n",
    "    ax1.set_xlabel(\"Date\")\n",
    "    ax1.xaxis.set_major_locator(dayLocator)\n",
    "    ax1.xaxis.set_minor_locator(hourLocator)\n",
    "    l = ax1.legend(loc='lower left')\n",
    "    for t in l.get_texts():\n",
    "        t.set_fontsize('x-small')\n",
    "    ax1.grid(axis='x')\n",
    "\n",
    "    ax2 = ax1.twinx()\n",
    "    p = ax2.plot(times, data['MeanWind'], linestyle='-', color='r', \n",
    "                 linewidth=2, label=\"Mean wind speed (m/s)\")\n",
    "    ax2.scatter(times, data['MeanWind'], marker='+', color='r')\n",
    "    \n",
    "    ax2.plot(times, data['WindSpeed'], linestyle='-', color='pink', \n",
    "             linewidth=2, label=\"Maximum wind speed (m/s)\")\n",
    "    ax2.scatter(times, data['WindSpeed'], marker='+', color='pink')\n",
    "    \n",
    "    ax2.plot(times, data['WindSpeed'] - 0.7*speed, linestyle=':', color='r',\n",
    "             linewidth=1, label='Modified maximum wind speed')\n",
    "    \n",
    "    ax2.plot(times,data['RadiusMaxWinds'], linestyle='-', color='b', \n",
    "             linewidth=2, label=\"Radius to maximum wind (km)\")\n",
    "\n",
    "    ax2.scatter(times,data['RadiusMaxWinds'], marker='+', color='b')\n",
    "    ax2.set_ylabel(\"Wind speed (m/s)\\nRadius (km)\")\n",
    "    ax2.xaxis.set_major_locator(dayLocator)\n",
    "    ax2.xaxis.set_minor_locator(hourLocator)\n",
    "    ax2.xaxis.set_major_formatter(dateFormat)\n",
    "    plt.xlim((times[0], times[-1]))\n",
    "    l = ax2.legend(loc='lower right')\n",
    "    for t in l.get_texts():\n",
    "        t.set_fontsize('x-small')\n",
    "    plt.grid(axis='both')\n",
    "    \n",
    "    if beta is not None:\n",
    "        pbeta = powell_beta(data['Latitude'], data['RadiusMaxWinds'])\n",
    "        #print(len(pbeta))\n",
    "        #print(len(times))\n",
    "        ax3 = plt.subplot(2,1,2)\n",
    "        ax3.plot(times, beta, linestyle='-', color='k', linewidth=2, label=r'$\\beta$')\n",
    "        ax3.scatter(times, beta, marker='+', color='k')\n",
    "        ax3.plot(times, pbeta, linestyle='-', color='0.75', linewidth=2, label=r\"Powell's $\\beta$\")\n",
    "        ax3.scatter(times, pbeta, marker='+', color='0.75')\n",
    "        ax3.set_ylabel(r'$\\beta$')\n",
    "        ax3.set_xlabel(\"Date\")\n",
    "        ax3.xaxis.set_major_locator(dayLocator)\n",
    "        ax3.xaxis.set_minor_locator(hourLocator)\n",
    "        ax3.xaxis.set_major_formatter(dateFormat)\n",
    "        l = ax3.legend(loc='upper left')\n",
    "        for t in l.get_texts():\n",
    "            t.set_fontsize('x-small')\n",
    "        ax3.grid()\n",
    "        plt.xlim((times[0], times[-1]))\n",
    "        \n",
    "    plt.show()\n",
    "\n",
    "        \n",
    "plottimes(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The modified maximum wind speed is reduced by 0.7 times the forward speed of the cyclone. This is to account for the wavenumber-1 asymmetry induced in the vortex flow by the forward motion. In reality, it is likely not a constant value, but some boundary layer models use this simple approximation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now to the task at hand. We're trying to fit the Holland profile as best as possible to the recorded wind radii. We perform this at each time step individually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import wind.windmodels as windmodels\n",
    "except ImportError:\n",
    "    print(\"Cannot import windmodels. Check that the TCRM folder is in the PYTHONPATH\")\n",
    "    \n",
    "holland = windmodels.HollandWindProfile\n",
    "powell = windmodels.PowellWindProfile\n",
    "willoughby = windmodels.WilloughbyWindProfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def get_wind(profile, r):\n",
    "    try:\n",
    "        V = profile.velocity(r)\n",
    "    except AssertionError:\n",
    "        print(\"Missing required fields\")\n",
    "        return np.zeros(len(r))\n",
    "    else:\n",
    "        return V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up functions to caculate the residuals of the model fit compared to the analysis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def residuals(beta, r, v, rmax, vmax, lat, lon, penv, pcentre):\n",
    "    \"\"\"\n",
    "    Return the residuals between the model fit and the analysis\n",
    "    \"\"\"\n",
    "    #assert len(r)==len(v)\n",
    "    robs = np.zeros(len(r)+1)\n",
    "    vobs = np.zeros(len(v)+1)\n",
    "    robs[1:] = r\n",
    "    robs[0] = rmax\n",
    "    vobs[1:] = v\n",
    "    vobs[0] = vmax\n",
    "    h = holland(lat, lon, penv, pcentre, rmax, beta)\n",
    "    try:\n",
    "        vmodel = h.velocity(robs)\n",
    "    except AssertionError:\n",
    "        return np.zeros(len(robs))\n",
    "    \n",
    "    err = abs(vmodel) - vobs\n",
    "    return err\n",
    "\n",
    "from scipy.optimize import leastsq\n",
    "\n",
    "def minimise(r, v, rmax, vmax, lat, lon, penv, pcentre, beta=1.5):\n",
    "    \"\"\"\n",
    "    Provide a functional interface to the leastsq function from scipy \n",
    "    and return the parameter value.\n",
    "    \n",
    "    \"\"\"\n",
    "    plsq = leastsq(residuals, [beta], args=(r, v, rmax, vmax, lat, lon, penv, pcentre))\n",
    "    return plsq[0]\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second panel here plots the optimum $\\beta$ value from Holland's profile (see Holland, 1980), estimated by fitting the profile to the analysed radii to various wind speeds (gales, strong gale, storm and hurricane-force wind speeds), as well as the radius to maximum winds.  We use least-squares minimisation to obtain the $\\beta$ value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta = 1.8*np.ones(len(times))\n",
    "\n",
    "\n",
    "for i, t in enumerate(times):\n",
    "    pc = pcentre[i]\n",
    "    pe = penv[i]\n",
    "    rm = rmax[i]\n",
    "    roci = data['RadiusOCI'][i]\n",
    "    vm = data['WindSpeed'][i] - 0.7*speed[i]\n",
    "    loni = lon[i]\n",
    "    lati = lat[i]\n",
    "\n",
    "    if np.isnan(data['R34'][i]):\n",
    "        r34 = np.max([data['NER34'][i], data['NWR34'][i],\n",
    "                       data['SWR34'][i], data['SER34'][i]])\n",
    "    else:\n",
    "        r34 = data['R34'][i]\n",
    "        \n",
    "    if np.isnan(data['R41'][i]):\n",
    "        r41 = np.max([data['NER41'][i], data['NWR41'][i],\n",
    "                       data['SWR41'][i], data['SER41'][i]])\n",
    "    else:\n",
    "        r41 = data['R41'][i]    \n",
    "    \n",
    "    r48 = data['R48'][i]\n",
    "    r64 = data['R64'][i]\n",
    "    \n",
    "    r = ma.array([r64, r48, r41, r34],\n",
    "                 mask=np.isnan([r64, r48, r41, r34]))\n",
    "    v = ma.array([64*1.4, 48*1.4, 41*1.4, 34*1.4],\n",
    "                 mask=np.isnan([r64, r48, r41, r34]))\n",
    "    \n",
    "    if np.isnan(rm) or np.isnan(pe) or np.isnan(pc) or np.isnan(vm):\n",
    "        beta[i] = beta[i]\n",
    "    else:\n",
    "        beta[i] = minimise(r, v, rm, vm, lati, loni, pe, pc)[0]\n",
    "\n",
    "\n",
    "plottimes(data, beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save the data to a formatted file, including the $\\beta$ values, for use in simulating the TC using TCRM. Note that you will have to manually edit the file to amend any ``nan`` values, either by removing the records containing them, or replacing them with estimated values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "track = np.core.records.fromarrays([np.ones(len(times)), timelist, age, lon, lat, speed, bearing, pcentre/100., penv/100., rmax, beta ],\n",
    "                 dtype=np.dtype({'names':(\"Index\", \"Datetime\", \"TimeElapsed\", \"Longitude\", \"Latitude\", \"Speed\", \"Bearing\", \n",
    "                                          \"CentralPressure\", \"EnvPressure\", \"Rmax\", \"beta\"),\n",
    "                                'formats':('i', 'object','f', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8', 'f8')}\n",
    "                                )\n",
    "                 ).T\n",
    "\n",
    "header = 'CycloneNumber,Datetime,TimeElapsed,Longitude,' + \\\n",
    "                 'Latitude,Speed,Bearing,' + \\\n",
    "                 'CentralPressure,EnvPressure,rMax,beta\\n'\n",
    "        \n",
    "fmt = '%i,%s,%7.3f,%8.3f,%8.3f,%6.2f,%6.2f,%7.2f,%7.2f,%6.2f,%6.3f'.split(',')\n",
    "with open(\"C:/temp/track.csv\", 'w') as fp:\n",
    "            fp.write('%' + header)\n",
    "            if len(track) > 0:\n",
    "                np.savetxt(fp, track, fmt=fmt,delimiter=',')        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we can interact with the data for a chosen time step. Select a time from the drop down list, then move the slider to visually find a fit to the data points. \n",
    "\n",
    "This highlights one issue with the analysis. There are some time steps where the radius to hurricane-force winds ($R_H$) is equal to the radius of maximum winds ($R_{max}$), but the maximum winds are greater than hurricane force. This is suspected to be because the radius to hurricane-force winds is only analysed when hurricane-force winds extend _around_ the cyclone, but the maximum winds are at a single location near the cyclone centre, and usually arises because of the asymmetry induced in the wind field due to the forward motion of the storm. \n",
    "\n",
    "To (partially) account for this, we presume that the forward speed of the storm contributes towards the maximum wind speed. Based on results from previous studies, we reduce the maximum wind speed ($V_m$) by $0.7 \\times V_{forward}$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def fitprofile(timestamp, beta=1.5):\n",
    "    \n",
    "    \n",
    "    idx = timelist.index(timestamp)\n",
    "    \n",
    "    pc = data['CentralPressure'][idx]\n",
    "    pe = data['PressureOCI'][idx]\n",
    "    rm = data['RadiusMaxWinds'][idx]\n",
    "    roci = data['RadiusOCI'][idx]\n",
    "    vm = data['WindSpeed'][idx] - 0.7*speed[idx]\n",
    "\n",
    "    if np.isnan(data['R34'][idx]):\n",
    "        r34 = np.mean([data['NER34'][idx], data['NWR34'][idx],\n",
    "                       data['SWR34'][idx], data['SER34'][idx]])\n",
    "    else:\n",
    "        r34 = data['R34'][idx]\n",
    "        \n",
    "    if np.isnan(data['R41'][idx]):\n",
    "        r41 = np.mean([data['NER41'][idx], data['NWR41'][idx],\n",
    "                       data['SWR41'][idx], data['SER41'][idx]])\n",
    "    else:\n",
    "        r41 = data['R41'][idx]    \n",
    "    \n",
    "    r48 = data['R48'][idx]\n",
    "    r64 = data['R64'][idx]\n",
    "    \n",
    "    r = np.arange(0, 201, 1., dtype=float)\n",
    "    h = holland(lat[idx], lon[idx], pe, pc, rm, beta)\n",
    "    p = powell(lat[idx], lon[idx], pe, pc, rm)\n",
    "    w = willoughby(lat[idx], lon[idx], pe, pc, rm)\n",
    "    try:\n",
    "        V = h.velocity(r)\n",
    "    except AssertionError:\n",
    "        print(\"Missing required fields\")\n",
    "    else:\n",
    "        fig = plt.figure(figsize=(16,12))\n",
    "        ax = plt.subplot(1,1,1)\n",
    "        ax.plot(r, abs(V), linewidth=2, color='b', label=\"Holland\")\n",
    "        Vp = p.velocity(r)\n",
    "        Vw = w.velocity(r)\n",
    "        ax.plot(r, abs(Vp), linewidth=1, color='r', label=\"Powell\")\n",
    "        ax.plot(r, abs(Vw), linewidth=1, color='g', label=\"Willoughby\")\n",
    "\n",
    "        ax.set_xlabel(\"Radial distance (km)\")\n",
    "        ax.set_ylabel(\"Wind speed (m/s)\")\n",
    "        ax.scatter(r34, convert(34*1.4, 'kts', 'mps'), s=30, marker='o')\n",
    "        ax.scatter(r41, convert(41*1.4, 'kts', 'mps'), s=30, marker='o')\n",
    "        ax.scatter(r48, convert(48*1.4, 'kts', 'mps'), s=30, marker='o')\n",
    "        ax.scatter(r64, convert(64*1.4, 'kts', 'mps'), s=30, marker='o')\n",
    "        ax.scatter(rm, vm - 0.7*speed[idx], s=30, marker='o', color='k')\n",
    "        ax.set_xlim((0,200))\n",
    "        ax.set_ylim((0, 80))\n",
    "        textstr = 'vmax = %.1f m/s \\nrmax = %.1f km\\nProfile vmax = %.1f m/s'%(vm,rm,abs(V).max())\n",
    "        props = dict(boxstyle='round', facecolor='gray', alpha=0.5)\n",
    "        ax.text(0.85, 0.95, textstr, transform=ax.transAxes, fontsize=12, verticalalignment='top', bbox=props)\n",
    "        ax.grid(True)\n",
    "        ax.legend(loc=7)\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a time that has all the required fields. If there aren't sufficient fields, then a little message will be printed, but doesn't indicate which fields are missing. You will need to look through the input file to check which fields are missing, or to identify a time that has the required fields. \n",
    "\n",
    "Use the dropdown to select a different time stamp to plot up. Once you have sufficient data to generate a plot, use the slider to adjust the beta value to visually get a good fit to the points. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "interact(fitprofile, \n",
    "         timestamp=widgets.DropdownWidget(values=timelist, value=timelist[0]), \n",
    "         beta=widgets.FloatSliderWidget(min=0.5, max=3.0, step=0.05, value=1.5, description=\"beta\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### References"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Holland, G. J. (1980): _An Analytic Model of the Wind and Pressure Profiles in Hurricanes_. Monthly Weather Review, __108__, 1212-1218\n",
    "2. Powell, M., G. Soukup, S. Cocke, S. Gulati, N. Morisseau-Leroy, S. Hamid, N. Dorst, & L. Axe (2005):  _State of Florida hurricane loss projection model: Atmospheric science component_. Journal of Wind Engineering and Industrial Aerodynamics, __93__, 651-674.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}