Note

This page was generated from a jupyter notebook.

Plotting grid data with Landlab

This tutorial illustrates how you can plot spatial data in Landlab, focusing in particular on Landlab’s imshow_grid plotter and associated functions. Landlab’s plotters are built onto the widely used Matplotlib Python package.

We start by importing the NumPy library, which we’ll use in producing some data to plot:

import numpy as np

From Landlab, we’ll need a grid on which to plot data, and a plotting function. We’ll start with just imshow_grid, but be aware that similar but more specifically named functions like imshow_grid_at_node are also available. These all wrap the same basic Landlab functionality, so we’re taking the most general method.

Note that you can use imshow_grid as a function or as a method of any landlab grid. For example, the following two usages can be used interchangably,

>>> grid.imshow(values)

>>> imshow_grid(grid, values)

from landlab import RadialModelGrid, RasterModelGrid

We’ll also need some functions from matplotlib proper to help us handle our graphical output:

import matplotlib.pyplot as plt

Plotting in 2D

The imshow plotter method is Landlab’s primary function for plotting data distributed across the grid. It’s pretty powerful, and comes with a fairly extensive suite of options to control the appearance of your output. You can see the full list of options in the imshow_grid documentation.

However, most simply, it just takes grid.imshow(data). Data can be either a field name string, or an array of the data itself.

%matplotlib inline
rmg = RasterModelGrid((50, 50), 1.0)
rmg.imshow(rmg.x_of_node)  # plot the x distances at nodes
plt.show()

Those units for the axis are taken from the grid property axis_units, which is a tuple that we can set. Alternatively, pass a tuple directly to the plotter with the keyword grid_units.

While we’re at it, let’s plot from a field instead of an array, and also mix up the default color scheme. The cmap keyword can take any input that you could also supply to matplotlib; see, e.g., http://matplotlib.org/examples/color/colormaps_reference.html.

rmg.axis_units = ("km", "km")
_ = rmg.add_field(
    "myfield", (rmg.x_of_node**2 + rmg.y_of_node**2) ** 0.5, at="node", clobber=True
)
rmg.imshow("myfield", cmap="bone")
plt.show()

The plotter works just fine with both raster grids and irregular grids. Name a plot with the var_name keyword.

radmg = RadialModelGrid(n_rings=10, spacing=10.0)
plt.subplot(121)
rmg.imshow(rmg.y_of_node, allow_colorbar=False, plot_name="regular grid")
plt.subplot(122)
radmg.imshow(radmg.x_of_node, allow_colorbar=False, plot_name="irregular grid")
plt.show()

Now, let’s look at some of the other more advanced options imshow can provide.

imshow offers plenty of keyword options for modifying the colorbar, including var_name, var_units, symmetric_cbar, vmin, vmax, and shrink. We’ve already seen allow_colorbar, which lets you suppress the bar entirely. Let’s see some in action.

radz = (radmg.x_of_node**2 + radmg.y_of_node**2) ** 0.5
radz = radz.max() - radz - 0.75 * radz.mean()
# let's plot these elevations truncated at radz >= 0
radmg.imshow(
    radz,
    grid_units=("m", "m"),
    vmin=0.0,
    shrink=0.75,
    var_name="radz",
    var_units="no units",
)
plt.show()

Now let’s explore color control. The grid takes the keyword color_for_background, which as you’d expect, colors any exposed part of the frame without cells over it. It knows the same color representations as matplotlib, e.g., (0., 0., 0.5), ‘0.3’, ‘b’, ‘yellow’.

radmg.imshow(radmg.y_of_node, color_for_background="0.3")
plt.show()

The plotter knows about boundary condition status, and we can control the colour of such nodes as well. This is useful if plotting an irregular watershed on a raster, for example. Here, None means transparent, as we will see in the next example.

rmg2 = RasterModelGrid((50, 50), (1.0, 2.0))
myvals = ((rmg2.x_of_node - 50.0) ** 2 + (rmg2.y_of_node - 25.0) ** 2) ** 0.5
rmg2.status_at_node[myvals > 30.0] = rmg2.BC_NODE_IS_CLOSED
rmg2.imshow(myvals, color_for_closed="blue", shrink=0.6)

Finally, note that the plotter recognises any masked node in a masked array as a closed node. This can be used as a convenient way to make grid overlays, as follows:

mymask_1stcondition = np.logical_or(rmg.x_of_node < 15, rmg.x_of_node > 35)
mymask_2ndcondition = np.logical_or(rmg.y_of_node < 15, rmg.y_of_node > 35)
mymask = np.logical_or(mymask_1stcondition, mymask_2ndcondition)
overlay_data = np.ma.array(rmg.y_of_node, mask=mymask, copy=False)
rmg.imshow(rmg.x_of_node)
rmg.imshow(overlay_data, color_for_closed=None, cmap="winter")
plt.show()

Plotting in 1D

Landlab basically lets you get on with it for yourself if plotting cross sections, or otherwise in 1D. We recommend the basic matplotlib plotting suite. Often plot() is totally adequate.

For a simple grid cross section, just reshape the data array back to a raster and take a slice:

# Note, in Landlab 2.0 a new component that will permit profiles based on endpoints
# will be added to the component library.

mg = RasterModelGrid((30, 30))
z = (mg.x_of_node**2 + mg.y_of_node**2) ** 0.5
z = z.max() - z
z_raster = z.reshape(mg.shape)
x_raster = mg.x_of_node.reshape(mg.shape)
for i in range(0, 30, 5):
    plt.plot(x_raster[i, :], z_raster[i, :])
plt.title("east-west cross sections though z")
plt.xlabel("x (m)")
plt.ylabel("z (m)")
plt.show()

Additionally, Landlab makes available a stream profiler tool. It finds the highest drainage area node in a landscape whenever it’s called, then follows the drainage structure back upstream from that node, always choosing the upstream node with the highest drainage area. This means we can do things like this:

from landlab.components import ChannelProfiler, FastscapeEroder, FlowAccumulator

mg = RasterModelGrid((100, 100), 1000.0)
mg.axis_units = ("m", "m")
z = mg.add_zeros("topographic__elevation", at="node")
z += np.random.rand(mg.number_of_nodes)  # roughen the initial surface

fr = FlowAccumulator(mg)
sp = FastscapeEroder(mg, K_sp=1.0e-5)
dt = 50000.0

for ndt in range(100):
    z[mg.core_nodes] += 10.0
    fr.run_one_step()
    sp.run_one_step(dt)
    if ndt % 5 == 0:
        print(ndt)

prf = ChannelProfiler(
    mg, number_of_watersheds=4, main_channel_only=False, minimum_channel_threshold=1e7
)
prf.run_one_step()

plt.figure(1)
prf.plot_profiles()
plt.show()

plt.figure(1)
prf.plot_profiles_in_map_view()
plt.show()

---

Generated by nbsphinx from a Jupyter notebook.