Terrain Following Mesh

Use a topographic surface to create a 3D terrain-following mesh.

Terrain following meshes are common in the environmental sciences, for instance in hydrological modelling (see Maxwell 2013 and ParFlow).

In this example, we domonstrate a simple way to make a 3D grid/mesh that follows a given topographic surface. In this example, it is important to note that the given digital elevation model (DEM) is structured (gridded and not triangulated): this is common for DEMs.

# sphinx_gallery_thumbnail_number = 3
import pyvista as pv
import numpy as np
from pyvista import examples

Download a gridded topography surface (DEM)

dem = examples.download_crater_topo()
HeaderData Arrays
N Cells1677401
N Points1680000
X Bounds1.810e+06, 1.831e+06
Y Bounds5.640e+06, 5.658e+06
Z Bounds0.000e+00, 0.000e+00
Dimensions1400, 1200, 1
Spacing1.500e+01, 1.500e+01, 0.000e+00
N Arrays1
NameFieldTypeN CompMinMax

Now let’s subsample and extract an area of interest to make this example simple (also the DEM we just load is pretty big). Since the DEM we loaded is a pyvista.UniformGrid mesh, we can use the pyvista.UniformGridFilters.extract_subset() filter:

subset = dem.extract_subset((500, 900, 400, 800, 0, 0), (5,5,1))
terrain mesh


[(1820500.0, 5649000.0, 16392.304845413266),
 (1820500.0, 5649000.0, 0.0),
 (0.0, 1.0, 0.0)]

Now that we have a region of interest for our terrain following mesh, lets make a 3D surface of that DEM:

terrain = subset.warp_by_scalar()
HeaderData Arrays
N Cells6400
N Points6561
X Bounds1.818e+06, 1.824e+06
Y Bounds5.646e+06, 5.652e+06
Z Bounds1.441e+03, 2.769e+03
Dimensions81, 81, 1
N Arrays1
NameFieldTypeN CompMinMax

terrain mesh


[(1830079.2876189426, 5658579.287618943, 11684.58760673553),
 (1820500.0, 5649000.0, 2105.2999877929688),
 (0.0, 0.0, 1.0)]

And now we have a 3D structured surface of the terrain! We can now extend that structured surface into a 3D mesh to form a terrain following grid. To do this, we first our cell spacings in the z-direction (these start from the terrain surface). Then we repeat the XYZ structured coordinates of the terrain mesh and decrease each Z level by our Z cell spacing. Once we have those structured coordinates, we can create a pyvista.StructuredGrid.

z_cells = np.array([25]*5 + [35]*3 + [50]*2 + [75, 100])

xx = np.repeat(terrain.x, len(z_cells), axis=-1)
yy = np.repeat(terrain.y, len(z_cells), axis=-1)
zz = np.repeat(terrain.z, len(z_cells), axis=-1) - np.cumsum(z_cells).reshape((1, 1, -1))

mesh = pv.StructuredGrid(xx, yy, zz)
mesh["Elevation"] = zz.ravel(order="F")
HeaderData Arrays
N Cells70400
N Points78732
X Bounds1.818e+06, 1.824e+06
Y Bounds5.646e+06, 5.652e+06
Z Bounds9.364e+02, 2.744e+03
Dimensions81, 81, 12
N Arrays1
NameFieldTypeN CompMinMax

cpos = [(1826736.796308761, 5655837.275274233, 4676.8405505181745),
 (1821066.1790519988, 5649248.765538796, 943.0995128226014),
 (-0.2797856225380979, -0.27966946337594883, 0.9184252809434081)]

mesh.plot(show_edges=True, lighting=False, cpos=cpos)
terrain mesh


[(1826736.796308761, 5655837.275274233, 4676.8405505181745),
 (1821066.1790519988, 5649248.765538796, 943.0995128226014),
 (-0.27978562253809786, -0.2796694633759488, 0.9184252809434079)]

Total running time of the script: ( 0 minutes 7.271 seconds)

Gallery generated by Sphinx-Gallery