Yavari_tributary_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": 5,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "dirname = '/Users/zoltan/Dropbox/Channels/Fluvial/Yavari/csv_files/'\n",
    "fnames,clxs,clys,rbxs,lbxs,rbys,lbys,curvatures,ages,widths,dates = ca.load_data(dirname)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['Yavari_19880715.csv',\n",
       " 'Yavari_20160813.csv',\n",
       " 'Yavari_tributary_1_19880715.csv',\n",
       " 'Yavari_tributary_1_20160813.csv']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fnames"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[datetime.datetime(1988, 7, 15, 0, 0),\n",
       " datetime.datetime(2016, 8, 13, 0, 0),\n",
       " datetime.datetime(1988, 7, 15, 0, 0),\n",
       " datetime.datetime(2016, 8, 13, 0, 0)]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Get migration rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "ts1 = 2 # first timestep\n",
    "ts2 = 3 # 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": 8,
   "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": 9,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# set intervals affected by cu=toffs to NaN - specific to Yavari tributary river segment\n",
    "migr_rate[1694:1744] = np.NaN\n",
    "migr_rate[2444:2665] = np.NaN\n",
    "migr_rate[3442:3595] = np.NaN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1272c07d0>]"
      ]
     },
     "execution_count": 10,
     "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": 11,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "df = pd.read_csv('Yavari_LT05_L1TP_005063_19880715_20170208_01_T1_tributary_1_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": 12,
   "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 = [66,81,83,85,98,105,118,120,124,130]\n",
    "cutoff_inds = [26,27,36,45,48,49,93,94,122,123]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plot curvature and migration rate series side-by-side"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "# plot curvature and migration rate along the channel\n",
    "\n",
    "W = np.nanmean(widths[2])\n",
    "\n",
    "fig, ax1 = plt.subplots(figsize=(25,4))\n",
    "plt.tight_layout()\n",
    "\n",
    "curv_scale = 0.4\n",
    "migr_scale = 2\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.4,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": 17,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "window_length = 500\n",
    "time_shifts = ca.get_time_shifts(migr_rate,curv,window_length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x12678df90>]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(time_shifts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-450.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# average lag\n",
    "25.0*np.round(np.mean(time_shifts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "406.98529411764707"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 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(25.0*(LZM-LZC))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate friction factor Cf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "migr_rate_t = migr_rate.copy()\n",
    "migr_rate_t[1694:1744] = 0\n",
    "migr_rate_t[2444:2665] = 0\n",
    "migr_rate_t[3442:3595] = 0\n",
    "\n",
    "curv_t = curv.copy()\n",
    "curv_t[1694:1744] = 0\n",
    "curv_t[2444:2665] = 0\n",
    "curv_t[3442:3595] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.00526796875\n"
     ]
    }
   ],
   "source": [
    "# this might take a while to run\n",
    "kl = 3.0 # preliminary kl value (guesstimate)\n",
    "k = 1\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,71,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": 32,
   "metadata": {
    "collapsed": true,
    "deletable": true,
    "editable": true
   },
   "outputs": [],
   "source": [
    "Cf_opt = 0.00526796875"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Estimate migration rate constant kl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "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",
    "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": 34,
   "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": 28,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "16.28"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "407/25.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": true,
    "editable": true
   },
   "source": [
    "## Plot actual migration rate against nominal migration rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11c0d6c90>"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# kernel density and scatterplot of actual vs. nominal migration rate\n",
    "w = np.nanmedian(widths[2]) \n",
    "curv_nodim = w*curv*kl_opt\n",
    "lag = 16\n",
    "plt.figure(figsize=(8,8))\n",
    "sns.kdeplot(curv_nodim[:-lag][np.isnan(migr_rate[lag:])==0], migr_rate[lag:][np.isnan(migr_rate[lag:])==0],\n",
    "           n_levels=20,shade=True,cmap='Blues',shade_lowest=False)\n",
    "plt.scatter(curv_nodim[:-lag][::20],migr_rate[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": 44,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.784497877311\n",
      "0.615436919505\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[lag:])==0],\n",
    "                                                                      migr_rate[lag:][np.isnan(migr_rate[lag:])==0])\n",
    "print r_value\n",
    "print r_value**2\n",
    "print p_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8668"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# number of data points used in analysis\n",
    "len(curv_nodim[:-lag][np.isnan(migr_rate[lag:])==0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "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(x,y)\n",
    "R1 = ca.get_predicted_migr_rate(curv,W=w,k=1,Cf=Cf_opt,D=D,kl=kl_opt,s=s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x128da49d0>]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plot actual and predicted migration rates\n",
    "plt.figure()\n",
    "plt.plot(s,migr_rate)\n",
    "plt.plot(s,R1,'r')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.773642177706\n",
      "0.598522219126\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "# get correlation coefficient for relationship between actual and predicted migration rate\n",
    "m_nonan = migr_rate[(np.isnan(R1)==0)&(np.isnan(migr_rate)==0)]\n",
    "R_nonan = R1[(np.isnan(R1)==0)&(np.isnan(migr_rate)==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": 41,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.7905648228134488"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 90th percentile of migration rate\n",
    "np.percentile(np.abs(m_nonan),90)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": false,
    "deletable": true,
    "editable": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x11c0d6990>"
      ]
     },
     "execution_count": 45,
     "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": 48,
   "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[LZC[i]:LZC[i+1]][::5],migr_rate[LZC[i]:LZC[i+1]][::5],c='r',s=15)\n",
    "for i in cutoff_inds:\n",
    "    plt.scatter(R1[LZC[i]:LZC[i+1]][::5],migr_rate[LZC[i]: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
}