summary_stats_with_pandas.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Back to Index](index.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Summary Statistics for Sacramento County"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Begin by importing all python modules we will need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# A 'magic' command to display plots inline\n",
    "%matplotlib inline\n",
    "\n",
    "import requests\n",
    "import pandas as pd\n",
    "from datetime import datetime\n",
    "import itertools\n",
    "import calendar\n",
    "import seaborn as sns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Get projections\n",
    "\n",
    "The code in the next cell contains a bunch of functions to:\n",
    "    - fetch annual averages of Max. Temp, Min. Temp and Precipitation for 4 GCMs and 2 scenarios for Sacramento County from the API\n",
    "    - convert units\n",
    "    - return a new dataframe that contains all the data\n",
    "    \n",
    "Just run the cell, you can go over the code at your leisure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def kelvin_to_F(val):\n",
    "    return  (val - 273.15) * 9/5 + 32\n",
    "\n",
    "def kgm2s_to_inchyear(row):\n",
    "    year = row.event.year\n",
    "    days = 366 if calendar.isleap(year) else 365\n",
    "    return (row.image * 86400) * 0.0393701 * days\n",
    "\n",
    "def get_projections():\n",
    "    # Create an empty list to hold dataframes\n",
    "    df_list = []\n",
    "\n",
    "    # Make a combined list of GCMs, scenarios, climate variables for looping\n",
    "    climvar = ['tasmax', 'tasmin', 'pr']\n",
    "    period = ['year']\n",
    "    models = ['CanESM2', 'CNRM-CM5', 'HadGEM2-ES', 'MIROC5']\n",
    "    scenarios = ['rcp45', 'rcp85']\n",
    "    zipped = itertools.product(climvar, period, models, scenarios)\n",
    "\n",
    "    # Request parameters\n",
    "    params = {'pagesize': 94, 'stat': 'mean', 'ref': '/api/counties/34/'}\n",
    "    \n",
    "    # Loop through zipped\n",
    "    for item in zipped:\n",
    "        # Create slug\n",
    "        slug = '_'.join(item)\n",
    "        url = 'http://api.cal-adapt.org/api/series/%s/rasters/' % slug\n",
    "        # Make request\n",
    "        response = requests.get(url, params=params)\n",
    "\n",
    "        # Get data\n",
    "        if response.ok:\n",
    "            print('Processing:', slug)\n",
    "            data = response.json()\n",
    "            # Create temp dataframe\n",
    "            df = pd.DataFrame(data['results'])\n",
    "            # Change format of `event` field to datetime\n",
    "            df['event'] = pd.to_datetime(df['event'], format='%Y-%m-%d')\n",
    "            # The data from API should be sorted, but sort ascending to be sure\n",
    "            df = df.sort_values('event')  \n",
    "            # Convert units\n",
    "            if 'tas' in slug:\n",
    "                df.image = df.image.apply(lambda x: kelvin_to_F(x))\n",
    "            else:\n",
    "                df.image = df.apply(kgm2s_to_inchyear, axis=1)\n",
    "            # Discard all columns except event and image\n",
    "            df = df[['event', 'image']]\n",
    "            df['climvar'] = item[0]\n",
    "            df['model'] = item[2]\n",
    "            df['scenario'] = item[3]\n",
    "            df_list.append(df)          \n",
    "        else:\n",
    "            print('Failed:', slug)\n",
    "    # Combine all the dataframes into one and return\n",
    "    return pd.concat(df_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute the `get_projections()` function. If all goes well, you should see a list of messages as the code process each timeseries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "projections = get_projections()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore the `projections` dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "projections = projections.reset_index(drop=True)\n",
    "projections.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Get historical observed data\n",
    "\n",
    "The code in the next cell contains a bunch of functions to:\n",
    "    - fetch annual averages of Max. Temp, Min. Temp and Precipitation for the livneh data for Sacramento County from the API\n",
    "    - convert units\n",
    "    - return a new dataframe that contains all the data\n",
    "    \n",
    "Just run the cell, you can go over the code at your leisure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def celsius_to_F(val):\n",
    "    return val * 9/5 + 32 \n",
    "\n",
    "def mmday_to_inchyear(row):\n",
    "    year = row.event.year\n",
    "    days = 366 if calendar.isleap(year) else 365\n",
    "    return row.image * 0.0393701 * days\n",
    "\n",
    "def get_observed():\n",
    "    # Create an empty list to hold dataframes\n",
    "    df_list = []\n",
    "\n",
    "    # Make a combined list of GCMs, scenarios, climate variables for looping\n",
    "    climvar = ['tasmax', 'tasmin', 'pr']\n",
    "    period = ['year']\n",
    "    zipped = itertools.product(climvar, period, ['livneh'])\n",
    "    \n",
    "    # Request parameters\n",
    "    params = {'pagesize': 64, 'stat': 'mean', 'ref': '/api/counties/34/'}\n",
    "\n",
    "    # Loop through zipped\n",
    "    for item in zipped:\n",
    "        # Create slug\n",
    "        slug = '_'.join(item)\n",
    "        url = 'http://api.cal-adapt.org/api/series/%s/rasters/' % slug \n",
    "        # Make request\n",
    "        response = requests.get(url, params=params)\n",
    "\n",
    "        # Get data\n",
    "        if response.ok:\n",
    "            print('Processing:', slug)\n",
    "            data = response.json()\n",
    "            # Create temp dataframe\n",
    "            df = pd.DataFrame(data['results'])\n",
    "            # Change format of `event` field to datetime\n",
    "            df.event = pd.to_datetime(df.event, format='%Y-%m-%d')\n",
    "            # The data from API should be sorted, but sort ascending to be sure\n",
    "            df = df.sort_values('event')     \n",
    "            # Convert units\n",
    "            if 'tas' in slug:\n",
    "                df.image = df.image.apply(lambda x: celsius_to_F(x))\n",
    "            else:\n",
    "                df.image = df.apply(mmday_to_inchyear, axis=1)\n",
    "            col_list = ['event', 'image']\n",
    "            df = df[col_list]\n",
    "            df['climvar'] = item[0]\n",
    "            df['model'] = 'livneh'\n",
    "            df['scenario'] = 'historical'\n",
    "            df_list.append(df)          \n",
    "        else:\n",
    "            print('Failed:', slug)\n",
    "    # Combine all the dataframes into one and return\n",
    "    return pd.concat(df_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Execute the `get_observed()` function. If all goes well, you should see a list of messages as the code process each timeseries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "observed = get_observed()\n",
    "observed = observed.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Explore the observed dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "observed.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Summary stats\n",
    "\n",
    "Combine projections and observed data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.concat([projections, observed])\n",
    "df = df.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate summary stats for baseline period 1961-1990"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "baseline = df[(df.event >= '1961-01-01') & (df.event <= '1990-01-01')]\n",
    "baseline = baseline.groupby(['climvar', 'scenario', 'model'])['image'].agg(['mean', 'max', 'min', 'std'])\n",
    "# groupby return a groupy object, convert it to a dataframe\n",
    "df_baseline = pd.DataFrame(baseline)\n",
    "df_baseline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate summary stats for 2020-2050"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "projections_2020_2050 = df[(df.event >= '2020-01-01') & (df.event <= '2050-01-01')]\n",
    "projections_2020_2050 = projections_2020_2050.groupby(['climvar', 'scenario', 'model'])['image'].agg(['mean', 'max', 'min', 'std'])\n",
    "# Convert groupby object to dataframe\n",
    "df_2020_2050 = pd.DataFrame(projections_2020_2050)\n",
    "df_2020_2050"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "projections_2070_2099 = df[(df.event >= '2070-01-01') & (df.event <= '2099-01-01')]\n",
    "projections_2070_2099 = projections_2070_2099.groupby(['climvar', 'scenario', 'model'])['image'].agg(['mean', 'max', 'min', 'std'])\n",
    "# Convert groupby object to dataframe\n",
    "df_2070_2099 = pd.DataFrame(projections_2070_2099)\n",
    "df_2070_2099"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate difference from baseline values for Tasmax and get quantiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Get a tasmax data for 2070-2099\n",
    "tmax = df[(df.event >= '2070-01-01') & (df.event <= '2099-01-01') & (df.climvar == 'tasmax')]\n",
    "# Create a new dataframe since we want to modify values\n",
    "tmax = pd.DataFrame(tmax)\n",
    "# Subtract baseline period mean from all values\n",
    "tmax.image = tmax.image - df_baseline.loc[('tasmax','historical', 'livneh'),'mean']\n",
    "# Groupby scenario and calculate quantiles. This uses data from all 4 models.\n",
    "tmax_quantiles = tmax.groupby(['scenario'])['image'].quantile([0.1, 0.25, 0.75, 0.9])\n",
    "# Convert groupby object to dataframe if you want to export later\n",
    "tmax_quantiles = pd.DataFrame(tmax_quantiles)\n",
    "tmax_quantiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "### 4. Export to Excel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "writer = pd.ExcelWriter('output.xlsx')\n",
    "df_baseline.to_excel(writer,'Baseline')\n",
    "df_2020_2050.to_excel(writer,'2020-2050')\n",
    "df_2070_2099.to_excel(writer,'2070-2099')\n",
    "writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "tmax_quantiles.to_csv('tmax.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Make charts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a linechart showing all timeseries for Max. Temperature and RCP 4.5. Export it to a png file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tasmax_rcp45 = df.loc[(df['climvar'] == 'tasmax') & (df['scenario'] != 'rcp85')]\n",
    "tasmax_rcp45 = tasmax_rcp45.pivot(index='event', columns='model', values='image')\n",
    "plot = tasmax_rcp45.plot()\n",
    "fig = plot.get_figure()\n",
    "fig.savefig('output.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Categorical [Scatterplots](https://seaborn.pydata.org/tutorial/categorical.html#categorical-scatterplots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tasmax_rcp85 = df.loc[(df['climvar'] == 'tasmax') & (df['scenario'] != 'rcp45')]\n",
    "sns.stripplot(x=\"model\", y=\"image\", data=tasmax_rcp85, jitter=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tasmax = df.loc[df['climvar'] == 'tasmax']\n",
    "sns.swarmplot(y=\"model\", x=\"image\", hue=\"scenario\", data=tasmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "precip = df.loc[df['climvar'] == 'pr']\n",
    "sns.swarmplot(y=\"model\", x=\"image\", hue=\"scenario\", data=precip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.boxplot(y=\"model\", x=\"image\", hue=\"scenario\", data=tasmax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "tas = df.loc[(df['climvar'] == 'tasmax') | (df['climvar'] == 'tasmin')]\n",
    "sns.factorplot(y=\"model\", x=\"image\", hue=\"scenario\", kind=\"box\", col=\"climvar\", data=tas)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}