6. Catchment analysis.¶

We start with “Swamp Mountain” from the previous notebooks. This is slightly modified so that there are no lakes / pits right at the boundary.

The catchments are identified by first finding all the outflow points of the mesh (local minima that correspond to the boundary mask) and then using the transpose of the downhill-propagation matrix \(D^T\) to run information (the unique ID of each outflow points) up to the top of the catchment.

The stopping condition is that no further change occurs.

Note in the context of multiple pathways, this operation produces a fuzzy catchment. The first thing we do in this notebook is to specify downhill_neighbours=1

import matplotlib.pyplot as plt
import numpy as np
from quagmire import tools as meshtools
%matplotlib inline

Construct the swamp-mountain mesh¶

This time we take care to avoid lakes etc on the boundaries as this makes the catchment analysis more complicated. Visualise the mesh to make sure that this works.

from quagmire import QuagMesh 
from quagmire import QuagMesh # all routines we need are within this class
from quagmire import QuagMesh

minX, maxX = -5.0, 5.0
minY, maxY = -5.0, 5.0,

spacingX = 0.05
spacingY = 0.05

x, y, simplices = meshtools.elliptical_mesh(minX, maxX, minY, maxY, spacingX, spacingY, 1.)

DM = meshtools.create_DMPlex(x, y, simplices)
DM = meshtools.refine_DM(DM, refinement_levels=2)

mesh = QuagMesh(DM, downhill_neighbours=1)

# Note ... this is what refinement does 
x = mesh.coords[:,0]
y = mesh.coords[:,1]

print( "\nNumber of points in the triangulation: {}".format(mesh.npoints))

Underlying Mesh type: TriMesh
- Delaunay triangulation 0.3174192500000004s
- Calculate node weights and area 0.02463904100000036s
- Find boundaries 0.018807584000001043s
- cKDTree 0.022422916000039095s
- Construct neighbour cloud arrays 0.6958790830000225s, (0.4476076249999892s + 0.24823249999997188s)
- Construct rbf weights 0.07766774999998916s

Number of points in the triangulation: 156925

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  -= height.min()

## Add smoothed random noise to make some "lakes" 

mesh._construct_rbf_weights(delta=mesh.delta*3.0)

randpts1 = np.where(np.random.random(height.shape)>0.995, -1.0, 0.0)
hrand1   = 20.0 * mesh.rbf_smoother(randpts1, iterations=10)
heightn = height + hrand1 * np.exp(-radius**2/15.0) 


with mesh.deform_topography():
    mesh.downhill_neighbours = 2
    mesh.topography.data = heightn


# let's use a rainfall proportional to height (any choice is ok)

rainfall_fn  = mesh.topography ** 2.0
flowrate_fn  = mesh.upstream_integral_fn(rainfall_fn)
str_power_fn = mesh.upstream_integral_fn(rainfall_fn)**2.0 * mesh.slope ** 2.0

0 - Build downhill matrices 3.9534281659999806s
0 - Build upstream areas 1.520386708999979s

Process the surface to fill swamps and lakes¶

## Quagmire also has a swamp filling algorithm

mesh1s = QuagMesh(DM)
with mesh1s.deform_topography():
    mesh1s.topography.data = mesh.topography.data
    
    
mesh1s.low_points_local_patch_fill(its=3, smoothing_steps=2, fraction=0.5)

for i in range(0,5):
    mesh1s.low_points_swamp_fill(ref_height=-0.01)
    
    # In parallel, we can't break if ANY processor has work to do (barrier / sync issue)
    low_points2 = mesh1s.identify_global_low_points()
    
    print("{} : {}".format(i,low_points2[0]))
    if low_points2[0] == 0:
        break
        

Underlying Mesh type: TriMesh
- Delaunay triangulation 0.2958527910000157s
- Calculate node weights and area 0.025907167000013942s
- Find boundaries 0.018733249999968393s
- cKDTree 0.020705749999933687s
- Construct neighbour cloud arrays 0.6505634170000576s, (0.42587704200002463s + 0.22465216600005533s)
- Construct rbf weights 0.05773158399995282s
- Build downhill matrices 2.4506870000000163s
- Build upstream areas 1.5300563330000614s
Low point local patch fill
- Build downhill matrices 2.43383074999997s
- Build downhill matrices 2.4526279579999937s
- Build downhill matrices 2.4396555829999897s
Low point local patch fill  10.608140625000033  seconds
- Build upstream areas 1.5366657079999868s
iterations, time =  1.2963485420000325
Build low point catchments -  1.2967541669999036  seconds
Sort spills -  0.0075850000000627915
Gather spill data -  6.412499999441934e-05
Sort all spills -  5.966699995951785e-05
- Build downhill matrices 2.4458806249999725s
Low point swamp fill  3.838560457999961  seconds
- Build upstream areas 1.7615719579999904s
: 72
iterations, time =  1.3446816670000317
Build low point catchments -  1.3451523749999978  seconds
Sort spills -  0.004237167000042064
Gather spill data -  6.758300003184559e-05
Sort all spills -  4.949999993186793e-05
- Build downhill matrices 2.4374321249999866s
Low point swamp fill  3.833339459000058  seconds
- Build upstream areas 1.972699499999976s
: 16
iterations, time =  1.536034083000004
Build low point catchments -  1.5364633750000394  seconds
Sort spills -  0.001867416999971283
Gather spill data -  6.645799999205337e-05
Sort all spills -  0.00011812500008545612
- Build downhill matrices 2.451272208999967s
Low point swamp fill  4.021867333000046  seconds
- Build upstream areas 2.0030419999999367s
: 7
iterations, time =  1.6317061669999475
Build low point catchments -  1.6321679579999682  seconds
Sort spills -  0.001334499999984473
Gather spill data -  6.400000006578921e-05
Sort all spills -  2.237499995771941e-05
- Build downhill matrices 2.4534335420000843s
Low point swamp fill  4.112114665999911  seconds
- Build upstream areas 1.980670000000032s
: 0

rainfall_fn_1s  = mesh1s.topography ** 2.0
flowrate_fn_1s  = mesh1s.upstream_integral_fn(rainfall_fn_1s)
str_power_fn_1s = mesh1s.upstream_integral_fn(rainfall_fn_1s)**2.0 * mesh1s.slope ** 2.0

Locating and viewing the outflow points¶

quagmire provides a mechanism to find the outflow points of a domain and return the node values. Note: in parallel these are the local node numbers and are not a unique ID. To do this, we can obtain the global ID from PETSc but it does also help to have the indices all be small numbers so we can map colours effectively.

outflows = mesh1s.identify_outflow_points()
print("Mesh has {} outflows".format(outflows.shape[0]))

Mesh has 21 outflows

import quagmire
print(quagmire.mesh.check_object_is_a_q_mesh_and_raise(mesh1s))

True

import k3d


low_points = mesh.identify_low_points()
lowsxyz = np.column_stack([mesh.tri.points[low_points], mesh.topography.data[low_points]])

xyz  = np.column_stack([mesh1s.tri.points, mesh1s.topography.data])
xyz2 = np.column_stack([mesh.tri.points, mesh.topography.data])

indices = mesh.tri.simplices.astype(np.uint32)

plot = k3d.plot(camera_auto_fit=False, grid_visible=False, lighting=0.75)

surface = k3d.mesh(xyz2, indices, wireframe=False, flat_shading=True, side="double",
                    color_map = k3d.colormaps.matplotlib_color_maps.Blues,
                    attribute=flowrate_fn_1s.evaluate(mesh),
                    color_range = [0.0,1.0]
                   )

lakes_surface = k3d.mesh(xyz-(0.0,0.0,0.005), indices, wireframe=False, flat_shading=True, side="double",
                    color=0x0010C0, opacity=0.5
                   )

lows    = k3d.points(lowsxyz, point_size=0.05, color=0xFF0000)


plot += surface + lows + lakes_surface

plot.display()

## Stream power / slope where the lakes / swamps are located:

Set neighbours to 1 and compute “uphill” connectivity¶

In serial, i.e. for this demonstration, we number the outflow points by their order in the local mesh numbering. We can then use the mesh.uphill_propagation routine to propagate this information from the outflow to the top of the catchment.

This routine returns the mesh data (for this processor) of a globally synchronised map of the information propagated from the selected points.

The routine is used in the computation of flood fills for the swamp algorithm and should be polished up a bit (sorry).

## Unique catchments requires the downhill matrix with downhill_neighbours=1

mesh1s.near_neighbours=1
# mesh1s.update_height(mesh1s.heightVariable.data)

## Need a unique ID that works in parallel ... global node number would work but 
## not that easy to map to colours in lavavu 

from petsc4py import PETSc
outflows
outflowID = mesh1s.lgmap_row.apply(outflows.astype(PETSc.IntType))

# But on 1 proc, this is easier / better:

outflowID = np.array(range(0, outflows.shape[0]))
ctmt = mesh1s.uphill_propagation(outflows,  outflowID, its=99999, fill=-999999).astype(np.int)

282  iterations, time =  1.5911248339999702

<ipython-input-20-a1316dff9ef9>:11: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  ctmt = mesh1s.uphill_propagation(outflows,  outflowID, its=99999, fill=-999999).astype(np.int)

Visualise the catchment information¶

import k3d

xyz  = np.column_stack([mesh1s.tri.points, mesh1s.topography.data])
xyz2 = np.column_stack([mesh.tri.points,   mesh.topography.data])

indices = mesh.tri.simplices.astype(np.uint32)

plot = k3d.plot(camera_auto_fit=False, grid_visible=False, lighting=0.75)

surface = k3d.mesh(xyz2, indices, wireframe=False, flat_shading=True, side="double",
                    color_map = k3d.colormaps.matplotlib_color_maps.Blues,
                    attribute=flowrate_fn_1s.evaluate(mesh),
                    color_range = [0.0,1.0]
                   )

lakes_surface = k3d.mesh(xyz-(0.0,0.0,0.01), indices, wireframe=False, flat_shading=True, side="double",
                    color=0x0010C0, opacity=0.5)


catchments = k3d.mesh(xyz2+(0.0,0.0,0.01), indices, wireframe=False, flat_shading=True, side="double",
                    color_map = k3d.colormaps.matplotlib_color_maps.Greens,
                    attribute=ctmt, opacity=0.5
                    # color_range = [0.0,1.0]
                   )




plot += surface + lows + lakes_surface + catchments

plot.display()

/Users/lmoresi/mambaforge/envs/underworld/lib/python3.9/site-packages/traittypes/traittypes.py:97: UserWarning: Given trait value dtype "float64" does not match required type "float32". A coerced copy has been created.
  warnings.warn(

Quagmire.tex