{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# What happens when you create a grid object?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " For instructions on how to run an interactive iPython notebook, click here: https://github.com/landlab/tutorials/blob/release/README.md
\n", "For the unexpanded version to download and run, click here: https://nbviewer.jupyter.org/github/landlab/tutorials/blob/release/grid_object_demo/grid_object_demo_unexpanded.ipynb
\n", "For more Landlab tutorials, click here: https://github.com/landlab/landlab/wiki/Tutorials\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Landlab supports a range of grid types. These include both rasters (with both square and rectangular cells), and a range of structured and unstructured grids based around the interlocking polygons and triangles of a Voronoi-Delaunay tesselation (radial, hexagonal, and irregular grids).\n", "\n", "Here, we look at some of the features of both of these types.\n", "\n", "We can create **grid** objects with the following lines of code." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from landlab import RasterModelGrid, VoronoiDelaunayGrid, HexModelGrid\n", "\n", "smg = RasterModelGrid((3, 4), 1.) # a square-cell raster, 3 rows x 4 columns, unit spacing\n", "rmg = RasterModelGrid((3, 4), xy_spacing=(1., 2.)) # a rectangular-cell raster\n", "hmg = HexModelGrid(shape=(3, 4), dx=1.)\n", "# ^a hexagonal grid with 3 rows, 4 columns from the base row, & node spacing of 1.\n", "x = np.random.rand(100) * 100.\n", "y = np.random.rand(100) * 100.\n", "vmg = VoronoiDelaunayGrid(x, y)\n", "# ^a Voronoi-cell grid with 100 randomly positioned nodes within a 100.x100. square" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "All these various `ModelGrid` objects contains various data items (known as *attributes*). These include, for example:\n", "* number nodes and links in the grid\n", "* *x* and *y* coordinates of each each node\n", "* starting (\"tail\") and ending (\"head\") node IDs of each link\n", "* IDs of links that are active\n", "* IDs of core nodes\n", "* etc.\n", "\n", "From here on we'll focus on the square raster grid as its geometry is a bit easier to think through, but all of the following applies to all grid types.\n", "\n", "## Understanding the topology of Landlab grids\n", "\n", "All grids consist of two interlocked sets of *points* joined by *lines* outlining *areas*. If we define data on the points we call **nodes**, then they are joined by **links**, which outline **patches**. Each node within the interior of the grid lies at the geometric center of the area of a **cell**. The cell's edges are **faces**, and the endpoints of the faces---which are also vertices of the cells---are **corners**.\n", "\n", "Note that this kind of scheme requires one set of features to be \"dominant\" over the other; i.e., either not every node has a cell, *or* not every link is crossed by a face. Both cannot be true, because one or other set of features has to define the edge of the grid. Landlab assumes that the node set is primary, so there are always more nodes than corners; more links than faces; and more patches than cells.\n", "\n", "Each of these sets of *\"elements\"* has its own set of IDs. These IDs are what allow us to index the various Landlab fields, which store spatial data. Each feature is ordered by **x, then y**. The origin is always at the bottom left node, unless you choose to move it (`grid.move_origin`)... except in the specific case of a radial grid, where logic and symmetry dictates it must be the central node.\n", "\n", "Whenever Landlab needs to order something rotationally (angles; elements around a different element type), it does so following the standard mathematical convention of **counterclockwise from east**. We'll see this in practical terms a bit later in this tutorial.\n", "\n", "The final thing to know is that **links and faces have directions**. This lets us record fluxes on the grid by associating them with, and mapping them onto, the links (or, much less commonly, the faces). All lines point into the **upper right half-space**. So, on our raster, this means the horizontal links point east and the vertical links point north.\n", "\n", "So, for reference, our raster grid looks like this:\n", "\n", "\n", " NODES: LINKS: PATCHES:\n", " 8 ----- 9 ---- 10 ---- 11 * -14-->* -15-->* -16-->* * ----- * ----- * ----- *\n", " | | | | ^ ^ ^ ^ | | | |\n", " | | | | 10 11 12 13 | 3 | 4 | 5 |\n", " | | | | | | | | | | | |\n", " 4 ----- 5 ----- 6 ----- 7 * --7-- * --8-- * --9-- * * ----- * ----- * ----- *\n", " | | | | ^ ^ ^ ^ | | | |\n", " | | | | 3 4 5 6 | 0 | 1 | 2 |\n", " | | | | | | | | | | | |\n", " 0 ----- 1 ----- 2 ----- 3 * --0-->* --1-->* --2-->* * ----- * ----- * ----- *\n", "\n", " CELLS: FACES: CORNERS:\n", " * ----- * ----- * ----- * * ----- * ----- * ----- * * ----- * ----- * ----- *\n", " | | | | | | | | | | | |\n", " | . ----- . ----- . | | . --5-->. --6-->. | | 3 ----- 4 ----- 5 |\n", " | | | | | | ^ ^ ^ | | | | | |\n", " * --| 0 | 1 |-- * * --2 3 4-- * * --| | |-- *\n", " | | | | | | | | | | | | | | |\n", " | . ----- . ----- . | | . --0-->. --1-->. | | 0 ----- 1 ----- 2 |\n", " | | | | | | | | | | | |\n", " * ----- * ----- * ----- * * ----- * ----- * ----- * * ----- * ----- * ----- *\n", "\n", "\n", "## Recording and indexing the values at elements\n", "\n", "Landlab lets you record values at any element you want. In practice, the most useful places to store data is on the primary elements of nodes, links, and patches, with the nodes being most useful for scalar values (e.g, elevations) and the links for fluxes with direction to them (e.g., velocity or discharge).\n", "\n", "In order to maintain compatibility across data types, *all* landlab data are stored in *number-of-elements*-long arrays. This includes both user-defined data and the properties of the nodes within the grid. This means that these arrays can be immediately indexed by their element ID. For example:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1., 1.])" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# what are the y-coordinates of the pair of nodes in the middle of our 3-by-4 grid?\n", "# the IDs of these nodes are 5 and 6, so:\n", "smg.y_of_node[[5, 6]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you're working with a raster, you can always reshape the value arrays back into two dimensions so you can take Numpy-style slices through it:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0., 1., 2., 3.])" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# what are the x-coordinates of nodes in the middle row?\n", "smg.x_of_node.reshape(smg.shape)[1, :]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This same data storage pattern is what underlies the Landlab **data fields**, which are simply one dimensional, number-of-elements-long arrays that store user defined spatial data across the grid, attached to the grid itself." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.add_zeros('node', 'elevation', noclobber=False)\n", "# ^Creates a new field of zero data associated with nodes\n", "smg.at_node['elevation'] # Note the use of dictionary syntax" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Or, equivalently, at links:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,\n", " 1., 1., 1., 1.])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.add_ones('link', 'slope', noclobber=False)\n", "# ^Creates a new array of data associated with links\n", "smg.at_link['slope']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The Landlab **components** use fields to share spatial information among themselves. See the *fields* and *components* tutorials for more information.\n", "\n", "\n", "## Getting this information from the grid object\n", "\n", "All of this topological information is recorded within our grid objects, and can be used to work with data arrays that are defined over the grid. The grid records the numbers of each element, their positions, and their relationships with one another. Let's take a look at some of this information for the raster:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "12" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.number_of_nodes" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "17" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.number_of_links" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The grid contains its geometric information too. Let's look at the *(x,y)* coordinates of the nodes:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0 0.0 0.0\n", "1 1.0 0.0\n", "2 2.0 0.0\n", "3 3.0 0.0\n", "4 0.0 1.0\n", "5 1.0 1.0\n", "6 2.0 1.0\n", "7 3.0 1.0\n", "8 0.0 2.0\n", "9 1.0 2.0\n", "10 2.0 2.0\n", "11 3.0 2.0\n" ] } ], "source": [ "for i in range(smg.number_of_nodes):\n", " print(i, smg.x_of_node[i], smg.y_of_node[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Link connectivity and direction is described by specifying the starting (\"tail\") and ending (\"head\") node IDs for each link (to remember this, think of an arrow: TAIL ===> HEAD)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Link 0 : node 0 ===> node 1\n", "Link 1 : node 1 ===> node 2\n", "Link 2 : node 2 ===> node 3\n", "Link 3 : node 0 ===> node 4\n", "Link 4 : node 1 ===> node 5\n", "Link 5 : node 2 ===> node 6\n", "Link 6 : node 3 ===> node 7\n", "Link 7 : node 4 ===> node 5\n", "Link 8 : node 5 ===> node 6\n", "Link 9 : node 6 ===> node 7\n", "Link 10 : node 4 ===> node 8\n", "Link 11 : node 5 ===> node 9\n", "Link 12 : node 6 ===> node 10\n", "Link 13 : node 7 ===> node 11\n", "Link 14 : node 8 ===> node 9\n", "Link 15 : node 9 ===> node 10\n", "Link 16 : node 10 ===> node 11\n" ] } ], "source": [ "for i in range(smg.number_of_links):\n", " print('Link',i,': node',smg.node_at_link_tail[i],'===> node',smg.node_at_link_head[i])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Boundary conditions are likewise defined on these elements (see also the full boundary conditions tutorial). Landlab is clever enough to ensure that the boundary conditions recorded on, say, the links get updated when you redefine the conditions on, say, the nodes.\n", "\n", "Nodes can be *core*, *fixed value*, *fixed gradient*, or *closed* (flux into or out of node is forbidden). Links can be *active* (can carry flux), *fixed* (always carries the same flux; joined to a fixed gradient node) or *inactive* (forbidden from carrying flux). \n", "\n", "Note that this boundary coding does not mean that a particular boundary condition is automatically enforced. It's up to the user to take advantage of these codes. For example, if you are writing a model that calculates flow velocity on links but wish the velocity to be zero at inactive links, you the programmer must ensure this, for instance by including a line like `my_velocity[grid.inactive_links] = 0.0`, or alternatively `my_velocity[grid.active_links] = ......`.\n", "\n", "Information on boundary coding is available from the grid:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5, 6])" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.core_nodes" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 4, 5, 7, 8, 9, 11, 12])" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.active_links" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 7, 8, 9, 11, 12])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# let's demonstrate the auto-updating of boundary conditions:\n", "from landlab import CLOSED_BOUNDARY\n", "smg.status_at_node[smg.nodes_at_bottom_edge] = CLOSED_BOUNDARY\n", "smg.active_links # the links connected to the bottom edge nodes are now inactive" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Element connectivity\n", "\n", "Importantly, we can also find out which elements are connected to which other elements. This allows us to do computationally vital operations involving mapping values defined at one element onto another, e.g., the net flux at a node; the mean slope at a patch; the node value at a cell.\n", "\n", "In cases where these relationships are one-to-many (e.g., `links_at_node`, `nodes_at_patch`), the shape of the resulting arrays is always (number_of_elements, max-number-of-connected-elements-across-grid). For example, on a raster, `links_at_node` is (nnodes, 4), because the cells are always square. On an irregular Voronoi-cell grid, `links_at_node` will be (nnodes, X) where X is the number of sides of the side-iest cell, and `nodes_at_patch` will be (npatches, 3) because all the patches are Delaunay triangles. And so on.\n", "\n", "Lets take a look. Remember, Landlab orders things **clockwise from east**, so for a raster the order will the EAST, NORTH, WEST, SOUTH." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 8, 11, 7, 4])" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.links_at_node[5]" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(12, 4)" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.links_at_node.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Undefined directions get recorded as `-1`:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([14, -1, -1, 10])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.links_at_node[8]" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, -1, -1, -1],\n", " [ 1, 0, -1, -1],\n", " [ 2, 1, -1, -1],\n", " [-1, 2, -1, -1],\n", " [ 3, -1, -1, 0],\n", " [ 4, 3, 0, 1],\n", " [ 5, 4, 1, 2],\n", " [-1, 5, 2, -1],\n", " [-1, -1, -1, 3],\n", " [-1, -1, 3, 4],\n", " [-1, -1, 4, 5],\n", " [-1, -1, 5, -1]])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.patches_at_node" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 5, 4, 0, 1],\n", " [ 6, 5, 1, 2],\n", " [ 7, 6, 2, 3],\n", " [ 9, 8, 4, 5],\n", " [10, 9, 5, 6],\n", " [11, 10, 6, 7]])" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.nodes_at_patch" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Where element-to-element mapping is one-to-one, you get simple, one dimensional arrays:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([5, 6])" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.node_at_cell # shape is (n_cells, )" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([-1, -1, -1, -1, -1, 0, 1, -1, -1, -1, -1, -1])" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.cell_at_node # shape is (n_nodes, ) with -1s as needed" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A bit of thought reveals that things get more complicated for links and faces, because they have direction. You'll need a convenient way to record whether a given flux (which is positive if it goes with the link's inherent direction, and negative if against) actually is travelling into or out of a given node. The grid provides `link_dirs_at_node` and `active_link_dirs_at_node` to help with this:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[-1, -1, 0, 0],\n", " [-1, -1, 1, 0],\n", " [-1, -1, 1, 0],\n", " [ 0, -1, 1, 0],\n", " [-1, -1, 0, 1],\n", " [-1, -1, 1, 1],\n", " [-1, -1, 1, 1],\n", " [ 0, -1, 1, 1],\n", " [-1, 0, 0, 1],\n", " [-1, 0, 1, 1],\n", " [-1, 0, 1, 1],\n", " [ 0, 0, 1, 1]], dtype=int8)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.link_dirs_at_node # all links; positive points INTO the node; zero where no link" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "True" ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# prove there are zeros where links are missing:\n", "np.all((smg.link_dirs_at_node == 0) == (smg.links_at_node == -1))" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[ 0, 0, 0, 0],\n", " [ 0, 0, 0, 0],\n", " [ 0, 0, 0, 0],\n", " [ 0, 0, 0, 0],\n", " [-1, 0, 0, 0],\n", " [-1, -1, 1, 0],\n", " [-1, -1, 1, 0],\n", " [ 0, 0, 1, 0],\n", " [ 0, 0, 0, 0],\n", " [ 0, 0, 0, 1],\n", " [ 0, 0, 0, 1],\n", " [ 0, 0, 0, 0]], dtype=int8)" ] }, "execution_count": 22, "metadata": {}, "output_type": "execute_result" } ], "source": [ "smg.active_link_dirs_at_node # in this one, inactive links get zero too" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Multiply the fluxes indexed by `links_at_node` and sum by axis=1 to have a very convenient way to calculate flux divergences at nodes:" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[-1. -1.]\n" ] } ], "source": [ "fluxes_at_node = smg.at_link['slope'][smg.links_at_node]\n", "# ^...remember we defined the slope field as ones, above\n", "fluxes_into_node = fluxes_at_node * smg.active_link_dirs_at_node\n", "flux_div_at_node = fluxes_into_node.sum(axis=1)\n", "print(flux_div_at_node[smg.core_nodes])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Why? Remember that earlier in this tutorial we already set the bottom edge to `CLOSED_BOUNDARY`. So each of our core nodes has a flux of +1.0 coming in from the left, but two fluxes of -1.0 leaving from both the top and the right. Hence, the flux divergence is -1. at each node.\n", "\n", "Note as well that Landlab offers the one-line grid method `calc_flux_div_at_node()` to perform this same operation. For more on this, see the **gradient_and_divergence** tutorial." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Click here for more Landlab tutorials" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.7" } }, "nbformat": 4, "nbformat_minor": 1 }