Jutai_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/Jutai/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": [
       "['Jutai_19890805.csv', 'Jutai_20170802.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(1989, 8, 5, 0, 0), datetime.datetime(2017, 8, 2, 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",
    "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": 110,
   "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": 111,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# set intervals affected by cu=toffs to NaN - specific to Jutai river \n",
    "migr_rate[1086:1293] = np.NaN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x142ddab10>]"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "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": 113,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('Jutai_LT05_L1TP_003063_19890805_20170202_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": 114,
   "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,115,117,119,163,189,191,204,218]\n",
    "cutoff_inds = [7,8,9,14,15,29,30,50,51,58,59,185,194,209,210]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plot curvature and migration rate series side-by-side"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot curvature and migration rate along the channel\n",
    "\n",
    "W = np.nanmean(widths[0]) # mean channel width\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(25,4))\n",
    "plt.tight_layout()\n",
    "\n",
    "curv_scale = 0.6\n",
    "migr_scale = 3\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": 116,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12aa06d10>]"
      ]
     },
     "execution_count": 116,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot widths and boundary between the two segments\n",
    "plt.figure()\n",
    "plt.plot(s,widths[0])\n",
    "plt.plot([s[9846],s[9846]],[0,800],'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "123.26877138353407"
      ]
     },
     "execution_count": 117,
     "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][:9846])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "260.77728739097813"
      ]
     },
     "execution_count": 118,
     "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][9846:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "309.73154362416108"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(25.0*(LZM[:149]-LZC[:149]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "552.08333333333337"
      ]
     },
     "execution_count": 120,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.mean(25.0*(LZM[149:]-LZC[149:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# First segment (Jutai A)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate friction factor Cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12b193c90>]"
      ]
     },
     "execution_count": 121,
     "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 = 1293\n",
    "i2 = 9846\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": 122,
   "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 = 3.0 # preliminary kl value (guesstimate)\n",
    "k = 1\n",
    "W = np.mean(widths[0][:9846])\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": 39,
   "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": 123,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\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(1,10):\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(1,10),errors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "kl_opt = 4.0 # the error is at minimum for kl = 4.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12.4"
      ]
     },
     "execution_count": 124,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "310/25.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13b9df190>]"
      ]
     },
     "execution_count": 125,
     "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": 126,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x12d919810>"
      ]
     },
     "execution_count": 126,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# kernel density and scatterplot of actual vs. nominal migration rate\n",
    "w = np.nanmean(widths[0][:9846]) \n",
    "\n",
    "curv_nodim = W*curv_t*kl_opt\n",
    "lag = 12\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 = 2.5\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": 127,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.860343681877\n",
      "0.740191250946\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": 128,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8541"
      ]
     },
     "execution_count": 128,
     "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": 129,
   "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": 130,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12dcb2290>]"
      ]
     },
     "execution_count": 130,
     "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": 131,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.853626448057\n",
      "0.728678112822\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": 132,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.746666174998416"
      ]
     },
     "execution_count": 132,
     "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 0x136ad7110>"
      ]
     },
     "execution_count": 70,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 2.5\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": 133,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 4.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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "# Second segment (Jutai B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate friction factor Cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12b3e0c50>]"
      ]
     },
     "execution_count": 82,
     "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 = 9846\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": 83,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00682734375\n"
     ]
    }
   ],
   "source": [
    "# this might take a while to run\n",
    "kl = 4.0 # preliminary kl value (guesstimate)\n",
    "k = 1\n",
    "W = np.mean(widths[0][9846:])\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": 84,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "Cf_opt = 0.00682734375"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate migration rate constant kl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n",
      "2\n",
      "3\n",
      "4\n",
      "5\n",
      "6\n",
      "7\n",
      "8\n",
      "9\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(1,10):\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(1,10),errors);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "kl_opt = 4.0 # the error is at minimum for kl = 4.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "22.08"
      ]
     },
     "execution_count": 87,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "552/25.0 # lag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13b92d950>]"
      ]
     },
     "execution_count": 88,
     "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": 90,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x13b8e7490>"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# kernel density and scatterplot of actual vs. nominal migration rate\n",
    "w = np.nanmean(widths[0][9846:]) \n",
    "\n",
    "curv_nodim = W*curv_t*kl_opt\n",
    "lag = 22\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 = 3.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": 91,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.717197986136\n",
      "0.514372951318\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": 92,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8045"
      ]
     },
     "execution_count": 92,
     "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": 93,
   "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": 94,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x13b8e7090>]"
      ]
     },
     "execution_count": 94,
     "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": 95,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.691878609838\n",
      "0.478696010751\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": 96,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.2412418121523747"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 90th percentile of migration rate\n",
    "np.percentile(np.abs(m_nonan),90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x13dcacc10>"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 3.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": 101,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot actual vs. predicted migration rate\n",
    "max_m = 5.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)"
   ]
  }
 ],
 "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
}