# 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.

```import numpy as np

import pyvista as pv
from pyvista import examples
```

```dem = examples.download_crater_topo()
dem
```
UniformGridInformation
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
scalar1of1Pointsfloat6417.339e+022.787e+03

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))
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:

```terrain = subset.warp_by_scalar()
terrain
```
StructuredGridInformation
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
scalar1of1Pointsfloat6411.441e+032.769e+03

```terrain.plot()
``` 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( * 5 +  * 3 +  * 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
```
StructuredGridInformation
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
ElevationPointsfloat6419.364e+022.744e+03

```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)
``` Total running time of the script: ( 0 minutes 1.594 seconds)

Gallery generated by Sphinx-Gallery