usgs_preprocessing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "import time\n",
    "\n",
    "\n",
    "def str_to_utc(s):\n",
    "    return int(time.mktime(dt.datetime.strptime(s, '%Y-%m-%d %H:%M:%S').utctimetuple()))\n",
    "\n",
    "def utc_to_str(t):\n",
    "    return dt.datetime.fromtimestamp(t).strftime('%Y-%m-%d %H:%M:%S')\n",
    "\n",
    "\n",
    "to_drop = ['agency_cd', 'site_no', 'tz_cd', 'GH_Inst_cd', 'year', 'month', 'week', 'yday', 'hour', 'byhour',\n",
    "           'byday', 'byweek', 'bymonth']\n",
    "\n",
    "data = pd.read_csv(\"gh.csv\")\n",
    "data['ts_tmp'] = data['dateTime'].apply(str_to_utc)\n",
    "data.drop(to_drop, axis=1, inplace=True)\n",
    "\n",
    "t = data['ts_tmp'].values\n",
    "idx = pd.DataFrame(t[0] + 900*np.arange((t[-1] - t[0]) // 900 + 1), columns=['timestamp'])\n",
    "D = idx.merge(data, how='left', left_on='timestamp', right_on='ts_tmp')\n",
    "D.drop('ts_tmp', axis=1, inplace=True)\n",
    "\n",
    "D.columns = ['timestamp', 'date', 'gh_raw']\n",
    "D['hour_id'] = D['timestamp'].values // 3600\n",
    "D['day_id']  = (D['timestamp'].values + 2*3600) // 86400\n",
    "D['week_id'] = (D['timestamp'].values - 86400*4 + 2*3600) // (86400*7)\\\n",
    "\n",
    "\n",
    "def find_skips(x):\n",
    "    skips = list()\n",
    "    start = -1\n",
    "    end = -1\n",
    "    for i in range(x.size):\n",
    "        if np.isnan(x[i]) and not np.isnan(x[i-1]):\n",
    "            start = i-1\n",
    "        if not np.isnan(x[i]) and np.isnan(x[i-1]):\n",
    "            end = i\n",
    "            skips.append([start, end])\n",
    "    return skips\n",
    "\n",
    "skips = find_skips(D['gh_raw'].values)\n",
    "\n",
    "\n",
    "def bezier_interpolator(x_a, y_a, dy_a, x_b, y_b, dy_b, x):\n",
    "    if x.size < 4:\n",
    "        By = 0.5*(y_a + y_b) * np.ones(x.size)\n",
    "    \n",
    "    elif dy_a * dy_b > 0:\n",
    "        # cubic bezier curve\n",
    "        x_1 = x[x.size // 3]  \n",
    "        x_2 = x[x.size - x.size // 3] \n",
    "        y_1 = y_a + dy_a*(x_1 - x_a)        \n",
    "        y_2 = y_b + dy_b*(x_2 - x_b)        \n",
    "        t = np.zeros(x.size)\n",
    "        for i in range(x.size):\n",
    "            r = np.roots([x_b - 3*x_2 + 3*x_1 - 1*x_a,\n",
    "                                3*x_2 - 6*x_1 - 3*x_a,\n",
    "                                        3*x_1 - 3*x_a,\n",
    "                                                1*x_a - x[i]])\n",
    "            idx = np.where(np.logical_and(r <= 1, r >= 0))[0][0]\n",
    "            t[i] = r[idx]\n",
    "        By = (1-t)**3 * y_a + 3*t*(1-t)**2 * y_1 + 3*(1-t)*t**2 * y_2 + t**3 * y_b\n",
    "        \n",
    "    else:\n",
    "        # quadratic bezier curve\n",
    "        x_1 = (y_b - y_a + dy_a*x_a - dy_b*x_b) / (dy_a - dy_b)\n",
    "        y_1 = y_a + dy_a*(x_1 - x_a)        \n",
    "        t = np.zeros(x.size)\n",
    "        for i in range(x.size):\n",
    "            r = np.roots([x_b - 2*x_1 - 1*x_a, \n",
    "                                2*x_1 - 2*x_a, \n",
    "                                        1*x_a - x[i]])\n",
    "            idx = np.where(np.logical_and(r <= 1, r >= 0))[0][0]\n",
    "            t[i] = r[idx]\n",
    "        By = (1-t)**2 * y_a + 2*t*(1-t)*y_1 + t**2 * y_b\n",
    "        \n",
    "    return By\n",
    "\n",
    "\n",
    "def impute(z):\n",
    "    y = z.copy()\n",
    "    dy = y[1:] - y[:-1]\n",
    "    for skip in skips:\n",
    "        a = skip[0]\n",
    "        x_a = 0\n",
    "        y_a = y[a]\n",
    "        dy_a = dy[a-16:a].mean()\n",
    "        b = skip[1]    \n",
    "        x_b = b - a\n",
    "        y_b = y[b]    \n",
    "        dy_b = dy[b+2:b+18].mean()\n",
    "        x = np.arange(b - a)\n",
    "        y[a:b] = np.round(100*bezier_interpolator(x_a, y_a, dy_a, x_b, y_b, dy_b, x)) / 100\n",
    "    return y\n",
    "  \n",
    "    \n",
    "D['gh'] = impute(D['gh_raw'].values)\n",
    "\n",
    "tmp = pd.DataFrame(D.groupby('hour_id')['gh'].mean())\n",
    "tmp.columns = ['gh_hour']\n",
    "D = D.merge(tmp, left_on='hour_id', right_index=True)\n",
    "\n",
    "tmp = pd.DataFrame(D.groupby('day_id')['gh'].mean())\n",
    "tmp.columns = ['gh_day']\n",
    "D = D.merge(tmp, left_on='day_id', right_index=True)\n",
    "\n",
    "tmp = pd.DataFrame(D.groupby('week_id')['gh'].mean())\n",
    "tmp.columns = ['gh_week']\n",
    "D = D.merge(tmp, left_on='week_id', right_index=True)\n",
    "\n",
    "D.to_csv('data.csv')"
   ]
  },
  {
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}