Tarauaca_migration_rates.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.signal import savgol_filter\n",
    "import cline_analysis as ca\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import datetime\n",
    "import os\n",
    "from scipy.signal import medfilt\n",
    "import functools\n",
    "from scipy.optimize import bisect\n",
    "from scipy import stats\n",
    "\n",
    "sns.set_style(\"whitegrid\")\n",
    "sns.set_style(\"ticks\")\n",
    "%matplotlib qt\n",
    "%config InlineBackend.figure_format = 'svg'\n",
    "plt.matplotlib.rcParams['svg.fonttype'] = 'svgfont' # fonts will be recognized by Adobe Illustrator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "dirname = '/Users/zoltan/Dropbox/Channels/Fluvial/Tarauaca/csv_files/'\n",
    "fnames,clxs,clys,rbxs,lbxs,rbys,lbys,curvatures,ages,widths,dates = ca.load_data(dirname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Tarauaca_20070729.csv',\n",
       " 'Tarauaca_20170724.csv',\n",
       " 'Tarauaca_Jurua_20070729.csv',\n",
       " 'Tarauaca_Jurua_20170724.csv']"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[datetime.datetime(2007, 7, 29, 0, 0),\n",
       " datetime.datetime(2017, 7, 24, 0, 0),\n",
       " datetime.datetime(2007, 7, 29, 0, 0),\n",
       " datetime.datetime(2017, 7, 24, 0, 0)]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Get migration rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ts1 = 0 # first timestep\n",
    "ts2 = 1 # second timestep\n",
    "\n",
    "d = dates[ts2]-dates[ts1]\n",
    "years = d.days/365.0\n",
    "\n",
    "x = np.array(clxs[ts1])\n",
    "y = np.array(clys[ts1])\n",
    "\n",
    "x = x[:12741] # get rid of extra length in x and y\n",
    "y = y[:12741]\n",
    "\n",
    "xn = np.array(clxs[ts2])\n",
    "yn = np.array(clys[ts2])\n",
    "\n",
    "migr_rate, migr_sign, p, q = ca.get_migr_rate(x,y,xn,yn,years,0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# qc\n",
    "plt.figure()\n",
    "plt.plot(x,y)\n",
    "plt.plot(xn,yn)\n",
    "plt.axis('equal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "migr_rate = medfilt(savgol_filter(migr_rate,41,3),kernel_size=5) # smoothing\n",
    "curv,s = ca.compute_curvature(x,y)\n",
    "curv = medfilt(savgol_filter(curv,71,3),kernel_size=5) # smoothing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# set intervals affected by cu=toffs to NaN - specific to Tarauaca river \n",
    "migr_rate[4672:4905] = np.NaN\n",
    "migr_rate[5461:5699] = np.NaN\n",
    "migr_rate[6812:6995] = np.NaN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x127e18290>]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# qc\n",
    "plt.figure()\n",
    "plt.plot(migr_rate)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Read 'valid' inflection points and corresponding points of zero migration from CSV file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('Tarauaca_LT05_L1TP_004065_20070729_20161112_01_T1_inflection_and_zero_migration_indices.csv')\n",
    "LZC = np.array(df['index of inflection point'])\n",
    "LZM = np.array(df['index of zero migration'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# indices of bends affected by low erodibility and cutoffs (these have been picked manually)\n",
    "erodibility_inds = [93,97,113,143,147]\n",
    "cutoff_inds = [47,48,66,67,142,159,160]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plot curvature and migration rate series side-by-side"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot curvature and migration rate along the channel\n",
    "\n",
    "W = np.nanmean(widths[0])\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(25,4))\n",
    "plt.tight_layout()\n",
    "\n",
    "curv_scale = 0.4\n",
    "migr_scale = 15\n",
    "y1 = curv_scale\n",
    "y2 = -3*curv_scale\n",
    "y3 = 3*migr_scale\n",
    "y4 = -migr_scale\n",
    "\n",
    "y5 = -2*curv_scale\n",
    "y6 = 2*migr_scale\n",
    "\n",
    "for i in range(0,len(LZC)-1,2):\n",
    "    xcoords = [s[LZC[i]],s[LZC[i+1]],s[LZC[i+1]],s[LZM[i+1]],s[LZM[i+1]],s[LZM[i]],s[LZM[i]],s[LZC[i]]]\n",
    "    ycoords = [y1,y1,0,y5,y2,y2,y5,0]\n",
    "    ax1.fill(xcoords,ycoords,facecolor=[0.85,0.85,0.85],edgecolor='k',zorder=0)\n",
    "\n",
    "deltas = 25.0\n",
    "ax1.fill_between(s, 0, curv*W)\n",
    "ax2 = ax1.twinx()\n",
    "ax2.fill_between(s, 0, migr_rate, facecolor='green')\n",
    "\n",
    "ax1.plot([0,max(s)],[0,0],'k--')\n",
    "ax2.plot([0,max(s)],[0,0],'k--')\n",
    "\n",
    "ax1.set_ylim(y2,y1)\n",
    "ax2.set_ylim(y4,y3)\n",
    "ax1.set_xlim(0,s[-1])\n",
    "\n",
    "for i in erodibility_inds:\n",
    "    xcoords = [s[LZC[i]],s[LZC[i+1]],s[LZC[i+1]],s[LZM[i+1]],s[LZM[i+1]],s[LZM[i]],s[LZM[i]],s[LZC[i]]]\n",
    "    ycoords = [y1,y1,0,y5,y2,y2,y5,0]\n",
    "    ax1.fill(xcoords,ycoords,facecolor=[1.0,0.85,0.85],edgecolor='k',zorder=0) \n",
    "    \n",
    "for i in cutoff_inds:\n",
    "    xcoords = [s[LZC[i]],s[LZC[i+1]],s[LZC[i+1]],s[LZM[i+1]],s[LZM[i+1]],s[LZM[i]],s[LZM[i]],s[LZC[i]]]\n",
    "    ycoords = [y1,y1,0,y5,y2,y2,y5,0]\n",
    "    ax1.fill(xcoords,ycoords,facecolor=[0.85,1.0,0.85],edgecolor='k',zorder=0) \n",
    "    \n",
    "for i in range(len(LZC)-1):\n",
    "    if np.sum(np.isnan(migr_rate[LZM[i]:LZM[i+1]]))>0:\n",
    "        xcoords = [s[LZC[i]],s[LZC[i+1]],s[LZC[i+1]],s[LZM[i+1]],s[LZM[i+1]],s[LZM[i]],s[LZM[i]],s[LZC[i]]]\n",
    "        ycoords = [y1,y1,0,y5,y2,y2,y5,0]\n",
    "        ax1.fill(xcoords,ycoords,color='w') \n",
    "        \n",
    "for i in range(len(LZC)-1):\n",
    "    if np.sum(np.isnan(migr_rate[LZM[i]:LZM[i+1]]))>0:\n",
    "        xcoords = [s[LZC[i]],s[LZC[i+1]],s[LZC[i+1]],s[LZM[i+1]],s[LZM[i+1]],s[LZM[i]],s[LZM[i]],s[LZC[i]]]\n",
    "        ycoords = [y3,y3,y6,0,y4,y4,0,y6]\n",
    "        ax2.fill(xcoords,ycoords,color='w') \n",
    "\n",
    "for i in range(0,len(LZC)-1,2):\n",
    "    ax1.text(s[LZC[i]],0.5,str(i),fontsize=12)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate lag between curvature and migration rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12631a490>]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot widths and boundary between the two segments\n",
    "plt.figure()\n",
    "plt.plot(s,widths[0][:12741])\n",
    "plt.plot([s[9475],s[9475]],[0,250],'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "94.58626630546298"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# first segment\n",
    "# average lag estimated from distances between inflection points and points of zero migration \n",
    "# (this is what was used in the paper)\n",
    "np.mean(widths[0][:9475])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "146.80917772320933"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# second segment\n",
    "# average lag estimated from distances between inflection points and points of zero migration \n",
    "# (this is what was used in the paper)\n",
    "np.mean(widths[0][9475:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "443.57142857142856"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(25.0*(LZM[:140]-LZC[:140]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "534.09090909090912"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(25.0*(LZM[140:]-LZC[140:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# First segment (Tarauaca A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate friction factor Cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12984d450>]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# first we need a continuous channel segment (e.g., no NaNs due to cutoffs)\n",
    "q=np.array(q)\n",
    "p=np.array(p)\n",
    "         \n",
    "i1 = 0\n",
    "i2 = 9475\n",
    "i1n = p[np.where(q==i1)[0][0]]\n",
    "i2n = p[np.where(q==i2)[0][0]]\n",
    "\n",
    "xt = x[i1:i2]\n",
    "yt = y[i1:i2]\n",
    "xnt = xn[i1n:i2n]\n",
    "ynt = yn[i1n:i2n]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(xt,yt)\n",
    "plt.plot(xnt,ynt)\n",
    "plt.axis('equal')\n",
    "\n",
    "migr_rate_t, migr_sign_t, pt, qt = ca.get_migr_rate(xt,yt,xnt,ynt,years,0)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(migr_rate_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x127edbc10>]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(migr_rate_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "migr_rate_t[4672:4905] = 0\n",
    "migr_rate_t[5461:5699] = 0\n",
    "migr_rate_t[6812:6995] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.003903515625\n"
     ]
    }
   ],
   "source": [
    "# this might take a while to run\n",
    "kl = 25.0 # preliminary kl value (guesstimate)\n",
    "k = 1\n",
    "W = np.mean(widths[0][:i2])\n",
    "D = (W/18.8)**0.7092 # depth in meters (from width)\n",
    "\n",
    "dx,dy,ds,s = ca.compute_derivatives(xt,yt)\n",
    "curv_t, s = ca.compute_curvature(xt,yt)\n",
    "curv_t = medfilt(savgol_filter(curv_t,71,3),kernel_size=5) # smoothing\n",
    "curv_t[4672:4905] = 0\n",
    "curv_t[5461:5699] = 0\n",
    "curv_t[6812:6995] = 0\n",
    "\n",
    "migr_rate_t = medfilt(savgol_filter(migr_rate_t,41,3),kernel_size=5)\n",
    "\n",
    "get_friction_factor_1 = functools.partial(ca.get_friction_factor,curvature=curv_t,migr_rate=migr_rate_t,\n",
    "                                          kl=kl,W=W, k=k, D=D, s=s)\n",
    "\n",
    "Cf_opt = bisect(get_friction_factor_1, 0.0002, 0.1)\n",
    "print Cf_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "Cf_opt = 0.003903515625"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate migration rate constant kl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n"
     ]
    }
   ],
   "source": [
    "# minimize the error between actual and predicted migration rates (using the 75th percentile)\n",
    "errors = []\n",
    "# curv_t, s = ca.compute_curvature(xt,yt)\n",
    "# curv_t = medfilt(savgol_filter(curv_t,71,3),kernel_size=5) # smoothing\n",
    "for i in np.arange(10,40):\n",
    "    print i\n",
    "    R1 = ca.get_predicted_migr_rate(curv_t,W=W,k=1,Cf=Cf_opt,D=D,kl=i,s=s)\n",
    "    errors.append(np.abs(np.percentile(np.abs(R1),75)-np.percentile(np.abs(migr_rate_t[1:-1]),75)))\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(np.arange(10,40),errors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "kl_opt = 31.0 # the error is at minimum for kl = 31.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "17.76"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "444/25.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x125cb8810>]"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(W*kl_opt*curv_t)\n",
    "plt.plot(migr_rate_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plot actual migration rate against nominal migration rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1301705d0>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "migr_rate_t[4672:4905] = np.NaN\n",
    "migr_rate_t[5461:5699] = np.NaN\n",
    "migr_rate_t[6812:6995] = np.NaN\n",
    "\n",
    "# kernel density and scatterplot of actual vs. nominal migration rate\n",
    "w = np.nanmean(widths[0][:9475]) \n",
    "\n",
    "curv_nodim = W*curv_t*kl_opt\n",
    "lag = 18\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.kdeplot(curv_nodim[:-lag][np.isnan(migr_rate_t[lag:])==0], migr_rate_t[lag:][np.isnan(migr_rate_t[lag:])==0],\n",
    "           n_levels=20,shade=True,cmap='Blues',shade_lowest=False)\n",
    "plt.scatter(curv_nodim[:-lag][::20],migr_rate_t[lag:][::20],c='k',s=15)\n",
    "max_x = 17.0\n",
    "plt.xlim(-max_x,max_x)\n",
    "plt.ylim(-max_x,max_x)\n",
    "plt.plot([-max_x,max_x],[-max_x,max_x],'k--') \n",
    "plt.xlabel('nominal migration rate (m/year)', fontsize=14)\n",
    "plt.ylabel('actual migration rate (m/year)', fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.645999757937\n",
      "0.417315687255\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "# get correlation coefficient for relationship between curvature and migration rate\n",
    "slope, intercept, r_value, p_value, slope_std_rror = stats.linregress(curv_nodim[:-lag][np.isnan(migr_rate_t[lag:])==0],\n",
    "                                                                      migr_rate_t[lag:][np.isnan(migr_rate_t[lag:])==0])\n",
    "print r_value\n",
    "print r_value**2\n",
    "print p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8803"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# number of data points used in analysis\n",
    "len(curv_nodim[:-lag][np.isnan(migr_rate_t[lag:])==0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# compute predicted migration rates\n",
    "D = (w/18.8)**0.7092 # depth in meters (from width)\n",
    "dx,dy,ds,s = ca.compute_derivatives(xt,yt)\n",
    "R1 = ca.get_predicted_migr_rate(curv_t,W=w,k=1,Cf=Cf_opt,D=D,kl=kl_opt,s=s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12e7ad0d0>]"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual and predicted migration rates\n",
    "plt.figure()\n",
    "plt.plot(s,migr_rate_t)\n",
    "plt.plot(s,R1,'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.636577276786\n",
      "0.40523062932\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "# get correlation coefficient for relationship between actual and predicted migration rate\n",
    "m_nonan = migr_rate_t[(np.isnan(R1)==0)&(np.isnan(migr_rate_t)==0)]\n",
    "R_nonan = R1[(np.isnan(R1)==0)&(np.isnan(migr_rate_t)==0)]\n",
    "\n",
    "slope, intercept, r_value, p_value, slope_std_rror = stats.linregress(R_nonan,m_nonan)\n",
    "print r_value\n",
    "print r_value**2\n",
    "print p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.710110374800637"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 90th percentile of migration rate\n",
    "np.percentile(np.abs(m_nonan),90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x128e54150>"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 17.0\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.kdeplot(R_nonan,m_nonan,n_levels=10,shade=True,cmap='Blues',shade_lowest=False)\n",
    "plt.plot([-max_m,max_m],[-max_m,max_m],'k--') \n",
    "plt.scatter(R_nonan[::20],m_nonan[::20],c='k',s=15)\n",
    "plt.xlim(-max_m,max_m)\n",
    "plt.ylim(-max_m,max_m)\n",
    "plt.xlabel('predicted migration rate (m/year)', fontsize=14)\n",
    "plt.ylabel('actual migration rate (m/year)', fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 20.0\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.kdeplot(R_nonan,m_nonan,n_levels=10,shade=True,cmap='Blues',shade_lowest=False)\n",
    "plt.plot([-max_m,max_m],[-max_m,max_m],'k--') \n",
    "plt.scatter(R_nonan[::20],m_nonan[::20],c='k',s=15)\n",
    "plt.xlim(-max_m,max_m)\n",
    "plt.ylim(-max_m,max_m)\n",
    "plt.xlabel('predicted migration rate (m/year)', fontsize=14)\n",
    "plt.ylabel('actual migration rate (m/year)', fontsize=14)\n",
    "\n",
    "# add points affected by cutoffs and low erodibility\n",
    "for i in erodibility_inds:\n",
    "    plt.scatter(R1[-i1+LZC[i]:-i1+LZC[i+1]][::10],migr_rate_t[-i1+LZC[i]:-i1+LZC[i+1]][::10],c='r',s=15)\n",
    "for i in cutoff_inds:\n",
    "    plt.scatter(R1[-i1+LZC[i]:-i1+LZC[i+1]][::10],migr_rate_t[-i1+LZC[i]:-i1+LZC[i+1]][::10],c='g',s=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Second segment (Tarauaca B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate friction factor Cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13454e410>]"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# first we need a continuous channel segment (e.g., no NaNs due to cutoffs)\n",
    "q=np.array(q)\n",
    "p=np.array(p)\n",
    "         \n",
    "i1 = 9475\n",
    "i2 = len(x)-1\n",
    "i1n = p[np.where(q==i1)[0][0]]\n",
    "i2n = p[np.where(q==i2)[0][0]]\n",
    "\n",
    "xt = x[i1:i2]\n",
    "yt = y[i1:i2]\n",
    "xnt = xn[i1n:i2n]\n",
    "ynt = yn[i1n:i2n]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(xt,yt)\n",
    "plt.plot(xnt,ynt)\n",
    "plt.axis('equal')\n",
    "\n",
    "migr_rate_t, migr_sign_t, pt, qt = ca.get_migr_rate(xt,yt,xnt,ynt,years,0)\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(migr_rate_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.004683203125\n"
     ]
    }
   ],
   "source": [
    "# this might take a while to run\n",
    "kl = 25.0 # preliminary kl value (guesstimate)\n",
    "k = 1\n",
    "W = np.mean(widths[0][i1:])\n",
    "D = (W/18.8)**0.7092 # depth in meters (from width)\n",
    "\n",
    "dx,dy,ds,s = ca.compute_derivatives(xt,yt)\n",
    "curv_t, s = ca.compute_curvature(xt,yt)\n",
    "curv_t = medfilt(savgol_filter(curv_t,71,3),kernel_size=5) # smoothing\n",
    "\n",
    "migr_rate_t = medfilt(savgol_filter(migr_rate_t,41,3),kernel_size=5)\n",
    "\n",
    "get_friction_factor_1 = functools.partial(ca.get_friction_factor,curvature=curv_t,migr_rate=migr_rate_t,\n",
    "                                          kl=kl,W=W, k=k, D=D, s=s)\n",
    "\n",
    "Cf_opt = bisect(get_friction_factor_1, 0.0002, 0.1)\n",
    "print Cf_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "Cf_opt = 0.004683203125"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate migration rate constant kl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10\n",
      "11\n",
      "12\n",
      "13\n",
      "14\n",
      "15\n",
      "16\n",
      "17\n",
      "18\n",
      "19\n",
      "20\n",
      "21\n",
      "22\n",
      "23\n",
      "24\n",
      "25\n",
      "26\n",
      "27\n",
      "28\n",
      "29\n",
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n"
     ]
    }
   ],
   "source": [
    "# minimize the error between actual and predicted migration rates (using the 75th percentile)\n",
    "errors = []\n",
    "for i in np.arange(10,40):\n",
    "    print i\n",
    "    R1 = ca.get_predicted_migr_rate(curv_t,W=W,k=1,Cf=Cf_opt,D=D,kl=i,s=s)\n",
    "    errors.append(np.abs(np.percentile(np.abs(R1),75)-np.percentile(np.abs(migr_rate_t[1:-1]),75)))\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(np.arange(10,40),errors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "kl_opt = 25.0 # the error is at minimum for kl = 25.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "21.36"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "534/25.0 # lag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x146fa6cd0>]"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(W*kl_opt*curv_t)\n",
    "plt.plot(migr_rate_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plot actual migration rate against nominal migration rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x126138310>"
      ]
     },
     "execution_count": 63,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# kernel density and scatterplot of actual vs. nominal migration rate\n",
    "w = np.nanmean(widths[0][i1:]) \n",
    "\n",
    "curv_nodim = W*curv_t*kl_opt\n",
    "lag = 21\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.kdeplot(curv_nodim[:-lag][np.isnan(migr_rate_t[lag:])==0], migr_rate_t[lag:][np.isnan(migr_rate_t[lag:])==0],\n",
    "           n_levels=20,shade=True,cmap='Blues',shade_lowest=False)\n",
    "plt.scatter(curv_nodim[:-lag][::20],migr_rate_t[lag:][::20],c='k',s=15)\n",
    "max_x = 15.0\n",
    "plt.xlim(-max_x,max_x)\n",
    "plt.ylim(-max_x,max_x)\n",
    "plt.plot([-max_x,max_x],[-max_x,max_x],'k--') \n",
    "plt.xlabel('nominal migration rate (m/year)', fontsize=14)\n",
    "plt.ylabel('actual migration rate (m/year)', fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.806775498977\n",
      "0.650886705749\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "# get correlation coefficient for relationship between curvature and migration rate\n",
    "slope, intercept, r_value, p_value, slope_std_rror = stats.linregress(curv_nodim[:-lag][np.isnan(migr_rate_t[lag:])==0],\n",
    "                                                                      migr_rate_t[lag:][np.isnan(migr_rate_t[lag:])==0])\n",
    "print r_value\n",
    "print r_value**2\n",
    "print p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3244"
      ]
     },
     "execution_count": 65,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# number of data points used in analysis\n",
    "len(curv_nodim[:-lag][np.isnan(migr_rate_t[lag:])==0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# compute predicted migration rates\n",
    "D = (w/18.8)**0.7092 # depth in meters (from width)\n",
    "dx,dy,ds,s = ca.compute_derivatives(xt,yt)\n",
    "R1 = ca.get_predicted_migr_rate(curv_t,W=w,k=1,Cf=Cf_opt,D=D,kl=kl_opt,s=s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x128335bd0>]"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual and predicted migration rates\n",
    "plt.figure()\n",
    "plt.plot(s,migr_rate_t)\n",
    "plt.plot(s,R1,'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.789850118104\n",
      "0.623863209069\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "# get correlation coefficient for relationship between actual and predicted migration rate\n",
    "m_nonan = migr_rate_t[(np.isnan(R1)==0)&(np.isnan(migr_rate_t)==0)]\n",
    "R_nonan = R1[(np.isnan(R1)==0)&(np.isnan(migr_rate_t)==0)]\n",
    "\n",
    "slope, intercept, r_value, p_value, slope_std_rror = stats.linregress(R_nonan,m_nonan)\n",
    "print r_value\n",
    "print r_value**2\n",
    "print p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10.223005238794242"
      ]
     },
     "execution_count": 69,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 90th percentile of migration rate\n",
    "np.percentile(np.abs(m_nonan),90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x129cf3850>"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 15.0\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.kdeplot(R_nonan,m_nonan,n_levels=10,shade=True,cmap='Blues',shade_lowest=False)\n",
    "plt.plot([-max_m,max_m],[-max_m,max_m],'k--') \n",
    "plt.scatter(R_nonan[::20],m_nonan[::20],c='k',s=15)\n",
    "plt.xlim(-max_m,max_m)\n",
    "plt.ylim(-max_m,max_m)\n",
    "plt.xlabel('predicted migration rate (m/year)', fontsize=14)\n",
    "plt.ylabel('actual migration rate (m/year)', fontsize=14)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 15.0\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.kdeplot(R_nonan,m_nonan,n_levels=10,shade=True,cmap='Blues',shade_lowest=False)\n",
    "plt.plot([-max_m,max_m],[-max_m,max_m],'k--') \n",
    "plt.scatter(R_nonan[::20],m_nonan[::20],c='k',s=15)\n",
    "plt.xlim(-max_m,max_m)\n",
    "plt.ylim(-max_m,max_m)\n",
    "plt.xlabel('predicted migration rate (m/year)', fontsize=14)\n",
    "plt.ylabel('actual migration rate (m/year)', fontsize=14)\n",
    "\n",
    "# add points affected by cutoffs and low erodibility\n",
    "for i in erodibility_inds:\n",
    "    plt.scatter(R1[-i1+LZC[i]:-i1+LZC[i+1]][::5],migr_rate_t[-i1+LZC[i]:-i1+LZC[i+1]][::5],c='r',s=15)\n",
    "for i in cutoff_inds:\n",
    "    plt.scatter(R1[-i1+LZC[i]:-i1+LZC[i+1]][::5],migr_rate_t[-i1+LZC[i]:-i1+LZC[i+1]][::5],c='g',s=15)"
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}