Tutorial that demonstrates how to use a landlab `LayerStack`.

layerstack_tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from landlab import RasterModelGrid\n",
    "from landlab.layers import LayerStack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creating a stack"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a single `LayerStack` with a particular base elevation, `z0`. This is the elevation of the bottom of the stack onto which sediment will be added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stack = LayerStack(z0=2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stack as a `base` and `top` property that gives the elevation to the bottom and top of the stack. Since it doesn't have any sediment in it yet, its base is the same elevation as its top. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.0, 2.0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.base, stack.top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sediment is added with the `add` method. Sediment layers are recorded in increments of `dz` and elevations to layer boundaries with the `z` property. The default `dz` is 1.0 (this can be changed with the `dz` keyword when instantiating a `LayerStack`) and, since there are no sediment layers yet, the `z` array records only the boundary at the base of the stack. Notice that `z` is elevation relative to the base of the stack, not its absolute elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.0, array([ 0.]))"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.dz, stack.z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding and removing sediment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll add some sediment. Enough to fill a single \"bin\" (or \"basket\")."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stack.add(1.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Having added enough sediment to fill 1.25 bins, we now have a partly filled bin on top of a completely filled bin. Thus, we have two layers and three boundaries - bottom of the first layer, top of the first layer (or bottom of the second layer) and the top of the second layer. The length of the `z` array will always be one more than the number of layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, array([ 0.  ,  1.  ,  1.25]))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.nlayers, stack.z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that the stack contains sediment the top of the stack is higher than the base by the amount of sediment that we added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.0, 3.25)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.base, stack.top"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `remove` method to remove sediment from the top of a stack."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "stack.remove(.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is now only one (partially filled) layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, array([ 0.  ,  0.75]))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.nlayers, stack.z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Moving the stack up and down"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stack can be moved up and down in two different ways. The first way is to change the value of its `base`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.5"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.base += .5\n",
    "stack.base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The second way is to change the value of its top."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.75"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.top = -1.\n",
    "stack.base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also some convience methods that move the stack up and down by an increment. Move the stack up with the `lift` method..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.75"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.lift(1.)\n",
    "stack.base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and down with the lower method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.25"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.lower(.5)\n",
    "stack.base"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting layers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `layer_at` method to get the layer that contains a particular elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.  ,  1.  ,  2.  ,  3.  ,  4.  ,  5.  ,  5.75])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.add(5.)\n",
    "stack.z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.layer_at(1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If an elevation is at a boundary, the default is to use the upper layer. You can use the `upper` keyword to control this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1, 0)"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.layer_at(1.), stack.layer_at(1., lower=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To extract the layer boundaries for a portion of a column, use the `extract` method. By default this methods gets all layers above a particular elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.75,  1.  ,  2.  ,  3.  ,  4.  ,  4.5 ])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.extract(.75)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If a second argument is provided, it gives the stop elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.75,  1.  ,  2.  ,  3.  ,  3.25])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.extract(.75, 3.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Memory managment"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Layer elevations (and other sediment properties, as we'll see later) are stored as numpy arrays. This means that as sediment is added to stacks, the arrays must be resized. Since resizing arrays requires allocating additional memory as well as copying data, this is a slow operation. To help with this, a `LayerStack` allocates more memory than is needed to store the current number of layers so that more memory needn't be allocated with the addition of every layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stack = LayerStack()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The size of the arrays can be obtained through the `allocated` property. An empty stack has allocated enough memory to hold four layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.allocated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we now add enough sediment to the stack that it would extend it beyond the number of buffered layers, more memory is allocated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(6, 13)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.add(5)\n",
    "stack.z.size, stack.allocated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adding even more sediment results in a larger buffer. This buffer is approximately 1/8 larger than the actual number of layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(506, 576)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.add(500)\n",
    "stack.z.size, stack.allocated"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding sediment properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is possible to keep track of sediment properties along with the position of layers. To do this, use the `fields` keyword when the `LayerStack` is instantiated. The following tells the `LayerStack` to keep track of a property called *age*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stack = LayerStack(fields=('age', ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The new object now has a property, `age`, that is an array that records the value of that property for each layer. Since there is no sediment in the new stack, it's just an empty array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([], dtype=float64)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.age"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To add new sediment to the stack, use the `add` method but now pass the *age* of the new sediment as a keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.add(.75, age=1.)\n",
    "stack.age"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When a layer is composed of sediment with different properties, the value for that layer is a weighted (by the amount of each sediment in the layer) average of the two types."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.25,  2.  ])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.add(.5, age=2.)\n",
    "stack.age"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding multiple grain types"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus far, our sediment contains only a single grain type. To tell the `LayerStack` to keep track of multiple grain types, use the `n_grains` keyword."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "stack = LayerStack(n_grains=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The stack now has an `f` property, which is the fraction of each grain type in the layer. Since the stack is empty, so is the `f` array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([], shape=(0, 2), dtype=float64)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Add more sediment to the stack with some composition of grain types. Notice that the `f` keyword is an array of length `n_grains - 1`. This is because the fractions necessarily add to 1., and so specify the fraction of the third grain type would be redundant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.2,  0.3]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.add(.75, f=(.2, .3))\n",
    "stack.f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As with other sediment properties, the recorded fractions within a single layer are calculated as a weighted mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.275,  0.3  ],\n",
       "       [ 0.5  ,  0.3  ],\n",
       "       [ 0.5  ,  0.3  ]])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "stack.add(1.5, f=(.5, .3))\n",
    "stack.f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following gets the fractions of the third grain type for each stack of the layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.425,  0.2  ,  0.2  ])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1. - stack.f.sum(axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Layers in ModelGrid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grid = RasterModelGrid((3, 4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are multiple ways we could attach layers stacks to a `ModelGrid`. Here I've just added them as a field called *layers*. It's nice because it's easy, doesn't require any changes to `ModelGrid`, and makes it clear where the layers are defined. Another way would be to add it as a propery so they could be accessed as `grid.layers` - short and sweet. I'm sure there are other ways."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "grid.at_cell['layers'] = [LayerStack() for _ in range(grid.number_of_cells)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "layers = grid.at_cell['layers']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "layers[1].add(2.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0. ,  1. ,  2. ,  2.5])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "layers[1].z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}