Fitting a GPD to TC wind observations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit Generalised Pareto Distribution to TC wind observations\n",
    "\n",
    "To verify the output from TCRM, we need to compare the synthetically-generated return period wind speeds with observations. Here, we extract maximum wind gust observations from Bureau of Meteorology (BoM) stations that are coincident with the passage of a tropical cyclone. \n",
    "\n",
    "We use a distance threshold of 600 km to select the observations, which is comparable to the 5$^{\\circ}$ margin around the synthetic TCs for which wind speeds are calculated within TCRM.\n",
    "\n",
    "Note that TCRM uses a Generalised Extreme Value (GEV) distribution fitted to annual maximum wind speeds, whereas we use a GPD in this analysis, as we have daily data. \n",
    "\n",
    "The observed data includes multiple observations for single TC events -- this may violate the requirement for independent observations. As a first pass, no attempt is made to filter these observations.\n",
    "\n",
    "### The Generalised Pareto Distribution\n",
    "The Generalised Pareto Distribution (GPD) is defined as:\n",
    "\n",
    "$H(y) = 1 - (1 + \\xi y / \\check{s}) ^{-1/\\xi}$\n",
    "\n",
    "where $\\check{s} = \\sigma + \\xi(u - \\mu)$ and $u = $ threshold value. $\\sigma$ and $\\mu$ are the scale and location parameters of a corresponding GEV distribution. If the data can be fitted to a GEV distribution, then values above the threshold can be fitted with a GPD.\n",
    "\n",
    "The fitted parameters are highly sensitive to the choice of threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import division\n",
    "import seaborn as sns\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from functools import wraps\n",
    "import time\n",
    "import os\n",
    "from os.path import join as pjoin\n",
    "from glob import glob\n",
    "\n",
    "from datetime import datetime\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from scipy.stats import genpareto\n",
    "import statsmodels.api as sm\n",
    "from statsmodels.nonparametric.kde import KDEUnivariate\n",
    "from statsmodels.distributions.empirical_distribution import ECDF\n",
    "\n",
    "import Utilities.metutils as metutils\n",
    "\n",
    "from return_period import returnLevels, empiricalReturnPeriod, returnPeriodUncertainty\n",
    "from distributions import fittedPDF\n",
    "\n",
    "from IPython.html.widgets import interact, fixed\n",
    "from IPython.html import widgets\n",
    "\n",
    "sns.set_context(\"poster\")\n",
    "sns.set_style(\"ticks\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def parseTime(yr, month, day, hour, minute):\n",
    "    \"\"\"\n",
    "    Parse year, month and day as strings and return a datetime.\n",
    "    \n",
    "    Handles the case of a missing time string (Pandas returns nan \n",
    "    if the field is empty).\n",
    "    \"\"\"\n",
    "    timestr = '{0}-{1:02d}-{2:02d} {3:02d}:{4:02d}'.format(yr, int(month), int(day), int(hour), int(minute))\n",
    "    \n",
    "    return datetime.strptime(timestr, '%Y-%m-%d %H:%M')\n",
    "\n",
    "def timer(func):\n",
    "    \"\"\"\n",
    "    Decorator to report execution time of a function/script.\n",
    "    \"\"\"\n",
    "    @wraps(func)\n",
    "    def wrap(*args, **kwargs):\n",
    "        t1 = time.time()\n",
    "        res = func(*args, **kwargs)\n",
    "\n",
    "        tottime = time.time() - t1\n",
    "        msg = \"%02d:%02d:%02d \" % \\\n",
    "          reduce(lambda ll, b : divmod(ll[0], b) + ll[1:],\n",
    "                        [(tottime,), 60, 60])\n",
    "\n",
    "        print \"Time for {0}: {1}\".format(func.func_name, msg) \n",
    "        return res\n",
    "\n",
    "    return wrap\n",
    "\n",
    "def find_nearest_index(array, value):\n",
    "    idx = (np.abs(array - value)).argmin()\n",
    "    return idx\n",
    "\n",
    "def calculateShape(mu, data):\n",
    "    \"\"\"\n",
    "    :param float mu: threshold parameter for the GPD distribution.\n",
    "    :param data: :class:`numpy.ndarray` of data values to fit.\n",
    "    \"\"\"\n",
    "    nobs = len(data)\n",
    "    nexc = len(data[data > mu])\n",
    "    rate = float(nexc)/float(nobs)\n",
    "    gpd = genpareto.fit(data[data > mu] - mu)\n",
    "\n",
    "    return gpd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plotDiagnostics(data, mu, xi, sigma):\n",
    "    \"\"\"\n",
    "    Create a 4-panel diagnostics plot of the fitted distribution.\n",
    "\n",
    "    :param data: :class:`numpy.ndarray` of observed data values (in units\n",
    "                 of metres/second).\n",
    "    :param float mu: Selected threshold value.\n",
    "    :param float xi: Fitted shape parameter.\n",
    "    :param float sigma: Fitted scale parameter.\n",
    "\n",
    "    \"\"\"\n",
    "    fig, ax = plt.subplots(2, 2)\n",
    "    axes = ax.flatten()\n",
    "    # Probability plots\n",
    "    sortedmax = np.sort(data[data > mu])   \n",
    "    gpdf = fittedPDF(data, mu, xi, sigma)\n",
    "    pp_x = sm.ProbPlot(sortedmax)\n",
    "    pp_x.ppplot(xlabel=\"Empirical\", ylabel=\"Model\", ax=axes[0], line='45')\n",
    "    axes[0].set_title(\"Probability plot\")\n",
    "\n",
    "    prplot = sm.ProbPlot(sortedmax, genpareto, distargs=(xi,),\n",
    "                         loc=mu, scale=sigma)\n",
    "    prplot.qqplot(xlabel=\"Model\", ylabel=\"Empirical\", ax=axes[1], line='45')\n",
    "    axes[1].set_title(\"Quantile plot\")\n",
    "\n",
    "    ax2 = axes[2]\n",
    "    rp = np.array([1, 2, 5, 10, 20, 50, 100, 200,\n",
    "                   500, 1000, 2000, 5000, 10000])\n",
    "    rate = float(len(sortedmax)) / float(len(data))\n",
    "    rval = returnLevels(rp, mu, xi, sigma, rate)\n",
    "\n",
    "    emprp = empiricalReturnPeriod(np.sort(data))\n",
    "    ax2.semilogx(rp, rval, label=\"Fitted RP curve\", color='r')\n",
    "    ax2.scatter(emprp[emprp > 1], np.sort(data)[emprp > 1],\n",
    "                color='b', label=\"Empirical RP\", s=100)\n",
    "    ax2.legend(loc=2)\n",
    "    ax2.set_xlabel(\"Return period\")\n",
    "    ax2.set_ylabel(\"Return level\")\n",
    "    ax2.set_title(\"Return level plot\")\n",
    "    ax2.grid(True)\n",
    "    maxbin = 4 * np.ceil(np.floor(data.max() / 4) + 1)\n",
    "    sns.distplot(sortedmax, bins=np.arange(mu, maxbin, 2),\n",
    "                 hist=True, axlabel='Wind speed (m/s)',\n",
    "                 ax=axes[3], kde_kws={\"label\":\"Empirical PDF\"})\n",
    "    axes[3].plot(sortedmax, gpdf, color='r', label=\"Fitted PDF\")\n",
    "    axes[3].set_title(\"Density plot\")\n",
    "    axes[3].legend(loc=1)\n",
    "    plt.tight_layout()\n",
    "    \n",
    "def plotFit(data, mu, xi, sigma, title):\n",
    "    \"\"\"\n",
    "    Plot a fitted distribution, with approximate 90% confidence interval\n",
    "    and empirical return period values.\n",
    "\n",
    "    :param data: :class:`numpy.ndarray` of observed data values.\n",
    "    :param float mu: Selected threshold value.\n",
    "    :param float xi: Fitted shape parameter.\n",
    "    :param float sigma: Fitted scale parameter.\n",
    "    :param str title: Title string for the plot.\n",
    "    :param str figfile: Path to store the file (includes image format)\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    rp = np.array([1, 2, 5, 10, 20, 50, 100, 200,\n",
    "                   500, 1000, 2000, 5000, 10000])\n",
    "    rate = float(len(data[data > mu])) / float(len(data))\n",
    "    rval = returnLevels(rp, mu, xi, sigma, rate)\n",
    "\n",
    "    emprp = empiricalReturnPeriod(data)\n",
    "    err = returnPeriodUncertainty(data, mu, xi, sigma, rp)\n",
    "\n",
    "    sortedmax = np.sort(data)\n",
    "    fig, ax1 = plt.subplots(1, 1)\n",
    "    ax1.semilogx(rp, rval, label=\"Fitted RP curve\")\n",
    "    ax1.semilogx(rp, rval + 1.96 * err, label=\"95% CI\",\n",
    "                 linestyle='--', color='0.5')\n",
    "    ax1.semilogx(rp, rval - 1.96 * err, linestyle='--', color='0.5')\n",
    "\n",
    "    ax1.scatter(emprp[emprp > 1], sortedmax[emprp > 1], s=100,\n",
    "                color='r', label=\"Empirical RP\")\n",
    "\n",
    "    title_str = (title + \"\\n\" +\n",
    "                 r\"$\\mu$ = {0:.3f}, $\\xi$ = {1:.5f}, $\\sigma$ = {2:.4f}\".\n",
    "                 format(mu, xi, sigma))\n",
    "    ax1.set_title(title_str)\n",
    "    ax1.legend(loc=2)\n",
    "    ax1.set_ylim((0, 100))\n",
    "    ax1.set_xlim((1, 10000))\n",
    "    ax1.set_ylabel('Wind speed (m/s)')\n",
    "    ax1.set_xlabel('Return period (years)')\n",
    "    ax1.grid(which='major')\n",
    "    ax1.grid(which='minor', linestyle='--', linewidth=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Data format\n",
    "\n",
    "Here we define the layout of the data within the files to enable easy reading via [`numpy.genfromtxt`](http://docs.scipy.org/doc/numpy-1.10.0/reference/generated/numpy.genfromtxt.html).\n",
    "In these files, wind speeds are recorded in km/h, so we convert to m/s for convenience.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "DTYPE = [('id', 'S8'), ('hm', 'S2'), ('StnNum', 'i'), ('Year', 'i'), ('Month', 'i'), \n",
    "         ('Day', 'i'), ('Hour', 'i'), ('Minute', 'i'), ('dtStdYear', 'i'), ('dtStdMonth', 'i'), \n",
    "         ('dtStdDay', 'i'), ('dtStdHour', 'i'), ('dtStdMinute', 'i'), ('Speed', 'f8'), \n",
    "         ('QSpeed', 'S1'), ('Dir', 'f8'), ('QDir', 'S1'), ('Gust', 'f8'), ('QGust', 'S1'), ('AWSFlag', 'S2'),\n",
    "         ('end', 'S1'), ('TCName', 'S10')]\n",
    "NAMES = [fields[0] for fields in DTYPE]\n",
    "CONVERT = {'Speed': lambda s: metutils.convert(float(s or 0), 'kmh', 'mps'),\n",
    "           'Gust': lambda s: metutils.convert(float(s or 0), 'kmh', 'mps')}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "@timer\n",
    "def selectThreshold(data, start=None):\n",
    "    \"\"\"\n",
    "    Select an appropriate threshold for fitting a generalised pareto\n",
    "    distribution. \n",
    "    \n",
    "    The only constraint placed on the selection is that the shape \n",
    "    parameter is negative (such that the distribution is bounded).\n",
    "    \n",
    "    :param data: :class:`numpy.ndarray` containing the observed values (with \n",
    "                 missing values removed).\n",
    "    :param start: starting point for the threshold value. If not given, \n",
    "                  defaults to the median of the ``data`` variable.\n",
    "    :returns: tuple of the shape, scale and threshold.\n",
    "    \"\"\"\n",
    "    \n",
    "    sh = []\n",
    "    sc = []\n",
    "    t = []\n",
    "    q1000list = []\n",
    "    q10000list = []\n",
    "    \n",
    "    eps = -0.01\n",
    "    datamax = data.max()\n",
    "    nobs = len(data)\n",
    "    if start:\n",
    "        startValue = start\n",
    "    else:\n",
    "        startValue = np.median(data)\n",
    "    for mu in np.arange(startValue, datamax, 0.01):\n",
    "        nexc = len(data[data > mu]) \n",
    "        rate = nexc / nobs\n",
    "        if nexc < 5:\n",
    "            break\n",
    "\n",
    "        pp = calculateShape(mu, data)\n",
    "        q1000, q10000 = returnLevels(np.array([1000, 10000]), mu, pp[0], pp[2], rate)\n",
    "        if np.isnan(q1000):\n",
    "            continue\n",
    "\n",
    "        if np.isnan(q10000):\n",
    "            continue\n",
    "\n",
    "        qdiff = np.abs(q10000 - q1000)\n",
    "        if pp[0] < eps and qdiff < 0.2*q10000 and qdiff > -eps: \n",
    "            t.append(mu)\n",
    "            sh.append(pp[0])\n",
    "            sc.append(pp[2])\n",
    "            q1000list.append(q1000)\n",
    "            q10000list.append(q10000)\n",
    "            \n",
    "    if len(t) == 0:\n",
    "        #print \"No suitable shape parameters identified\"\n",
    "        return 0, 0, 0\n",
    "    Av1000 = np.mean(np.array(q1000list))\n",
    "    Av10000 = np.mean(np.array(q10000list))\n",
    "    Av1000 = np.ceil(Av1000 + 0.05*Av1000)\n",
    "    Av10000 = np.ceil(Av10000 + 0.05*Av10000)\n",
    "\n",
    "    idx1000 = find_nearest_index(np.array(q1000list), Av1000)\n",
    "    idx10000 = find_nearest_index(np.array(q10000list), Av10000)\n",
    "    \n",
    "    u1000 = t[idx1000]\n",
    "    u10000 = t[idx10000]\n",
    "\n",
    "    if u1000 > u10000:\n",
    "        shmax = sh[idx1000]\n",
    "        scmax = sc[idx1000]\n",
    "    else:\n",
    "        shmax = sh[idx10000]\n",
    "        scmax = sc[idx10000]\n",
    "\n",
    "    return shmax, scmax, u1000    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "input_path = 'N:\\\\climate_change\\\\CHARS\\\\B_Wind\\\\data\\\\derived\\\\obs\\\\tc\\\\ibtracs'\n",
    "#input_path = \"C:\\\\WorkSpace\\\\data\\\\Derived\\\\obs\\\\tc\\\\ibtracs\\\\\"\n",
    "basename = 'bom_{0:06d}.csv'\n",
    "stnNum = \"4032\"\n",
    "stnName = \"Port Hedland\"\n",
    "\n",
    "fname = pjoin(input_path, basename.format(int(stnNum)))\n",
    "if os.path.exists(fname):\n",
    "    df = pd.read_csv(fname, skipinitialspace=True, skiprows=1, names=NAMES, \n",
    "                     parse_dates=[['Year', 'Month', 'Day', 'Hour', 'Minute']], \n",
    "                     date_parser=parseTime, index_col=False, converters=CONVERT)\n",
    "    df.describe()\n",
    "else:\n",
    "    print \"{0} does not exist\".format(fname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "plt.title('TC-related gust wind speeds for {0}'.format(stnName))\n",
    "plt.xlabel('Year')\n",
    "plt.ylabel('Gust wind speed (m/s)')\n",
    "x = [idx for idx in df.Year_Month_Day_Hour_Minute]\n",
    "y = df.Gust\n",
    "plt.scatter(x,y)\n",
    "plt.axhline(np.median(y), linewidth=1, linestyle='--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Filter out any suspect observations using the quality flag in the data. Also remove any null observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "quality = df['QGust'].fillna(\"X\").map(lambda x: x in ['Y','N','X',' ', np.nan])\n",
    "dmax = df['Gust'][df['Gust'].notnull() & quality]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we perform the fitting to a GPD. The function returns the location, shape and scale parameters of the GPD, which we can use to plot the fitted GPD against the empirical data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "xi, sigma, mu = selectThreshold(dmax, start=np.min(dmax))\n",
    "print xi, sigma, mu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a dummy vector that contains the TC-related observations, plus a bunch of zero values to extend the vector to the number of days between the start and end of the dataset. This gives us the chance to plot the empirical return periods along with the fitted return periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tdelta = df.Year_Month_Day_Hour_Minute.max().year - df.Year_Month_Day_Hour_Minute.min().year\n",
    "dummydata = np.zeros(int(tdelta+1)*365)\n",
    "ndata = len(dmax)\n",
    "dummydata[-ndata:] = dmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plotFit(dummydata, mu, xi, sigma, stnName)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plotDiagnostics(dummydata, mu, xi, sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "@timer\n",
    "def selThreshold(dmax):\n",
    "    \"\"\"\n",
    "    Select the best fitting threshold that maximises the return period values, but minimises the $R^$ value\n",
    "    when fitted against the observed distribution.\n",
    "    \"\"\"\n",
    "    eps = -0.01\n",
    "    datamax = data.max()\n",
    "    nobs = len(data)\n",
    "    mu = np.median(data)\n",
    "    while mu < datamax:\n",
    "        nexc = len(data[data > mu])\n",
    "        exceed = data[data > mu]\n",
    "        rate = nexc / nobs\n",
    "        if nexc < 10:\n",
    "            break\n",
    "        pp = calculateShape(mu, data)\n",
    "        \n",
    "        if pp[0] > eps:\n",
    "            break\n",
    "            \n",
    "        emppdf = empiricalPDF(exceed)\n",
    "        \n",
    "        try:\n",
    "            popt, pcov = curve_fit(lambda x, xi, sigma: \\\n",
    "                                   genpareto.pdf(x, xi, loc=mu, scale=sigma),\n",
    "                                   np.sort(exceed), emppdf, (xi, sigma))\n",
    "        except RuntimeError as e:\n",
    "            return 0.\n",
    "        sd = np.sqrt(np.diag(pcov))"
   ]
  }
 ],
 "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
}