partway through adding comments

[skip ci]

Author: kbarnhart

rockblock/introduction_to_rockblock.ipynb

@@ -4,7 +4,13 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Introduction to the RockBlock and LayeredRockBlock objects\n"
+    "# 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",
+    "\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).  "
    ]
   },
   {
@@ -31,7 +37,21 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "## Part 1: Eroding layered rock\n",
+    "\n",
+    "First we will create an instance of a LayedRockBlock. Both LayeredRockBlock and RockBlock work closely with a landlab ModelGrid and stores information about rock type at each node. \n",
+    "\n",
+    "To create a LayeredRock you need the following information:\n",
+    "\n",
+    "1. A dictionary of rock property attributes. \n",
+    "2. A model grid that has the field `'topographic__elevation'` already created.\n",
+    "3. A list of elevations that your layers will go through at specified anchor point (default value for the anchor point is (0, 0). \n",
+    "4. A functional form in two variables (x and y) that defines the shape of your surface. \n",
+    "\n",
+    "The use of this function form makes it possible for any function of x and y to be passed to Layered Rock Block.\n",
+    "In this example the option for an anticline is currently uncommented, though options for steep and shallow dipping layers are also provided. "
+   ]
   },
   {
    "cell_type": "code",
@@ -54,7 +74,7 @@
     "ids = np.tile([0,1], 20) \n",
     "\n",
     "# Anticline\n",
-    "rb = LayeredRockBlock(mg, z0s, ids, x0=5000, y0=5000, function = lambda x, y : ((0.001*x)**2+(0.003*y)**2), attrs=attrs)\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",
     "\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",
@@ -66,14 +86,16 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "Now that we've created the LayeredRockBlock, model grid fields for each of the LayeredRockBlock attributes exist and have been set to the values of the rock exposed at the surface. \n",
+    "\n",
+    "Here we plot the value of `'K_sp'` as a function of the model grid. "
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "imshow_grid(mg, 'K_sp')"
@@ -82,7 +104,15 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "As you can see (in the default anticline option) we have concentric elipses of stronger and weaker rock. \n",
+    "\n",
+    "Next, lets instantiate a FlowAccumulator and a FastscapeEroder to create a simple landscape evolution model. \n",
+    "\n",
+    "We will point the FastscapeEroder to the model grid field `'K_sp'` so that it will respond to the spatially variable erodabilities created by the LayeredRockBlock. \n",
+    "\n",
+    "We will also instatiate an xarray dataset used to store the output of our model through time for visualization. "
+   ]
   },
   {
    "cell_type": "code",
@@ -108,13 +138,17 @@
     "                           'rock_type__id': (('y', 'x', 't'), np.empty((mg.shape[0], mg.shape[1], nts)))},\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)})"
+    "                        't': dt*np.arange(nts)/1e6})"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "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",
+    "Within each timestep we save information about the model for plotting. "
+   ]
   },
   {
    "cell_type": "code",
@@ -139,30 +173,34 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "Now that the model has run, lets start by plotting the resulting topography. "
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
-    "imshow_grid(mg, 'topographic__elevation')"
+    "imshow_grid(mg, 'topographic__elevation', cmap='viridis')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "The layers of rock clearly influence the form of topography. \n",
+    "\n",
+    "Next we will use HoloViews to visualize the topography and rock type together. \n",
+    "\n",
+    "To start, we create a HoloViewDataset from our xarray datastructure. "
+   ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {
-    "collapsed": true
-   },
+   "metadata": {},
    "outputs": [],
    "source": [
     "hvds = hv.Dataset(ds)\n",
@@ -172,7 +210,9 @@
   {
    "cell_type": "markdown",
    "metadata": {},
-   "source": []
+   "source": [
+    "Next we specify that we want two images, one showing rock type and one showing topographic elevation. A slider bar "
+   ]
   },
   {
    "cell_type": "code",
@@ -236,7 +276,7 @@
     "                            '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)})\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",