1    Terrain Visualization

In this activity, you will work with creating and manipulating 3D surface meshes using PyVista, a Python interface for the Visualization Toolkit  (VTK). VTK is a powerful open-source library for computer graphics, visualization, and image processing.  You can learn more about both tools through these references: - https://docs.pyvista.org/ - https://vtk.org/

We will also be using the itkwidgets library, which provides interactive Jupyter widgets for plot- ting, to visualize our meshes.

The outline of this activity will be:  1.  Creating a 3D surface mesh 2. Writing code to coarsen the mesh 3.  Writing code to visualize the error in elevation between the original mesh and the coarse mesh

2    Part 1:  Creating a 3D Surface Mesh

We will start by using 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).

Below, 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.

[1]: #grade (DO NOT DELETE THIS LINE)

import math

import numpy as np

import pylab as plt

import pyvista as pv

from pyvista import examples

Matplotlib created a temporary config/cache directory at /tmp/matplotlib-

ao5ci0qn because the default path (/tmp/cache/matplotlib) is not a writable

directory; it is highly recommended to set the MPLCONFIGDIR environment variable

to a writable directory, in particular to speed up the import of Matplotlib and

to better support multiprocessing.

[2]: %matplotlib widget

from pyvista import set_plot_theme

set_plot_theme( 'document ')

pv .set_jupyter_backend( 'pythreejs ')

Download a gridded topography surface (DEM) using one of the examples provided by PyVista.

[3]: dem = examples .download_crater_topo()

dem

[3]: UniformGrid (0x7fc34bdac460)

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

[4]: subset = dem .extract_subset((572, 828, 472, 728, 0, 0), (1,1,1))

subset

[4]: UniformGrid (0x7fc34d373e20)

Let’s plot the area we just extracted to see what it looks like.

[5]:  pv .plot(subset)                                                                        

Renderer(camera=PerspectiveCamera(aspect=1 .3333333333333333,␣

,→ children=(DirectionalLight(intensity=0 .25, positi…

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

[6]: terrain = subset .warp_by_scalar() #Warp into a 3D surface mesh (without volume)

terrain

[6]: StructuredGrid (0x7fc34ba76a00)

N Cells:      65536

N Points:     66049

X Bounds:     1 .819e+06, 1 .822e+06

Y Bounds:     5 .647e+06, 5 .651e+06

Z Bounds:     1 .777e+03, 2 .787e+03

Dimensions:   257, 257, 1

N Arrays:     1

We can see that our terrain is now a pyvista .StructuredGrid mesh. Now let’s plot our terrain. [7]:  pv .plot(terrain)                                                                      

Renderer(camera=PerspectiveCamera(aspect=1 .3333333333333333,␣

,→ children=(DirectionalLight(intensity=0 .25, positi…

2.1    Part 2:  Coarsening the Mesh  (and Writing Code)

In this section, you will write code to generate a new coarse mesh from our terrain mesh.  Coarse meshes generally provide less accurate solutions, but are computationally faster.

Your new mesh should be a StructuredGrid, just like the original mesh, but with a lower resolution. This means you will need to redefine the (x, y, z) coordinate points of your mesh. We will explain how to redefine your coordinates a little later on.

First,  let’s  start  with  understanding  how  to  generate  a  new  mesh.    You  can  initialize  a  new StructuredGrid object  directly  from  the  three  point  arrays  that  each  contain  the x,  y,  and  z coordinates of all points in the mesh, respectively. Note: Each array is a 3D array with dimensions M x N x 1 (with the z-axis always being of length 1).

You     will     find     the     following     reference     helpful:         https://docs.pyvista.org/core/point- grids.html#pyvista.StructuredGrid.

Let’s look at the example below for initializing a new StructuredGrid.

[9]: StructuredGrid (0x7fc34af9bf40)

N Cells:      729

N Points:     1000

X Bounds:     -1 .000e+01, 8 .000e+00

Y Bounds:     -1 .000e+01, 8 .000e+00

Z Bounds:     -1 .000e+01, 8 .000e+00

Dimensions:   10, 10, 10

N Arrays:     0

Now, let’s follow the same general steps as in the above example to generate our new coarse mesh from our previously created terrain mesh.

We can coarsen the mesh by merging every 2f quads/cells into one and dropping the center point, where f is your sampling factor aka the factor by which you want to reduce the resolution. In other words, we can produce a reduced version of the mesh by sampling one out of every f points along each axis of the mesh.

Write code to coarsen terrain by a factor of 2.  In other words, we will be converting the mesh from a  257x257-point mesh to a  129x129-point mesh  (or equivalently,  a  256x256-cell mesh to a 128x128-cell mesh).

In the code block below, complete the implementation of coarsen(terrain) by defining three new point arrays, xnew, ynew, and znew and composing them into a new StructuredGrid object named coarse. Your function should return a tuple containing  (xnew, ynew, znew, coarse).

[ ]: ### Tests for coarsenMesh.

### Please DO NOT hard-code the answers as we will also be using hidden test␣ ,→ cases when grading your submission.

print(f"Coarsening from {terrain .dimensions[0]} to {math .ceil(terrain .

,→dimensions[0]/2)}...")

xnew, ynew, znew, coarse = coarsen(terrain)

assert xnew .shape == (129,129,1)

assert ynew .shape == (129,129,1)

np .testing .assert_allclose(xnew[0][0][0],1818580, rtol=1e-7)    np .testing .assert_allclose(xnew[5][120][0],1818730, rtol=1e-7)  np .testing .assert_allclose(ynew[128][120][0],5650680, rtol=1e-7) np .testing .assert_allclose(znew[12][12][0],1880.53, rtol=1e-5)

We can plot the z-values of our new coarsened mesh by adding an additional attribute values to our mesh, which will contain a normalized, column-major flattened representation of the z-axis values of our grid.

See       the       following       reference       for       more       information       on       array       flattening: https://numpy.org/doc/stable/reference/generated/numpy.ndarray.flatten.html

[ ]: #Plot the z-values using the viridis (default) color map

coarse[ 'values '] = pv .plotting .normalize(coarse .z .flatten("F"))

pv .plot(coarse, scalars= 'values ')

2.2    Part 3:  Visualizing Error Values

Now that we have generated our coarse mesh, we can visualize the error in elevation between our coarse mesh and our original mesh. More specifically, we want to compute the error value for each point between the new  (bilinearly interpolated) center point elevation and the original.  We will then visualize the error as a scalar field on the original mesh.

Since we will need to bilinearly interpolate the center point elevation (i.e. the z-value) of each point in our coarse mesh in order to match the dimensions of our original mesh, let’s define a function to do just that.

Define  the  function  bilin() to  bilinearly  interpolate  the  value  at  coordinates  (x,y) within  a rectilinear grid of points.

The  parameters  of your  function  are:  - x = x-coordinate of point whose value we wish to interpolate - y = y-coordinate of point whose value we wish to interpolate - points = a list of four triplets of the form  (xc, yc, val), where val denotes the function value associated with coordinates  (xc, yc)

This function should return a bilinearly interpolated value associated with coordinate  (x,y) w.r.t the rectilinear grid formed by points.

Hints: - You may assume the four triplets within points form a valid rectangle - You may assume x and y fall within the rectangle formed by the points parameter - You should NOT assume the four triplets within points are in any specific order

[ ]: ### Tests for bilin(x, y, points) function.

### Please DO NOT hard-code the answers as we will also be using hidden test␣

,→ cases when grading your submission.

testing_points = [(1,1,3), (3,1,6), (1,3,7), (3,3,9)]

result = bilin(2,2,testing_points)

np .testing .assert_allclose(result,6.25, rtol=1e-2)

result = bilin(2.5,2.5,testing_points)

np .testing .assert_allclose(result,7.6875, rtol=1e-3)

result = bilin(1.1,1.1,testing_points)

np .testing .assert_allclose(result,3.3475, rtol=1e-3)

Now, using your bilin() function, create a new mesh or StructuredGrid object named intmesh,

reconstructed from coarse using bilinear interpolation, with the same dimensions as our original mesh terrain. Your new mesh should contain the interpolated z-values for each point in terrain.

As a starting point, we have defined some of the variables that will be helpful to you for creating your new interpolated mesh.  Specifically, we will be checking the values in errz and intz, as defined below:  -  intz:  a  3D  array with the  same  shape  as terrain .z that will  contain the bilinearly interpolated z-values from the coarsened mesh.Note:  intz is a 3D M x N x  1 array where the last  dimension  contains  the  z-values.   You  should  note  this  when  assigning  elements  to  intz. The  interpolated  mesh  is  a  3D  M  by  N  by  1  array,  so tests will  fail  if you  assign values to  it as if it were just a  2D M by N array,  bypassing the z-axis entirely.   So,  if your code looks like intz[x][y] = bilin( . . .) the fix is to change it to something like intz[x][y] = [ bilin( . . .) ] or intz[x][y][0] = bilin( . . .).

•  errz:  a list of scalar values.   This should contain the absolute error values between each z-value in the original mesh and each interpolated z-value in the new returned mesh

Just like how we added the attribute values to our coarse mesh in order to plot the z-values of the mesh, you should add an additional attribute errors to intmesh in order to plot the absolute error values between the z-values in the original mesh and the interpolated z-values in our new returned mesh.

In  the  code  block  below,  complete  the  implementation  of  reconstructMesh(terrain,coarse),

which should return a tuple containing  (intz, errz, intmesh).

[ ]:  intmesh                                                                                 Now we can visualize the error values that we computed! We recommend adjusting the color map