Mamore_migration_rates.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Import packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/Users/zoltan/Dropbox/Python\n"
     ]
    }
   ],
   "source": [
    "cd /Users/zoltan/Dropbox/Python"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "u'/Users/zoltan/Dropbox/Python'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pwd"
   ]
  },
  {
   "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": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<module 'cline_analysis' from 'cline_analysis.py'>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reload(ca)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "dirname = '/Users/zoltan/Dropbox/Channels/Fluvial/Bolivia/Mamore/csv_files_new/'\n",
    "fnames,clxs,clys,rbxs,lbxs,rbys,lbys,curvatures,ages,widths,dates = ca.load_data(dirname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Mamore_19840811.csv',\n",
       " 'Mamore_19861105.csv',\n",
       " 'Mamore_19890708.csv',\n",
       " 'Mamore_19900524.csv']"
      ]
     },
     "execution_count": 17,
     "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(1984, 8, 11, 0, 0),\n",
       " datetime.datetime(1986, 11, 5, 0, 0),\n",
       " datetime.datetime(1989, 7, 8, 0, 0),\n",
       " datetime.datetime(1990, 5, 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 = 1 # first timestep\n",
    "ts2 = 2 # 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",
    "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": 10,
   "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": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12063ce10>]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "curv,s = ca.compute_curvature(x,y)\n",
    "plt.figure()\n",
    "plt.plot(s,curv)\n",
    "curv = medfilt(savgol_filter(curv,71,3),kernel_size=5) # smoothing\n",
    "plt.plot(s,curv);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x126fa5510>]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(s,migr_rate)\n",
    "migr_rate = medfilt(savgol_filter(migr_rate,41,3),kernel_size=5) # smoothing\n",
    "plt.plot(s,migr_rate);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# set intervals affected by cu=toffs to NaN - specific to Mamore river \n",
    "migr_rate[907:1519] = np.NaN\n",
    "migr_rate[12053:12424] = np.NaN\n",
    "migr_rate[16809:17128] = np.NaN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1289baf90>]"
      ]
     },
     "execution_count": 45,
     "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('MamoreLT05_L1TP_232071_19861105_20170216_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 = [69,77,149,153,157,159]\n",
    "cutoff_inds = [61,100,101,107,116,147,148]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plot curvature and migration rate series side-by-side"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot curvature and migration rate along the channel\n",
    "\n",
    "W = np.nanmean(widths[1])\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(25,4))\n",
    "plt.tight_layout()\n",
    "\n",
    "curv_scale = 0.6\n",
    "migr_scale = 60\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate lag between curvature and migration rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1216cc490>]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot widths and boundary between the two segments\n",
    "plt.figure()\n",
    "plt.plot(s,widths[1])\n",
    "plt.plot([s[10516],s[10516]],[0,800],'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "226.15607678950275"
      ]
     },
     "execution_count": 22,
     "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[1][:10516])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "314.2718470671267"
      ]
     },
     "execution_count": 23,
     "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[1][10516:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "647.2772277227723"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(25.0*(LZM[:101]-LZC[:101]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "911.69354838709683"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(25.0*(LZM[101:]-LZC[101:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# First segment (Mamore A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate friction factor Cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x122b4d810>]"
      ]
     },
     "execution_count": 27,
     "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 = 1519\n",
    "i2 = 10516\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": 28,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00760703125\n"
     ]
    }
   ],
   "source": [
    "# this might take a while to run\n",
    "kl = 60.0 # preliminary kl value (guesstimate)\n",
    "k = 1\n",
    "W = np.mean(widths[1][i1: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",
    "\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": 29,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "Cf_opt = 0.00760703125"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate migration rate constant kl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30\n",
      "31\n",
      "32\n",
      "33\n",
      "34\n",
      "35\n",
      "36\n",
      "37\n",
      "38\n",
      "39\n",
      "40\n",
      "41\n",
      "42\n",
      "43\n",
      "44\n",
      "45\n",
      "46\n",
      "47\n",
      "48\n",
      "49\n",
      "50\n",
      "51\n",
      "52\n",
      "53\n",
      "54\n",
      "55\n",
      "56\n",
      "57\n",
      "58\n",
      "59\n",
      "60\n",
      "61\n",
      "62\n",
      "63\n",
      "64\n",
      "65\n",
      "66\n",
      "67\n",
      "68\n",
      "69\n",
      "70\n",
      "71\n",
      "72\n",
      "73\n",
      "74\n",
      "75\n",
      "76\n",
      "77\n",
      "78\n",
      "79\n"
     ]
    }
   ],
   "source": [
    "# minimize the error between actual and predicted migration rates (using the 75th percentile)\n",
    "errors = []\n",
    "for i in np.arange(30,80):\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(30,80),errors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "kl_opt = 60.0 # the error is at minimum for kl = 60.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "25.88"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "647/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": 36,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x1277e14d0>"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# kernel density and scatterplot of actual vs. nominal migration rate\n",
    "w = np.mean(widths[1][i1:i2]) \n",
    "\n",
    "curv_nodim = W*curv_t*kl_opt\n",
    "lag = 26\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 = 50.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": 37,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.778991262865\n",
      "0.606827387619\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": 38,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8971"
      ]
     },
     "execution_count": 38,
     "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": 39,
   "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": 40,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1273efed0>]"
      ]
     },
     "execution_count": 40,
     "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": 41,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.771376417234\n",
      "0.595021577065\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": 42,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "38.974032537357473"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 90th percentile of migration rate\n",
    "np.percentile(np.abs(m_nonan),90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x129b4a750>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 50.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": 44,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 70.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",
    "    if (i>=8) and (i<101):\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",
    "    if (i>=8) and (i<101):\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 (Mamore B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate friction factor Cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12ffbee50>]"
      ]
     },
     "execution_count": 50,
     "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 = 10516\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": 51,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x132d89ed0>]"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "migr_rate_t[12053-i1:12424-i1] = 0\n",
    "migr_rate_t[16809-i1:17128-i1] = 0\n",
    "plt.figure()\n",
    "plt.plot(migr_rate_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00448828125\n"
     ]
    }
   ],
   "source": [
    "# this might take a while to run\n",
    "kl = 100.0 # preliminary kl value (guesstimate)\n",
    "k = 1\n",
    "W = np.mean(widths[1][i1: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[12053-i1:12424-i1] = 0\n",
    "curv_t[16809-i1:17128-i1] = 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": 53,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "Cf_opt = 0.00448828125"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "227.33329569767974"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate migration rate constant kl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "120\n",
      "121\n",
      "122\n",
      "123\n",
      "124\n",
      "125\n",
      "126\n",
      "127\n",
      "128\n",
      "129\n",
      "130\n",
      "131\n",
      "132\n",
      "133\n",
      "134\n",
      "135\n",
      "136\n",
      "137\n",
      "138\n",
      "139\n",
      "140\n",
      "141\n",
      "142\n",
      "143\n",
      "144\n",
      "145\n",
      "146\n",
      "147\n",
      "148\n",
      "149\n",
      "150\n",
      "151\n",
      "152\n",
      "153\n",
      "154\n",
      "155\n",
      "156\n",
      "157\n",
      "158\n",
      "159\n"
     ]
    }
   ],
   "source": [
    "# minimize the error between actual and predicted migration rates (using the 75th percentile)\n",
    "errors = []\n",
    "for i in np.arange(120,160):\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(120,160),errors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "kl_opt = 128.0 # the error is at minimum for kl = 128.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "36.48"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "912/25.0 # lag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x134f25e50>]"
      ]
     },
     "execution_count": 60,
     "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": 68,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x12f8b5e10>"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "migr_rate_t[12053-i1:12424-i1] = np.NaN\n",
    "migr_rate_t[16809-i1:17128-i1] = np.NaN\n",
    "\n",
    "# kernel density and scatterplot of actual vs. nominal migration rate\n",
    "w = np.nanmean(widths[1][i1:i2]) \n",
    "\n",
    "curv_nodim = W*curv_t*kl_opt\n",
    "lag = 36\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 = 100.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": 69,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.641711040854\n",
      "0.411793059954\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": 70,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7341"
      ]
     },
     "execution_count": 70,
     "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": 71,
   "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": 72,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1229afb10>]"
      ]
     },
     "execution_count": 72,
     "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": 73,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.625278959609\n",
      "0.39097377733\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": 74,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "82.04226937039796"
      ]
     },
     "execution_count": 74,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 90th percentile of migration rate\n",
    "np.percentile(np.abs(m_nonan),90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x136733910>"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 100.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": 80,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 120.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",
    "    if i>=101:\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",
    "    if i>=101:\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
}