Contents

2. Meshes with Topography

This notebook introduces the QuagMesh object, which builds the following data structures:

  • hill slope

  • downhill propagation matrices

  • upstream area

in addition to the data structures inherited from QuagMesh. These form the necessary structures to propagate information from higher to lower elevations. Derivatives are computed on the mesh to calculate the height field, smoothing operators are available to reduce short wavelength features and artefacts.

In this notebook we setup a height field and calculate its derivatives on an unstructued mesh. We smooth the derivatives using the radial-basis function (RBF) smoothing kernel.

Note: The API for the structured mesh is identical

from quagmire.tools import meshtools
from quagmire import QuagMesh, QuagMesh
from quagmire import function as fn

import matplotlib.pyplot as plt
import numpy as np
# from scipy.ndimage import imread
# from quagmire import tools as meshtools
# from quagmire import QuagMesh
%matplotlib inline

minX, maxX = -5.0, 5.0
minY, maxY = -5.0, 5.0,
dx, dy = 0.02, 0.02

x, y, simplices = meshtools.elliptical_mesh(minX, maxX, minY, maxY, dx, dy)

DM = meshtools.create_DMPlex_from_points(x, y, bmask=None)

mesh = QuagMesh(DM, downhill_neighbours=1)

print ("Triangulation has {} points".format(mesh.npoints))

Underlying Mesh type: TriMesh
0 - Delaunay triangulation 0.11218824999999555s
0 - Calculate node weights and area 0.020123624999996537s
0 - Find boundaries 0.007983166000002484s
0 - cKDTree 0.006998583000004999s
0 - Construct neighbour cloud arrays 0.2953688329999977s, (0.19133074999999877s + 0.10399454199999525s)
0 - Construct rbf weights 0.028541665999995303s
Triangulation has 62234 points

Height field

We generate a cylindrically symmetry domed surface and add multiple channels incised along the boundary. The height and slope fields reside as attributes on the QuagMesh instance.

radius  = np.sqrt((x**2 + y**2))
theta   = np.arctan2(y,x) + 0.1

height  = np.exp(-0.025*(x**2 + y**2)**2) + 0.25 * (0.2*radius)**4  * np.cos(5.0*theta)**2 ## Less so
height  += 0.5 * (1.0-0.2*radius)
height  += np.random.random(height.size) * 0.01 # random noise

# This fails because the topography variable is locked
mesh.topography.data = height

# This unlocks the variable and rebuilds the necessary downhill data structures
with mesh.deform_topography():
    print("Update topography data array (automatically rebuilds matrices)")
    mesh.topography.data = height
    print("Update topography data array (automatically rebuilds matrices ONCE ONLY)")
    mesh.topography.data = height + 0.01

quagmire.MeshVariable: h(x,y) - is locked
Update topography data array (automatically rebuilds matrices)
Update topography data array (automatically rebuilds matrices ONCE ONLY)
0 - Build downhill matrices 1.1500198329999947s
0 - Build upstream areas 0.5090339589999999s

mesh.topography.data

array([1.50868874, 1.51634837, 1.50855902, ..., 0.70499342, 0.66836758,
       0.61934834])

s = mesh.slope
s
\[\displaystyle \left| \nabla h(x,y) \right|\]

Derivatives and slopes

The slope of the topography is defined through a built in lazy-evaluate function mesh.slope (which was described in the Functions notebook). Other gradients are available through the usual quagmire mathematics functions.

If you want more control of the underlying operations, derivatives can also be evaluated on the mesh using the inbuilt routine in the stripy object. It employs automatically selected tension factors to preserve shape properties of the data and avoid overshoot and undershoot associated with steep gradients. Note: In parallel it is wise to check if this tensioning introduces artefacts near the boundaries.

dfdx, dfdy = mesh.derivative_grad(f, nit=10, tol=1e-8):

where nit and tol control the convergence criteria.

Smoothing

We have included the capacity to build (Gaussian) Radial Basis Function kernels on the mesh that can be used for smoothing operations. Radial-basis function (RBF) smoothing kernel works by setting up a series of gaussian functions based on the distance \(d\) between neighbouring nodes and a scaling factor, \(\Delta\):

\[ W_i = \frac{\exp \left( \frac{d_i}{\Delta} \right)^2}{\sum_{i} \left( \frac{d_i}{\Delta} \right)^2} \]

delta is set to the mean distance between nodes by default, but it may be changed to increase or decrease the smoothness:

rbf1  = mesh.build_rbf_smoother(1.0, 1)
rbf01 = mesh.build_rbf_smoother(0.1, 1)
rbf001 = mesh.build_rbf_smoother(0.01, 1)

print(rbf1.smooth_fn(rainfall, iterations=1).evaluate(0.0,0.0))
print(rbf1.smooth_fn(height, iterations=1).evaluate(0.0,0.0))
print(rbf01.smooth_fn(rainfall, iterations=1).evaluate(0.0,0.0))

rbf005 = mesh.build_rbf_smoother(0.05, 1)
rbf010 = mesh.build_rbf_smoother(0.10, 1)
rbf050 = mesh.build_rbf_smoother(0.50, 1)

rbf_slope005 = rbf005.smooth_fn(mesh.slope).evaluate(mesh)
rbf_slope010 = rbf010.smooth_fn(mesh.slope).evaluate(mesh)
rbf_slope050 = rbf050.smooth_fn(mesh.slope).evaluate(mesh)

NOTE - Building the RBF smoothing machinery is expensive and cannot be reused if the kernel properties are changed. We therefore have a two-stage implementation which builds and caches the smoothing matrices and defines a quagmire function that can be used in the usual way.

import k3d

plot = k3d.plot(grid_visible=False)
indices = mesh.tri.simplices.astype(np.uint32)
points  = np.column_stack([mesh.tri.points, height]).astype(np.float32)

plt_mesh = k3d.mesh(points, indices, wireframe=False, flat_shading=False, side="double",
                    color_map = k3d.colormaps.basic_color_maps.CoolWarm,
                    attribute=mesh.slope.evaluate(mesh),
                    # color_range = [-1.1,2.01]
                   )

plot += plt_mesh
plot.display()