Analysis of permutations of vulnerability.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from __future__ import print_function, division\n",
    "import warnings\n",
    "warnings.simplefilter(action = \"ignore\", category = FutureWarning)\n",
    "import re\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as stats\n",
    "from os.path import join as pjoin\n",
    "sns.set_style('darkgrid')\n",
    "\n",
    "datapath = \"C:/WorkSpace/data/Derived/impact/cairns/bysuburb\"\n",
    "filename = \"cairns_aggimpact_classified_local_shield.csv\"\n",
    "\n",
    "data = pd.read_csv(pjoin(datapath, filename), header=0, skiprows=[1, 2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the full output file, and do some initial tidying up (i.e. remove a bunch of the extraneous metadata fields). We also take the chance to calculate the structural loss for each permutation of the vulnerability models. This bit can take some time (!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filename = \"cairns_impact_classified_local_shield.csv\"\n",
    "\n",
    "data = pd.read_csv(pjoin(datapath, filename), header=0)\n",
    "slr_cols = [col for col in data.keys() if col.startswith('structural_loss_ratio') and col[-1].isdigit()]\n",
    "slrdf = data[slr_cols]\n",
    "slrdf.loc[:, '0.2s gust at 10m height m/s'] = data['0.2s gust at 10m height m/s']\n",
    "slrdf.loc[:, 'SUBURB_2015'] = data['SUBURB_2015']\n",
    "slrdf.loc[:, 'LID'] = data['LID']\n",
    "slrdf.loc[:, 'REPLACEMENT_VALUE'] = data['REPLACEMENT_VALUE']\n",
    "slrdf.loc[:, 'M4'] = data['M4']\n",
    "slrdf.loc[:, 'YEAR_BUILT'] = data['YEAR_BUILT']\n",
    "sldf = slrdf[['0.2s gust at 10m height m/s', 'SUBURB_2015', 'LID', 'REPLACEMENT_VALUE', 'M4', 'YEAR_BUILT']]\n",
    "sl_cols = []\n",
    "for col in slr_cols:\n",
    "    slcolname = col.replace(\"_ratio\", \"\")\n",
    "    sl_cols.append(slcolname)\n",
    "    sldf.loc[:, slcolname] = slrdf[col] * sldf['REPLACEMENT_VALUE']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(16, 6))\n",
    "meansl = np.mean(sldf.loc[:, sl_cols], axis=1).values\n",
    "sumsl = np.sum(sldf.loc[:, sl_cols]).values\n",
    "meanslmean = np.mean(meansl)\n",
    "sumslmean = np.mean(sumsl)\n",
    "sns.distplot(meansl, bins=np.arange(0, 600000, 30000),\n",
    "             fit=stats.lognorm, axlabel=\"Mean structural loss ($)\", ax=ax0, \n",
    "             kde_kws={'label':'KDE', 'gridsize':10000, 'bw':10000},\n",
    "             hist_kws={'label':'Probability'},\n",
    "             fit_kws={'color':'red', 'label':'Fitted lognormal'})\n",
    "sns.distplot(sumsl, fit=stats.norm, axlabel=\"Total structural loss ($ million)\", ax=ax1, \n",
    "             kde_kws={'label':'KDE'},\n",
    "             hist_kws={'label':'Probability'},\n",
    "             fit_kws={'color':'red', 'label':'Fitted normal'})\n",
    "#ax0.axvline(meanslmean, label=\"Mean: ${0:0.0f}\".format(meanslmean))\n",
    "ax1.axvline(sumslmean, label=\"Mean: ${0:0.1f} million\".format(sumslmean/10**6))\n",
    "\n",
    "ax0.legend()\n",
    "ax0.set_xlim((0, 600000))\n",
    "ax0.set_yscale(\"log\", nonposy='clip')\n",
    "ax1.legend(loc=2)\n",
    "\n",
    "meanslcov = np.std(meansl)/np.mean(meansl)\n",
    "sumslcov = np.std(sumsl)/np.mean(sumsl)\n",
    "print(meanslcov, sumslcov)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The left panel is the mean distribution of structural losses across the assets (y-scale is logarithmic). The greatest proportion of losses are less than \\$10,000 (the mean is somewhere around \\$10,300). The red curve is a fitted lognormal distribution. I've tried a number of other distributions (exponential, weibull, pareto) and none fit as well as the lognormal.\n",
    "\n",
    "The right panel shows the distribution of total structural losses for the community. The red curve is a fitted normal distribution. \n",
    "\n",
    "In both cases, even though there's only 100 realisations of the losses, the distributions appear to be converging towards the respective fitted distributions. We can script this post-processing up and run a larger number of simulations (e.g. $10^5$) to confirm this convergence. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "pal = sns.color_palette(\"cubehelix\",5)\n",
    "sns.pairplot(data, vars=[\"0.2s gust at 10m height m/s\", \n",
    "                         \"structural_loss_ratio\", \n",
    "                         \"structural_loss\", \n",
    "                         'M4', \"REPLACEMENT_VALUE\"],\n",
    "             diag_kind='kde', diag_kws=dict(shade=True),  \n",
    "             hue='YEAR_BUILT', palette=pal, size=2.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A scatter plot of the various continuous variables in the output (I've ignored contents value, since our wind vulneraility models don't consider contents). In the first instance, they're colour-coded based on the AS4055 site classification (i.e. slope,  terrain and shielding factors), which are linked to the address, not the building atttributes. \n",
    "\n",
    "Highest gust speeds are recorded in locations where M4 values are around 1.1.\n",
    "\n",
    "Highest structural loss ratios occur at slightly lower M4 values, but near the peak gust wind speeds. The highest gusts don't produce the highest structural loss (or loss ratio). \n",
    "\n",
    "The loss vs loss ratio plot is curiously linear. I suspect that most of the data points are close to (0, 0), and we are seeing the extrema."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 3, figsize=(16, 18), sharey=True)\n",
    "ax = sns.stripplot(x='structural_loss_ratio', y='SUBURB_2015', data=data, jitter=True, ax=axes[0])\n",
    "ax = sns.stripplot(x='structural_loss', y='SUBURB_2015', data=data, jitter=True, ax=axes[1])\n",
    "ax = sns.stripplot(x='0.2s gust at 10m height m/s', y='SUBURB_2015', data=data, jitter=True, ax=axes[2])\n",
    "\n",
    "axes[0].set_xlim((0,1))\n",
    "axes[1].set_xlim((0, 600000))\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "slr = slrdf[slr_cols].values\n",
    "slrm = np.mean(slr, axis=1)\n",
    "sl = sldf[sl_cols].values\n",
    "slm = np.mean(sl, axis=1)\n",
    "fig, (ax0, ax1) = plt.subplots(2, 1, figsize=(12,16), sharex=True)\n",
    "sns.regplot(np.where(slrdf['M4']==-9999, 1, slrdf['M4']), slrm, ax=ax0)\n",
    "sns.regplot(sldf['M4'], slm, ax=ax1)\n",
    "#plt.scatter(slrm, slrdf['M4'])\n",
    "ax0.set_ylim((0, 1))\n",
    "ax1.set_xlim((0.5, 1.5))\n",
    "ax1.set_ylim((0, 600000))\n",
    "ax0.set_ylabel(\"Mean structural loss ratio\")\n",
    "ax1.set_ylabel(\"Mean structural loss ($)\")\n",
    "ax1.set_xlabel(\"M4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2,figsize=(16,8))\n",
    "ax = sns.distplot(slm, fit=stats.lognorm, label='Mean loss per asset', kde_kws=dict(bw=5000, label=\"KDE\"),\n",
    "                  fit_kws={'color':'red', 'label':'Lognorm fit'}, ax=axes[0])\n",
    "ax.set_xlim((0, 600000))\n",
    "ax.set_yscale(\"log\", nonposy='clip')\n",
    "ax.set_xlabel(\"Mean structural loss ($)\")\n",
    "ax.legend()\n",
    "\n",
    "ax = sns.distplot(slrm, fit=stats.lognorm, label='Mean loss ratio', kde_kws=dict(bw=0.025, label=\"KDE\"),\n",
    "                  fit_kws={'color':'red', 'label':'Lognorm fit'}, ax=axes[1])\n",
    "ax.set_xlim((0, 1))\n",
    "ax.set_yscale(\"log\", nonposy='clip')\n",
    "ax.set_xlabel(\"Mean loss ratio\")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "suburb_sl = sldf.groupby('SUBURB_2015')\n",
    "suburb_sl.agg([sum, np.mean, np.std, np.max])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.jointplot('structural_loss_ratio', 'structural_loss', data, kind='kde', xlim=(0,0.2), ylim=(0,50000) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pd.pivot_table(data, index='SUBURB_2015', columns='AS4055_CLASS', values='structural_loss_ratio', aggfunc=np.max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.distplot(data['structural_loss_ratio'])"
   ]
  },
  {
   "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
}