{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Terrain Following Mesh {#terrain_following_mesh_example}\n\nUse a topographic surface to create a 3D terrain-following mesh.\n\nTerrain following meshes are common in the environmental sciences, for\ninstance in hydrological modelling (see [Maxwell\n2013](https://www.sciencedirect.com/science/article/abs/pii/S0309170812002564)\nand [ParFlow](https://parflow.org)).\n\nIn this example, we demonstrate a simple way to make a 3D grid/mesh that\nfollows a given topographic surface. In this example, it is important to\nnote that the given digital elevation model (DEM) is structured (gridded\nand not triangulated): this is common for DEMs.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n\nimport pyvista as pv\nfrom pyvista import examples"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Download a gridded topography surface (DEM)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dem = examples.download_crater_topo()\ndem"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now let\\'s subsample and extract an area of interest to make this\nexample simple (also the DEM we just load is pretty big). Since the DEM\nwe loaded is a `pyvista.ImageData`{.interpreted-text role=\"class\"} mesh,\nwe can use the\n`pyvista.ImageDataFilters.extract_subset`{.interpreted-text role=\"func\"}\nfilter:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "subset = dem.extract_subset((500, 900, 400, 800, 0, 0), (5, 5, 1))\nsubset.plot(cpos=\"xy\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now that we have a region of interest for our terrain following mesh,\nlets make a 3D surface of that DEM:\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "terrain = subset.warp_by_scalar()\nterrain"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "terrain.plot()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And now we have a 3D structured surface of the terrain. We can now\nextend that structured surface into a 3D mesh to form a terrain\nfollowing grid. To do this, we first our cell spacings in the\nz-direction (these start from the terrain surface). Then we repeat the\nXYZ structured coordinates of the terrain mesh and decrease each Z level\nby our Z cell spacing. Once we have those structured coordinates, we can\ncreate a `pyvista.StructuredGrid`{.interpreted-text role=\"class\"}.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "z_cells = np.array([25] * 5 + [35] * 3 + [50] * 2 + [75, 100])\n\nxx = np.repeat(terrain.x, len(z_cells), axis=-1)\nyy = np.repeat(terrain.y, len(z_cells), axis=-1)\nzz = np.repeat(terrain.z, len(z_cells), axis=-1) - np.cumsum(z_cells).reshape((1, 1, -1))\n\nmesh = pv.StructuredGrid(xx, yy, zz)\nmesh[\"Elevation\"] = zz.ravel(order=\"F\")\nmesh"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "cpos = [\n    (1826736.796308761, 5655837.275274233, 4676.8405505181745),\n    (1821066.1790519988, 5649248.765538796, 943.0995128226014),\n    (-0.2797856225380979, -0.27966946337594883, 0.9184252809434081),\n]\n\nmesh.plot(show_edges=True, lighting=False, cpos=cpos)"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.2"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}