5. Preprocessing a surface¶

Pit filling and swamp filling …

Take the previous mesh with random noise

Quagmire allows the user to specify the number of downhill pathways to model landscape evolution. This is set using:

mesh.downhill_neighbours = 1
mesh.update_height(height)

where an integer specifies the number of downhill neighbour nodes (receipients) that will receive a packet of information from a donor node. The QuagMesh object can also be initialised with:

mesh = QuagMesh(DM, downhill_neighbours=1)

to specify the number of downhill neighbours (default is 2).

In this notebook we use a landscape function with many outflow points to examine the effect of varying the number of recipient nodes on catchment area, stream lengths, and outflow fluxes.

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

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, refinement_levels=2)

mesh = QuagMesh(DM)

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.46267416599999933s
- Calculate node weights and area 0.043066750000001264s
- Find boundaries 0.02290166599999921s
- cKDTree 0.02247120799999891s
- Construct neighbour cloud arrays 0.6619863330000015s, (0.44069049999999876s + 0.2212610419999983s)
- Construct rbf weights 0.05931295799999958s

Number of points in the triangulation: 156925

# topography

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 some pits 

hrand0 = np.where(np.random.random(height.shape)>0.995, -0.1, 0.0)


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

rbf_smoother = mesh.build_rbf_smoother(0.05, iterations=3)
h0 = mesh.add_variable(name="h0")
h0.data = np.where(np.random.random(height.shape)>0.995, -0.33, 0.0)

hrand1 = 25.0 * rbf_smoother.smooth_fn(h0, iterations=25).evaluate(mesh)


# randpts1 = np.where(np.random.random(height.shape)>0.995, -1.0, 0.0)
# hrand1   = 10.0 * rbf_smoother.smooth(randpts1, iterations=10)

heightn = np.maximum(0.0, height + hrand0 + hrand1)

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

low_points = mesh.identify_low_points()
low_point_coords = mesh.coords[low_points] 
print(low_points.shape)

0 - Build downhill matrices 3.843434333000001s
0 - Build upstream areas 1.5636943749999972s
(886,)

h0.data = mesh.topography.data

mesh

import k3d

lowsxyz = np.column_stack([mesh.tri.points[low_points], heightn[low_points]])
xyz  = np.column_stack([mesh.tri.points, height])
xyz2 = np.column_stack([mesh.tri.points, heightn])

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

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

surface = k3d.mesh(xyz2, indices, wireframe=False, flat_shading=True, side="double",
                    color_map = k3d.colormaps.matplotlib_color_maps.Cividis,
                    attribute=mesh.topography.data,
                    color_range = [0.0,1.5]
                   )

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

plot += surface + lows

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(

rainfall_fn = mesh.topography**2
flowrate_fn = mesh.upstream_integral_fn(rainfall_fn)
stream_power_fn = flowrate_fn ** 2.0 * mesh.slope ** 2.0 * mesh.mask

cumulative_rain_n1 = mesh.upstream_integral_fn(rainfall_fn).evaluate(mesh)

import k3d

lowsxyz = np.column_stack([mesh.tri.points[low_points], heightn[low_points]])
xyz  = np.column_stack([mesh.tri.points, height])
xyz2 = np.column_stack([mesh.tri.points, heightn])
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=stream_power_fn.evaluate(mesh),
                    color_range = [0.0,0.01]
                   )

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

plot += surface + lows

plot.display()

## Pit filling algorithm in quagmire

mesh1p = QuagMesh(DM)

rainfall_fn_1p = mesh1p.topography**2
cumulative_rain_n1p = mesh1p.upstream_integral_fn(rainfall_fn_1p)
stream_power_n1p    = cumulative_rain_n1p ** 2.0 * mesh1p.slope ** 2.0 * mesh1p.mask

with mesh1p.deform_topography():
    mesh1p.topography.data = mesh.topography.data

Underlying Mesh type: TriMesh
- Delaunay triangulation 0.3040237499999421s
- Calculate node weights and area 0.0241118330000063s
- Find boundaries 0.01836920799996733s
- cKDTree 0.019008834000032948s
- Construct neighbour cloud arrays 0.6636431250001351s, (0.43505416700008936s + 0.2285542079998777s)
- Construct rbf weights 0.06518099999993865s
- Build downhill matrices 2.430833333999999s
- Build upstream areas 1.5887439580001228s

mesh1p.low_points_local_patch_fill(its=2, smoothing_steps=1, fraction=0.25)

Low point local patch fill
0 - Build downhill matrices 2.4737774160000754s
0 - Build downhill matrices 2.4994586660000095s
Low point local patch fill  7.116479208999863  seconds
0 - Build upstream areas 1.5854269160001877s

mesh1p.identify_low_points().shape

(193,)

import k3d

low_points2 = mesh1p.low_points
height2 = mesh1p.topography.data

lowsxyz = np.column_stack([mesh1p.tri.points[low_points], heightn[low_points]])
lows2xyz = np.column_stack([mesh1p.tri.points[low_points2], height2[low_points2]])

xyz  = np.column_stack([mesh.tri.points, mesh1p.topography.data])
xyz2 = np.column_stack([mesh.tri.points, heightn])
indices = mesh.tri.simplices.astype(np.uint32)

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


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

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


plot += surface + lows + lows2 

plot.display()

## Quagmire also has a swamp filling algorithm
## NOTE this is much more efficient if it follows the pit filling

mesh1s = QuagMesh(DM)
rainfall_fn_1s = mesh1s.topography**2

with mesh1s.deform_topography():
    mesh1s.topography.data = mesh1p.topography.data
    
print("There are still {} low points".format(mesh1s.identify_low_points().shape))

Underlying Mesh type: TriMesh
- Delaunay triangulation 0.3330004590000044s
- Calculate node weights and area 0.0643369580000126s
- Find boundaries 0.018808749999834617s
- cKDTree 0.0196499579999454s
- Construct neighbour cloud arrays 0.6992000420000295s, (0.44694904200014207s + 0.25221525000006295s)
- Construct rbf weights 0.058850083999914204s
- Build downhill matrices 2.4509212079999543s
- Build upstream areas 1.571919041000001s
(193,)

for i in range(0,50):
    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

iterations, time =  1.2702137920000496
Build low point catchments -  1.2706606249998913  seconds
Sort spills -  0.010811040999897159
Gather spill data -  0.0019107919999896694
Sort all spills -  0.00026254099998368474
- Build downhill matrices 2.4735215000000608s
Low point swamp fill  3.997965583999985  seconds
- Build upstream areas 1.8586511250000513s
: 96
iterations, time =  1.3629859999998644
Build low point catchments -  1.3634352080000554  seconds
Sort spills -  0.005808292000210713
Gather spill data -  7.35829999030102e-05
Sort all spills -  9.187499995277904e-05
- Build downhill matrices 2.4238212090001525s
Low point swamp fill  3.935164916999838  seconds
- Build upstream areas 1.8560250830000768s
: 70
iterations, time =  1.4846040420000008
Build low point catchments -  1.485055124999917  seconds
Sort spills -  0.004419750000124623
Gather spill data -  6.962499992368976e-05
Sort all spills -  8.179099995686556e-05
- Build downhill matrices 2.443070875000103s
Low point swamp fill  4.039379707999842  seconds
- Build upstream areas 2.0338651250001476s
: 62
iterations, time =  1.5380792920000204
Build low point catchments -  1.5385104590000083  seconds
Sort spills -  0.003909999999905267
Gather spill data -  0.00012487499998314888
Sort all spills -  7.504199993491056e-05
- Build downhill matrices 2.5028409589999683s
Low point swamp fill  4.142744708000009  seconds
- Build upstream areas 2.005432916000018s
: 70
iterations, time =  1.7582492910000838
Build low point catchments -  1.758841208999911  seconds
Sort spills -  0.006341291000126148
Gather spill data -  0.00023295899995900982
Sort all spills -  0.00010766600007627858
- Build downhill matrices 3.0342551250000724s
Low point swamp fill  4.935335333000012  seconds
- Build upstream areas 2.779530582999996s
: 64
iterations, time =  2.100522917000035
Build low point catchments -  2.10303458400017  seconds
Sort spills -  0.003463125000052969
Gather spill data -  0.0029077920000872837
Sort all spills -  9.208300002683245e-05
- Build downhill matrices 3.149582374999909s
Low point swamp fill  5.405693749999955  seconds
- Build upstream areas 2.578506583000035s
: 46
iterations, time =  1.7416526659999363
Build low point catchments -  1.742177791000131  seconds
Sort spills -  0.002930999999989581
Gather spill data -  7.379199996648822e-05
Sort all spills -  5.487500015988189e-05
- Build downhill matrices 2.9200388749998183s
Low point swamp fill  4.733386958999972  seconds
- Build upstream areas 2.319822375000058s
: 38
iterations, time =  1.8350086250000004
Build low point catchments -  1.8355307499998617  seconds
Sort spills -  0.0033648750002157612
Gather spill data -  7.483300009880622e-05
Sort all spills -  5.600000008598727e-05
- Build downhill matrices 2.569668374999992s
Low point swamp fill  4.467106000000058  seconds
- Build upstream areas 2.469780124999943s
: 30
iterations, time =  1.621527874999856
Build low point catchments -  1.6220797500000117  seconds
Sort spills -  0.0024131250002028537
Gather spill data -  8.974999991551158e-05
Sort all spills -  4.429099999470054e-05
- Build downhill matrices 2.603554832999862s
Low point swamp fill  4.277322541999865  seconds
- Build upstream areas 2.067419541999925s
: 28
iterations, time =  1.6253727079999862
Build low point catchments -  1.6258912910000163  seconds
Sort spills -  0.002284459000065908
Gather spill data -  9.158300008493825e-05
Sort all spills -  3.85000000733271e-05
- Build downhill matrices 2.548146124999903s
Low point swamp fill  4.217911333000075  seconds
- Build upstream areas 2.6503998749999482s
: 15
iterations, time =  1.6346717500000523
Build low point catchments -  1.6351971250001043  seconds
Sort spills -  0.0018603329999677953
Gather spill data -  6.962499992368976e-05
Sort all spills -  0.000473916999908397
- Build downhill matrices 2.594370875000095s
Low point swamp fill  4.272328667000011  seconds
- Build upstream areas 2.288963542000147s
: 12
iterations, time =  1.8006609159999698
Build low point catchments -  1.8014930829999685  seconds
Sort spills -  0.0017905830000017886
Gather spill data -  7.483399986085715e-05
Sort all spills -  0.00014474999989033677
- Build downhill matrices 3.254716791999954s
Low point swamp fill  5.0976602920000005  seconds
- Build upstream areas 2.1872506670001712s
: 10
iterations, time =  1.6066588749999937
Build low point catchments -  1.6071492080000098  seconds
Sort spills -  0.0014944170000035228
Gather spill data -  9.075000002667366e-05
Sort all spills -  3.300000003036985e-05
- Build downhill matrices 2.5504909999999654s
Low point swamp fill  4.190495874999897  seconds
- Build upstream areas 2.1012557080000533s
: 5
iterations, time =  1.5801827910001975
Build low point catchments -  1.580588834000082  seconds
Sort spills -  0.0013571249999131396
Gather spill data -  6.670800007668731e-05
Sort all spills -  2.212499998677231e-05
- Build downhill matrices 2.490343208000013s
Low point swamp fill  4.101697958000159  seconds
- Build upstream areas 2.1326142080001773s
: 7
iterations, time =  1.7845102500000394
Build low point catchments -  1.7850002919999497  seconds
Sort spills -  0.0023967919998995058
Gather spill data -  7.13749998340063e-05
Sort all spills -  2.375000008214556e-05
- Build downhill matrices 2.4884266659998957s
Low point swamp fill  4.30308137499992  seconds
- Build upstream areas 2.1146000000001095s
: 5
iterations, time =  1.5713676249999935
Build low point catchments -  1.5717985830001453  seconds
Sort spills -  0.0011000420001892053
Gather spill data -  0.00017608300004212651
Sort all spills -  3.087499999310239e-05
- Build downhill matrices 2.5157320420000815s
Low point swamp fill  4.114825457999814  seconds
- Build upstream areas 2.0815340420001576s
: 3
iterations, time =  1.5753376669999852
Build low point catchments -  1.5758753329998854  seconds
Sort spills -  0.0010460419998707948
Gather spill data -  6.816700010858767e-05
Sort all spills -  1.9374999965293682e-05
- Build downhill matrices 2.4807477500000914s
Low point swamp fill  4.085463207999965  seconds
- Build upstream areas 2.1508202499999243s
: 1
iterations, time =  1.5369610840000405
Build low point catchments -  1.5374124999998457  seconds
Sort spills -  0.0006748330001755676
Gather spill data -  7.76249999034917e-05
Sort all spills -  1.9416999975874205e-05
- Build downhill matrices 2.540103290999923s
Low point swamp fill  4.099138374999939  seconds
- Build upstream areas 2.0256681250000383s
: 0

cumulative_rain_n1s = mesh1s.upstream_integral_fn(rainfall_fn_1s)
stream_power_n1s    = cumulative_rain_n1s ** 2.0 * mesh1s.slope ** 2.0 * mesh1s.mask

import k3d

low_points2 = mesh1p.low_points
height2 = mesh1p.topography.data

lowsxyz = np.column_stack([mesh1p.tri.points[low_points], heightn[low_points]])
lows2xyz = np.column_stack([mesh1p.tri.points[low_points2], height2[low_points2]])

xyz  = np.column_stack([mesh1s.tri.points, mesh1s.topography.data])
xyz2 = np.column_stack([mesh1p.tri.points, mesh1p.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=cumulative_rain_n1s.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)
lows2    = k3d.points(lows2xyz, point_size=0.05, color=0x0000FF)


plot += surface + lows + lows2 + lakes_surface

plot.display()

Quagmire.tex

5. Preprocessing a surface¶