Covariance and correlation testing of TC track behaviour.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Covariance and correlation testing of TC track behaviour\n",
    "\n",
    "To demonstrate the autoregresive correlation and covariance of various track parameters, we extract those parameters from the track data and plot up the covariance and correlation matrices. \n",
    "\n",
    "The intention here is to address reviewer concerns about the use of a lag-1 autoregressive model for the intensity of TC tracks. We examine the covariance and correlation of the parameters for the first 7 days of the life of a TC. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from __future__ import print_function, division\n",
    "\n",
    "import io\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "from os.path import join as pjoin\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import ImageGrid\n",
    "\n",
    "from scipy.stats import pearsonr, betai\n",
    "\n",
    "\n",
    "from Utilities.config import ConfigParser\n",
    "from Utilities.loadData import loadTrackFile\n",
    "from Utilities.config import ConfigParser\n",
    "from Utilities.track import Track\n",
    "from Utilities.metutils import convert\n",
    "\n",
    "import seaborn as sns\n",
    "sns.set_context(\"notebook\")\n",
    "sns.set_style(\"whitegrid\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "configstr = \"\"\"\n",
    "[DataProcess]\n",
    "InputFile=C:/WorkSpace/data/TC/Allstorms.ibtracs_wmo.v03r06.csv\n",
    "Source=IBTRACS\n",
    "StartSeason=1961\n",
    "FilterSeasons=True\n",
    "\n",
    "[Region]\n",
    "; Domain for windfield and hazard calculation\n",
    "gridLimit={'xMin':90.,'xMax':180.,'yMin':-30.0,'yMax':-5.0}\n",
    "gridSpace={'x':1.0,'y':1.0}\n",
    "gridInc={'x':1.0,'y':0.5}\n",
    "\n",
    "[TrackGenerator]\n",
    "NumSimulations=5000\n",
    "YearsPerSimulation=10\n",
    "SeasonSeed=68876543\n",
    "TrackSeed=334825\n",
    "TimeStep=1.0\n",
    "\n",
    "[Input]\n",
    "landmask = C:/WorkSpace/tcrm/input/landmask.nc\n",
    "mslpfile = C:/WorkSpace/data/MSLP/slp.day.ltm.nc\n",
    "datasets = IBTRACS,LTMSLP\n",
    "\n",
    "[Output]\n",
    "Path=C:/WorkSpace/data/TC/aus\n",
    "\n",
    "[Hazard]\n",
    "Years=2,5,10,20,25,50,100,200,250,500,1000,2000,2500,5000\n",
    "MinimumRecords=10\n",
    "CalculateCI=False\n",
    "\n",
    "[Logging]\n",
    "LogFile=C:/WorkSpace/data/TC/aus/log/aus.log\n",
    "LogLevel=INFO\n",
    "Verbose=False\n",
    "\n",
    "[IBTRACS]\n",
    "; Input data file settings\n",
    "url = ftp://eclipse.ncdc.noaa.gov/pub/ibtracs/v03r06/wmo/csv/Allstorms.ibtracs_wmo.v03r06.csv.gz\n",
    "path = C:/WorkSpace/data/TC/\n",
    "filename = Allstorms.ibtracs_wmo.v03r06.csv\n",
    "columns = tcserialno,season,num,skip,skip,skip,date,skip,lat,lon,skip,pressure\n",
    "fielddelimiter = ,\n",
    "numberofheadinglines = 3\n",
    "pressureunits = hPa\n",
    "lengthunits = km\n",
    "dateformat = %Y-%m-%d %H:%M:%S\n",
    "speedunits = kph\n",
    "\n",
    "[LTMSLP]\n",
    "; MSLP climatology file settings\n",
    "URL = ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.derived/surface/slp.day.1981-2010.ltm.nc\n",
    "path = C:/WorkSpace/data/MSLP\n",
    "filename = slp.day.ltm.nc\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "config = ConfigParser()\n",
    "config.readfp(io.BytesIO(configstr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the track file and pull out those tracks that are within the specified domain. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "trackFile = config.get('DataProcess', 'InputFile')\n",
    "source = config.get('DataProcess', 'Source')\n",
    "\n",
    "print(\"Track file: {0}\".format(trackFile))\n",
    "print(\"Track format: {0}\".format(source))\n",
    "alltracks = loadTrackFile(configstr, trackFile, source)\n",
    "#tracks = parseTracks(configstr, trackFile, source, 1.)\n",
    "\n",
    "gridLimit=config.geteval('Region', 'gridLimit')\n",
    "\n",
    "tracks = [track for track in alltracks if track.inRegion(gridLimit)]\n",
    "\n",
    "print(\"There are {0:d} tracks in the input dataset\".format(len(tracks)))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to find the indices of all 6-hourly timesteps (but not 3-hourly, etc.). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def getIndices(elapsedTime, delta=6, maxtime=120):\n",
    "    filter1 = (np.mod(elapsedTime, delta)==0)\n",
    "    filter2 = (elapsedTime<=maxtime)\n",
    "    idx = np.nonzero(filter1 & filter2)\n",
    "    return idx\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def corrcoef(matrix):\n",
    "    rows, cols = matrix.shape[0], matrix.shape[1]\n",
    "    r = np.ones(shape=(rows, rows))\n",
    "    p = np.ones(shape=(rows, rows))\n",
    "    for i in range(rows):\n",
    "        for j in range(i+1, rows):\n",
    "            r_, p_ = pearsonr(matrix[i], matrix[j])\n",
    "            r[i, j] = r[j, i] = r_\n",
    "            p[i, j] = p[j, i] = p_\n",
    "    return r, p"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we parse each track and pull out the parameters of interest (pressure, speed, bearing, pressure rate of change) at each 6-hourly timestep. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ntracks = len(tracks)\n",
    "ndays = 7\n",
    "nhours = ndays * 24.\n",
    "interval = 6\n",
    "steps = range(interval, int(nhours + 1), interval)\n",
    "prsteps = range(0, int(nhours + 1), interval)\n",
    "alldp = ma.zeros((nhours/interval, ntracks))\n",
    "allprs = ma.zeros((nhours/interval + 1, ntracks))\n",
    "allspd = ma.zeros((nhours/interval + 1, ntracks))\n",
    "allbear = ma.zeros((nhours/interval + 1, ntracks))\n",
    "for i, track in enumerate(tracks):\n",
    "    idx = getIndices(track.TimeElapsed, delta=interval, maxtime=nhours)\n",
    "    tstep = track.TimeElapsed[idx]\n",
    "    tstepidx = [steps.index(int(t)) for t in tstep if int(t)!=0]\n",
    "    prtstepidx = [prsteps.index(int(t)) for t in tstep]\n",
    "\n",
    "\n",
    "    dt = np.diff(track.TimeElapsed[idx]) \n",
    "    dp = np.diff(track.CentralPressure[idx])\n",
    "    prs = track.CentralPressure[idx]\n",
    "    spd = track.Speed[idx]\n",
    "    bear = track.Bearing[idx]\n",
    "    pmask = (track.CentralPressure[idx[0][1:]]!=0)\n",
    "    prs = ma.array(prs, mask=(track.CentralPressure[idx[0]]!=0))\n",
    "    spd = ma.array(spd)\n",
    "    bear = ma.array(bear)\n",
    "    dpdt = dp/interval\n",
    "    dpmask = (np.abs(dpdt) > 10)\n",
    "\n",
    "    dpdt = ma.array(dp/interval, mask=(pmask & dpmask))\n",
    "\n",
    "    if len(dpdt)!=len(tstepidx):\n",
    "        pass\n",
    "    else:\n",
    "        alldp[tstepidx, i] = dpdt\n",
    "    allprs[prtstepidx, i] = prs\n",
    "    allspd[prtstepidx, i] = spd\n",
    "    allbear[prtstepidx, i] = bear\n",
    "    \n",
    "alldp[np.where(np.abs(alldp) > 10)] = ma.masked\n",
    "\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "cm = ax.pcolor(alldp, cmap='RdBu', vmin=-5, vmax=5)\n",
    "plt.colorbar(cm)\n",
    "\n",
    "ax.set_yticklabels(range(0,int(nhours)+1, 24))\n",
    "\n",
    "corr, dpp = corrcoef(alldp)\n",
    "cov = np.cov(alldp)\n",
    "pcor, pp = corrcoef(allprs)\n",
    "pcov = np.cov(allprs)\n",
    "scor, sp = corrcoef(allspd)\n",
    "scov = np.cov(allspd)\n",
    "bcor, bp = corrcoef(allbear)\n",
    "bcov = np.cov(allbear)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(1, (18, 6))\n",
    "grid = ImageGrid(fig, 111, nrows_ncols=(1,3), axes_pad=0.2,\n",
    "                 add_all=True,cbar_location=\"top\", cbar_mode=\"each\",\n",
    "                 cbar_size=\"5%\", cbar_pad=\"2%\")\n",
    "im0 = grid[0].pcolor(np.arange(0, nhours, interval), np.arange(0, nhours, interval), cov, cmap='gray_r',)\n",
    "grid[0].cax.colorbar(im0)\n",
    "grid[0].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].set_yticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].autoscale()\n",
    "\n",
    "im1 = grid[1].pcolor(np.arange(0, nhours, interval), np.arange(0, nhours, interval), corr, cmap='seismic', vmin=-1, vmax=1)\n",
    "grid[1].cax.colorbar(im1)\n",
    "grid[1].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[1].autoscale()\n",
    "\n",
    "im2 = grid[2].pcolor(np.arange(0, nhours, interval), np.arange(0, nhours, interval), dpp, cmap='gray_r', vmax=0.05)\n",
    "grid[2].cax.colorbar(im2)\n",
    "grid[2].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[2].autoscale()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On the left is the covariance matrix, centre the correlation matrix and right is the p-values for the Peason correlation for the pressure rate of change every six hours ($p > 0.05$ are censored). As expected, the diagonal elements are all unity, representing self-correlation. The off-diagonal elements represent auto-regressive correlations, with the first off-diagonal elements representing the lag-1 autocorrelation. There is generally near-zero correlation off the diagonal, except for $48 \\lt t \\lt 96$ hours. In this time range, there is negative lag-1 autocorrelation, which may be interpreted as a period when fluctuating intensity is more common. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(1, (18, 6))\n",
    "grid = ImageGrid(fig, 111, nrows_ncols=(1,3), axes_pad=0.2,\n",
    "                 add_all=True,cbar_location=\"top\", cbar_mode=\"each\",\n",
    "                 cbar_size=\"5%\", cbar_pad=\"2%\")\n",
    "im0 = grid[0].pcolor(np.arange(0, nhours + 1, interval), \n",
    "                     np.arange(0, nhours + 1, interval), \n",
    "                     pcov[1:, 1:], cmap='gray_r',)\n",
    "grid[0].cax.colorbar(im0)\n",
    "grid[0].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].set_yticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].autoscale()\n",
    "im1 = grid[1].pcolor(np.arange(0, nhours + 1, interval), \n",
    "                     np.arange(0, nhours + 1, interval), \n",
    "                     pcor[1:, 1:], cmap='seismic', vmin=-1, vmax=1)\n",
    "grid[1].cax.colorbar(im1)\n",
    "grid[1].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[1].autoscale()\n",
    "\n",
    "im2 = grid[2].pcolor(np.arange(0, nhours, interval), \n",
    "                     np.arange(0, nhours, interval), \n",
    "                     pp[1:, 1:], cmap='gray_r', vmax=0.05)\n",
    "grid[2].cax.colorbar(im2)\n",
    "grid[2].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[2].autoscale()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, the covariance and correlation matrices are derived from the pressure observations. There is clearly a very strong autocorrelation, with the off-diagonal terms all displaying strongly positive correlations. In fact, the strength of the autocorrelation increases for $t>48$ hours.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(1, (18, 6))\n",
    "grid = ImageGrid(fig, 111, nrows_ncols=(1,3), axes_pad=0.2,\n",
    "                 add_all=True,cbar_location=\"top\", cbar_mode=\"each\",\n",
    "                 cbar_size=\"5%\", cbar_pad=\"2%\")\n",
    "im0 = grid[0].pcolor(np.arange(0, nhours+1, interval), \n",
    "                     np.arange(0, nhours+1, interval), \n",
    "                     scov[1:, 1:], cmap='gray_r',)\n",
    "grid[0].cax.colorbar(im0)\n",
    "grid[0].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].set_yticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].autoscale()\n",
    "im1 = grid[1].pcolor(np.arange(0, nhours + 1, interval), \n",
    "                     np.arange(0, nhours + 1, interval), \n",
    "                     scor[1:, 1:], cmap='seismic', vmin=-1, vmax=1)\n",
    "grid[1].cax.colorbar(im1)\n",
    "grid[1].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[1].autoscale()\n",
    "\n",
    "im2 = grid[2].pcolor(np.arange(0, nhours+1, interval), \n",
    "                     np.arange(0, nhours+1, interval), \n",
    "                     sp[1:, 1:], cmap='gray_r', vmax=0.05)\n",
    "grid[2].cax.colorbar(im2)\n",
    "grid[2].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[2].autoscale()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(1, (18, 6))\n",
    "grid = ImageGrid(fig, 111, nrows_ncols=(1,3), axes_pad=0.2,\n",
    "                 add_all=True,cbar_location=\"top\", cbar_mode=\"each\",\n",
    "                 cbar_size=\"5%\", cbar_pad=\"2%\")\n",
    "im0 = grid[0].pcolor(np.arange(0, nhours + 1, interval), \n",
    "                     np.arange(0, nhours + 1, interval), bcov[1:, 1:], cmap='gray_r',)\n",
    "grid[0].cax.colorbar(im0)\n",
    "grid[0].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].set_yticks(np.arange(24, nhours + 1, 24))\n",
    "grid[0].autoscale()\n",
    "im1 = grid[1].pcolor(np.arange(0, nhours + 1, interval), \n",
    "                     np.arange(0, nhours + 1, interval), \n",
    "                     bcor[1:, 1:], cmap='seismic', vmin=-1, vmax=1)\n",
    "grid[1].cax.colorbar(im1)\n",
    "grid[1].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[1].autoscale()\n",
    "\n",
    "im2 = grid[2].pcolor(np.arange(0, nhours + 1, interval), \n",
    "                     np.arange(0, nhours + 1, interval), \n",
    "                     bp[1:, 1:], cmap='gray_r', vmax=0.05)\n",
    "grid[2].cax.colorbar(im2)\n",
    "grid[2].set_xticks(np.arange(24, nhours + 1, 24))\n",
    "grid[2].autoscale()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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
}