Note
Go to the end to download the full example code.
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 demonstrate 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.
from __future__ import annotations
import numpy as np
import pyvista as pv
from pyvista import examples
Download a gridded topography surface (DEM)
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.ImageData
mesh, we can use
the pyvista.ImageDataFilters.extract_subset()
filter:
subset = dem.extract_subset((500, 900, 400, 800, 0, 0), (5, 5, 1))
subset.plot(cpos='xy')

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

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')
mesh

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