tutorial now commented

kbarnhart · kbarnhart · commit d858566f59e3 · 2018-05-01T09:49:23.000-06:00
diff --git a/rockblock/introduction_to_rockblock.ipynb b/rockblock/introduction_to_rockblock.ipynb
@@ -6,11 +6,11 @@
    "source": [
     "# Introduction to the RockBlock and LayeredRockBlock objects\n",
     "\n",
-    "The RockBlock and LayeredRockBlock are two landlab objects meant to make it easier to work with spatially variable lithology that produces spatially variable parameter values (e.g. stream power erodability and diffusivity). \n",
+    "The RockBlock and LayeredRockBlock are two landlab objects meant to make it easier to work with spatially variable lithology that produces spatially variable parameter values (e.g. stream power erodability or diffusivity). \n",
     "\n",
     "In this tutorial we will first use the LayeredRockBlock to erode either dipping layeres or an anticline. Then we will use RockBlock to create inverted topography. \n",
     "\n",
-    "To start, we will import the necessary modules. A note: this tutorial uses the [HoloViews package](http://holoviews.org) for visualization. This package is a great tool for dealing with multidimentional annotated data (e.g. an xarray dataset). If you get an error on import, consider updating dask (this is what the author needed to do in April 2018).  "
+    "To start, we will import the necessary modules. A note: this tutorial uses the [HoloViews package](http://holoviews.org) for visualization. This package is a great tool for dealing with multidimentional annotated data (e.g. an xarray dataset). If you get an error on import, consider updating dask (this is what the author needed to do in April 2018). You will also need to have the [Bokeh](https://bokeh.pydata.org/en/latest/) and [Matplotlib](https://matplotlib.org) packages installed."
    ]
   },
   {
@@ -30,7 +30,7 @@
     "\n",
     "from landlab import RasterModelGrid\n",
     "from landlab.layers import RockBlock, LayeredRockBlock\n",
-    "from landlab.components import FlowAccumulator, FastscapeEroder\n",
+    "from landlab.components import FlowAccumulator, FastscapeEroder, LinearDiffuser\n",
     "from landlab.plot import imshow_grid"
    ]
   },
@@ -64,7 +64,7 @@
     "attrs = {'K_sp': {0: 0.0003,\n",
     "                  1: 0.0001}}\n",
     "\n",
-    "mg = RasterModelGrid((50, 50), 100)\n",
+    "mg = RasterModelGrid((100, 100), 200)\n",
     "z = mg.add_zeros('node', 'topographic__elevation') \n",
     "random_field = 0.01*np.random.randn(mg.size('node'))\n",
     "z += random_field - random_field.min()\n",
@@ -74,7 +74,7 @@
     "ids = np.tile([0,1], 20) \n",
     "\n",
     "# Anticline\n",
-    "rb = LayeredRockBlock(mg, z0s, ids, x0=2500, y0=2500, function = lambda x, y : ((0.005*x)**2+(0.005*y)**2), attrs=attrs)\n",
+    "rb = LayeredRockBlock(mg, z0s, ids, x0=10000, y0=10000, function = lambda x, y : ((0.0015*x)**2+(0.0007*y)**2), attrs=attrs)\n",
     "\n",
     "# Shallow dips\n",
     "#rb = LayeredRockBlock(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.001*x)+(0.003*y)), attrs=attrs)\n",
@@ -117,25 +117,21 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "nts = 400\n",
+    "nts = 500\n",
     "U = 0.001\n",
     "dt = 1000\n",
     "\n",
     "fa = FlowAccumulator(mg)\n",
     "sp = FastscapeEroder(mg, K_sp='K_sp')\n",
     "\n",
     "out_fields = ['topographic__elevation',\n",
-    "              'K_sp', \n",
     "              'rock_type__id']\n",
     "\n",
-    "ds = xr.Dataset(data_vars={'topographic__elevation' : (('y', 'x', 't'), np.empty((mg.shape[0], mg.shape[1], nts))),\n",
-    "                           'K_sp' : (('y', 'x', 't'), np.empty((mg.shape[0], mg.shape[1], nts))),\n",
-    "                           'rock_type__id': (('y', 'x', 't'), np.empty((mg.shape[0], mg.shape[1], nts)))},\n",
+    "ds = xr.Dataset(data_vars={'topographic__elevation' : (('t', 'y', 'x'), np.empty((nts, mg.shape[0], mg.shape[1]))),\n",
+    "                           'rock_type__id': (('t', 'y', 'x'), np.empty((nts, mg.shape[0], mg.shape[1])))},\n",
     "                coords={'x': mg.x_of_node.reshape(mg.shape)[0,:],\n",
     "                        'y': mg.y_of_node.reshape(mg.shape)[:, 1],\n",
     "                        't': dt*np.arange(nts)/1e6})"
@@ -147,6 +143,8 @@
    "source": [
     "Here we run the model. In each time step we first run the FlowAccumulator to direct flow and accumulatate drainage area. Then the FastscapeEroder erodes the topography based on the stream power equation using the erodability value in the field`'K_sp'`. We create an uplift field that uplifts only the model grid's core nodes. After uplifting these core nodes, we update the LayeredRockBlock. Importantly, we must tell the LayeredRockBlock how it has been advected upward by uplift. \n",
     "\n",
+    "`rb.run_one_step` has an optional argument `rock_id` to use when some material may be deposited. Since we are using the FastscapeEroder which is fully detachment limited, we don't need to set this. \n",
+    "\n",
     "Within each timestep we save information about the model for plotting. "
    ]
   },
@@ -167,7 +165,7 @@
     "    rb.run_one_step(dz_advection = dz_ad)\n",
     "    \n",
     "    for of in out_fields:\n",
-    "        ds[of][:,:,i] = mg['node'][of].reshape(mg.shape)\n"
+    "        ds[of][i,:,:] = mg['node'][of].reshape(mg.shape)"
    ]
   },
   {
@@ -211,19 +209,16 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Next we specify that we want two images, one showing rock type and one showing topographic elevation. A slider bar "
+    "Next we specify that we want two images, one showing rock type and one showing topographic elevation. A slider bar shows us model time in millions of years.  "
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "%opts Image (cmap='viridis')\n",
-    "\n",
     "topo = hvds.to(hv.Image, ['x', 'y'], ['topographic__elevation'])\n",
     "rock = hvds.to(hv.Image, ['x', 'y'], ['rock_type__id'])\n",
     "\n",
@@ -236,64 +231,81 @@
     "collapsed": true
    },
    "source": [
-    "Make inverted topography by filling RockBlock in with resistant material only where the DA is large. "
+    "We can see the form of the anticline advecting through the topography. Cool!\n",
+    "\n",
+    "\n",
+    "## Part 2 Creation of Inverted Topography\n",
+    "\n",
+    "Here we will explore making inverted topography by eroding a RockBlock with constant properties for half of the model evaluation time, and then filling RockBlock in with resistant material only where the drainage area is large. This is meant as a simple example of filling in valleys with volcanic material. \n",
+    "\n",
+    "All of the details of the options for creating a [RockBlock](http://landlab.readthedocs.io/en/release/landlab.layers.rockblock.html) can be found here. \n",
+    "\n",
+    "In the next code block we make a new model and run it. Note that because we are using the LinearDiffuser which may result in depositing material we must specify a rock_id in the RockBlock `run_one_step` method.\n",
+    "\n",
+    "We also are handling the model grid boundary conditions differently than in the last example, setting the boundaries on the left and right sides to closed. "
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "mg2 = RasterModelGrid((50, 50), 100)\n",
+    "mg2 = RasterModelGrid((100, 100), 200)\n",
+    "mg2.set_closed_boundaries_at_grid_edges(True, False, True, False)\n",
     "z2 = mg2.add_zeros('node', 'topographic__elevation') \n",
     "random_field = 0.01*np.random.randn(mg2.size('node'))\n",
     "z2 += random_field - random_field.min()\n",
     "\n",
-    "thicknesses2 = [1000]\n",
+    "thicknesses2 = [10000]\n",
     "ids2 =[0]\n",
     "\n",
     "attrs2 = {'K_sp': {0: 0.0001,\n",
-    "                   1: 0.00001}}\n",
+    "                   1: 0.00001},\n",
+    "          'D': {0: 0.4,\n",
+    "                1: 0.001}}\n",
     "\n",
     "rb2 = RockBlock(mg2, thicknesses2, ids2, attrs=attrs2)\n",
     "\n",
-    "nts = 400\n",
-    "U = 0.001\n",
+    "nts = 500\n",
+    "U = 0.005\n",
     "dt = 1000\n",
     "\n",
     "fa2 = FlowAccumulator(mg2)\n",
     "sp2 = FastscapeEroder(mg2, K_sp='K_sp')\n",
+    "ld2 = LinearDiffuser(mg2, linear_diffusivity='D')\n",
     "\n",
     "out_fields = ['topographic__elevation',\n",
-    "              'K_sp', \n",
     "              'rock_type__id']\n",
     "\n",
-    "ds2 = xr.Dataset(data_vars={'topographic__elevation' : (('y', 'x', 't'), np.empty((mg2.shape[0], mg2.shape[1], nts))),\n",
-    "                            'K_sp' : (('y', 'x', 't'), np.empty((mg2.shape[0], mg2.shape[1], nts))),\n",
-    "                            'rock_type__id': (('y', 'x', 't'), np.empty((mg2.shape[0], mg2.shape[1], nts)))},\n",
-    "                 coords={'x': mg2.x_of_node.reshape(mg2.shape)[0,:],\n",
-    "                         'y': mg2.y_of_node.reshape(mg2.shape)[:, 1],\n",
-    "                         't': dt*np.arange(nts)/1e6})\n",
+    "out_fields = ['topographic__elevation',\n",
+    "              'rock_type__id']\n",
+    "\n",
+    "ds2 = xr.Dataset(data_vars={'topographic__elevation' : (('t', 'y', 'x'), np.empty((nts, mg2.shape[0], mg2.shape[1]))),\n",
+    "                           'rock_type__id': (('t', 'y', 'x'), np.empty((nts, mg2.shape[0], mg2.shape[1])))},\n",
+    "                coords={'x': mg2.x_of_node.reshape(mg2.shape)[0,:],\n",
+    "                        'y': mg2.y_of_node.reshape(mg2.shape)[:, 1],\n",
+    "                        't': dt*np.arange(nts)/1e6})\n",
     "half_nts = int(nts/2)\n",
     "for i in range(half_nts):\n",
     "    fa2.run_one_step()\n",
-    "    sp2.run_one_step(dt = dt)\n",
+    "    sp2.run_one_step(dt=dt)\n",
+    "    ld2.run_one_step(dt=dt)\n",
     "    dz_ad2 = np.zeros(mg2.size('node'))\n",
     "    dz_ad2[mg2.core_nodes] = U * dt\n",
     "    z2 += dz_ad2\n",
-    "    rb2.run_one_step(dz_advection = dz_ad2)\n",
+    "    rb2.run_one_step(dz_advection = dz_ad2, rock_id=0)\n",
     "    \n",
     "    for of in out_fields:\n",
-    "        ds2[of][:,:,i] = mg2['node'][of].reshape(mg2.shape)\n"
+    "        ds2[of][i,:,:] = mg2['node'][of].reshape(mg2.shape)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "After the first half of run time, let's look at the topography. "
+   ]
   },
   {
    "cell_type": "code",
@@ -307,52 +319,64 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "We can see that we have developed ridges and valleys as we'd expect from a model with stream power erosion and linear diffusion. \n",
+    "\n",
+    "Next we will create some volcanic deposits that fill the channels in our model."
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "volcano_deposits = np.zeros(mg2.size('node'))\n",
-    "da_big_enough = mg2['node']['drainage_area']>1e4\n",
-    "volcano_deposits[da_big_enough] = 30.0\n",
-    "volcano_deposits[mg2.boundary_nodes] = 0.0\n",
-    "rb2.add_layer(volcano_deposits, rock_id=1)\n",
+    "volcanic_deposits = np.zeros(mg2.size('node'))\n",
+    "da_big_enough = mg2['node']['drainage_area']>5e4\n",
     "\n",
-    "imshow_grid(mg2, da_big_enough)"
+    "topo_difference_from_top = mg2['node']['topographic__elevation'].max() - mg2['node']['topographic__elevation']\n",
+    "\n",
+    "volcanic_deposits[da_big_enough] = 0.5 * topo_difference_from_top[da_big_enough]\n",
+    "volcanic_deposits[mg2.boundary_nodes] = 0.0\n",
+    "\n",
+    "z2 += volcanic_deposits\n",
+    "rb2.run_one_step(rock_id=1)\n",
+    "\n",
+    "imshow_grid(mg2, volcanic_deposits)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "We should expect that the locations of our valleys and ridges change as the river system encouters the much stronger volcanic rock. "
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "for i in range(half_nts, nts):\n",
     "    fa2.run_one_step()\n",
-    "    sp2.run_one_step(dt = dt)\n",
-    "    dz_ad = np.zeros(mg2.size('node'))\n",
-    "    dz_ad[mg2.core_nodes] = U * dt\n",
-    "    z2 += dz_ad\n",
-    "    rb2.run_one_step(dz_advection = dz_ad)\n",
+    "    sp2.run_one_step(dt=dt)\n",
+    "    ld2.run_one_step(dt=dt)\n",
+    "    dz_ad2 = np.zeros(mg2.size('node'))\n",
+    "    dz_ad2[mg2.core_nodes] = U * dt\n",
+    "    z2 += dz_ad2\n",
+    "    rb2.run_one_step(dz_advection = dz_ad2, rock_id=0)\n",
     "    \n",
     "    for of in out_fields:\n",
-    "        ds2[of][:,:,i] = mg2['node'][of].reshape(mg2.shape)"
+    "        ds2[of][i,:,:] = mg2['node'][of].reshape(mg2.shape)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "Now that the model has run, let's plot the final elevation"
+   ]
   },
   {
    "cell_type": "code",
@@ -366,7 +390,9 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "And now a HoloView Plot that lets us explore the time evolution of the topography"
+   ]
   },
   {
    "cell_type": "code",
@@ -384,22 +410,24 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {
     "collapsed": true
    },
-   "outputs": [],
-   "source": []
+   "source": [
+    "Sure enough, the volcanic deposits impact the location of the ridges and valleys. The old valleys become ridges because it takes so much time for them to be eroded. \n",
+    "\n",
+    "You can explore how this changes as the thickness of the deposit changes and as the relative erodabilities change. \n",
+    "\n",
+    "\n",
+    "## The end.\n",
+    "\n",
+    "Nice work getting to the end of the tutorial!\n",
+    "\n",
+    "For more detailed information about the [RockBlock](http://landlab.readthedocs.io/en/release/landlab.layers.rockblock.html) and [LayeredRockBlock](http://landlab.readthedocs.io/en/release/landlab.layers.layeredrockblock.html) objects, check out their detailed documentation. \n",
+    "\n",
+    "# **Click [here](https://github.com/landlab/landlab/wiki/Tutorials) for more Landlab tutorials**"
+   ]
   }
  ],
  "metadata": {
