Note

This page was generated from a jupyter notebook.

Viewing Landlab output in ParaView

(GE Tucker, CU Boulder, June 2023)

ParaView is a powerful open-source visualization engine developed by kitware. It provides fast 3D rendering of surfaces and volumes, and is a great way to inspect output from a Landlab-built model. There are currently two ways to write Landlab grid and field data in a format that ParaView can read. The first, for raster grids, is to use the NetCDF file format. The second, and still under development as of this writing, is to use the Legacy VTK file format, which works for hex grid and in principle other grid types as well. This tutorial covers the first method: writing Landlab output to NetCDF file(s) for input into ParaView. It explains, with a worked example, how to create NetCDF output, how to read this into ParaView, how to convert the flat images into 3D surfaces, how to color them according to different fields, and how to create animations.

Creating a Landlab model as an example

In the example, we’ll create a simple landscape evolution model in which the landscape consists of a layer of soil on top of a layer of rock. Weathering of rock to soil will be calculated using an inverse exponential function, as implemented by the ExponentialWeatherer component. Soil creep will be represented with a nonlinear diffusion model, using the DepthDependentTaylorDiffuser component (Barnhart et al., 2019). Runoff routing will be modeled with a steepest descent “D8” algorithm using the FlowAccumulator component (Barnhart et al., 2020). Fluvial erosion, transport, and deposition will be modeled using the SPACE component (Shobe et al., 2017). To learn more about the theory, math, and numerics behind these, see the references listed.

Begin with some imports:

import numpy as np

from landlab import RasterModelGrid, imshow_grid
from landlab.components import (
    DepthDependentTaylorDiffuser,
    ExponentialWeatherer,
    FlowAccumulator,
    SpaceLargeScaleEroder,
)
from landlab.io.netcdf import write_raster_netcdf

# Parameters
nrows = 100  # number of node rows
ncols = 160  # number of node columns
dx = 10.0  # grid node spacing, m
max_soil_prod_rate = 0.001  # maximum soil production rate, m/y
soil_prod_decay_depth = 0.5  # decay depth for soil production, m
soil_transport_velocity = 0.02  # transport coefficient for soil creep, m/y
slope_crit = 1.0  # threshold slope factor for soil cree, -
soil_transport_decay_depth = 0.5  # decay depth for transport rate, m
nterms = 2  # number of terms for diffusion equation
K_sed = 0.001  # erosion coefficient for sediment
K_br = 0.0001  # erosion coefficient for rock
nsteps = 200  # number of time steps
dt = 10.0  # time-step duration, years
base_output_name = "eroding_landscape"  # base name for output files
output_interval = 100.0  # interval for output, y
fields_to_output = [
    "topographic__elevation",
    "soil__depth",
    "bedrock__elevation",
    "soil_production__rate",
    "drainage_area",
    "surface_water__discharge",
    "topographic__steepest_slope",
    "sediment__influx",
    "sediment__outflux",
]

# Create the grid
grid = RasterModelGrid((nrows, ncols), xy_spacing=dx)

# Create input fields
elev = grid.add_zeros("topographic__elevation", at="node")
soil = grid.add_zeros("soil__depth", at="node")
rock = grid.add_zeros("bedrock__elevation", at="node")

# Instantiate components
weatherer = ExponentialWeatherer(grid, max_soil_prod_rate, soil_prod_decay_depth)
diffuser = DepthDependentTaylorDiffuser(
    grid,
    soil_transport_velocity=soil_transport_velocity,
    slope_crit=slope_crit,
    soil_transport_decay_depth=soil_transport_decay_depth,
    nterms=nterms,
)
router = FlowAccumulator(grid, flow_director="FlowDirectorD8")
space = SpaceLargeScaleEroder(grid, K_sed=K_sed, K_br=K_br)

# Setup initial topography in the fashion of a couple of strath terraces
high_terrace = 20.0
low_terrace = 10.0
elev[:] = high_terrace
high_terr_cutout_amp = 380.0
high_terr_cutout_pd = 1300.0
low_terr_cutout_amp = 90.0
low_terr_cutout_pd = 700.0
init_soil_high_terr = 2.0
init_soil_low_terr = 1.0
init_soil_base = 0.5

# Make a sinusoidal cutout into the high terrace
trace_y = high_terr_cutout_amp * np.sin(
    2.0 * np.pi * grid.x_of_node / high_terr_cutout_pd
)
elev[grid.y_of_node < trace_y] = low_terrace
elev[grid.y_of_node > trace_y + 2 * high_terr_cutout_amp] = low_terrace

# ...and the low terrace
trace_y = low_terr_cutout_amp * np.sin(
    2.0 * np.pi * grid.x_of_node / low_terr_cutout_pd
)
elev[grid.y_of_node < trace_y + low_terr_cutout_amp] = 0.0
elev[grid.y_of_node > trace_y + 10 * low_terr_cutout_amp] = 0.0

# add some random noise and smooth it with a moving average
elev[:] += np.random.rand(grid.number_of_nodes)

soil[:] = init_soil_low_terr
soil[elev > low_terrace + 1.0] = init_soil_high_terr
soil[elev < low_terrace] = init_soil_base
rock[:] = elev - soil

imshow_grid(grid, elev)

imshow_grid(grid, soil)

# Setup for output, and write first file
router.run_one_step()  # to generate some flow accumulation for the 0th output
weatherer.calc_soil_prod_rate()  # ditto for soil prod rate
frame_number = 0
write_raster_netcdf(
    base_output_name + str(frame_number).zfill(4) + ".nc",
    grid,
    names=fields_to_output,
)
next_output = output_interval

# Run model
for i in range(1, nsteps + 1):
    router.run_one_step()
    weatherer.calc_soil_prod_rate()
    diffuser.run_one_step(dt)
    space.run_one_step(dt)
    if i * dt >= next_output:
        frame_number += 1
        write_raster_netcdf(
            base_output_name + str(frame_number).zfill(4) + ".nc",
            grid,
            names=fields_to_output,
        )
        next_output += output_interval

Installing ParaView

The first step is to download and install the ParaView application from https://www.paraview.org/download/.

Reading output files into ParaView

1. Open ParaView

2. Click the open-file icon in the upper left, or choose File -> Open from the menu bar

3. Choose NetCDF Reader

4. Click the blue Apply button

Inspecting different fields

The drop-down menu toward the upper left gives you a list of fields to choose from. Try selecting different ones.

Running animations

Look for the playback controls at the top of the window. Use the right arrow to play an animation.

Viewing in 3D

So far the images are just that: 2D images. To make the landscape 3D, we have to tell ParaView which field to use for the 3rd dimension, and then “extrude” it. The command to do this is actually a plugin “filter”:

1. From the menu bar, select Filters -> Alphabetical -> (scroll waaay down) Warp By Scalar

2. Around the center left of the main window, look for a pop-up menu called “Scalars”, and select “topographic__elevation”

3. Click the Apply button. You should notice that the terrain image now has some shading to it.

4. To view it in 3D, look for the little button just above the main view panel called “2D”. Click it and it will switch to saying “3D”, meaning you are now in 3D view mode.

You should now be able to use the mouse to rotate the 3D image and zoom in or out.

Closing thoughts

This tutorial just gives a small taste of what’s possible using ParaView. Check out their documentation to learn more. And see the upcoming tutorial on using Legacy VTK file output as an alternative way to get Landlab output into ParaView, which also works for Hex grids.

References

Barnhart, K. R., Glade, R. C., Shobe, C. M., and Tucker, G. E. (2019) Terrainbento 1.0: a Python package for multi-model analysis in long-term drainage basin evolution. Geosci. Model Dev., v. 12, p. 1267-1297, doi:10.5194/gmd-12-1267-2019.

Barnhart, K.R., Hutton, E.W.H., Tucker, G.E., Gasparini, N.M., Istanbulluoglu, E., Hobley, D.E.J., Lyons⁠, N.J., Mouchene, M., Nudurupati, S.S., Adams, J.M., and Bandaragoda, C. (2020) Short communication: Landlab 2.0: A software package for Earth surface dynamics. Earth Surface Dynamics, 8, 379–397, doi:10.5194/esurf-8-379-2020.

Shobe, C.M., Tucker, G.E., and Barnhart, K.R. The SPACE 1.0 model: a Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution. Geoscientific Model Development, 2017, doi:10.5194/gmd-10-4577-2017.

---

Generated by nbsphinx from a Jupyter notebook.