Mixed-element Meshing - Toy Problem#
This workflow provides a testbed for testing features of mixed-element mesh wotkflow. The objects of river tree and huc boundary are a created to mimic what we get from NHDPlus. This example tests the case of two rivers, which is rare but still needs to be considered to make the workflow general
import watershed_workflow.ui
watershed_workflow.ui.setup_logging(1)
## Importing all the package needed to run the workflow
import os
import numpy as np
from matplotlib import pyplot as plt
import shapely
import logging
import copy
import geopandas as gpd
import sys
import watershed_workflow
import watershed_workflow.utils
import watershed_workflow.mesh
import watershed_workflow.resampling
import watershed_workflow.condition
crs = None
## NHDPlus input type objects
# creating watershed polygon
ws1 = shapely.geometry.Polygon(100*np.array([[0,0],[1,0],[2,0],[3,0],[4,0],[4,1],[4,2],[4,3],[4,4],[3,4.5],[2,5],[1,4.5], [0,4], [0,3], [0,2] ,[0,1]], 'd'))
my_hucs = gpd.GeoDataFrame()
my_hucs['geometry'] = [ws1,]
# creating reaches miminc NHDPlus flowlines
reach1=shapely.geometry.LineString([(200,200),(200,0)])
reach2=shapely.geometry.LineString([ (50,300), (100,300), (100,200),(200,200)])
reach3=shapely.geometry.LineString([ (350,400),(350,300),(300,300),(300,200),(200,200)])
reach4=shapely.geometry.LineString([(100,400), (200,300)])
reaches = gpd.GeoDataFrame()
reaches.geometry = [reach1, reach2, reach3, reach4]
def createCopy(watershed, rivers):
"""To compare before/after, we often want to create copies. Note in real workflows most things are done in-place without copies."""
return watershed.deepcopy(), [r.deepcopy() for r in rivers]
# creating watershed_workflow.split_hucs.SplitHUCs object for demo watershed
watershed = watershed_workflow.split_hucs.SplitHUCs(my_hucs)
2026-02-09 11:42:21,212 - root - INFO: Removing holes on 1 polygons
2026-02-09 11:42:21,212 - root - INFO: -- removed interior
2026-02-09 11:42:21,212 - root - INFO: -- union
2026-02-09 11:42:21,212 - root - INFO: Parsing 1 components for holes
2026-02-09 11:42:21,213 - root - INFO: -- complete
rivers = watershed_workflow.river_tree.createRivers(reaches, method='geometry')
# copy the originals
watershed_orig, rivers_orig = watershed, rivers
# simplify the geometry for a "nicer" discrete object, using a uniform sampling strategy
watershed1, rivers1 = createCopy(watershed_orig, rivers_orig)
watershed_workflow.simplify(watershed1, rivers1, 10, keep_points=True)
2026-02-09 11:42:21,222 - root - INFO:
2026-02-09 11:42:21,222 - root - INFO: Simplifying
2026-02-09 11:42:21,222 - root - INFO: ------------------------------
2026-02-09 11:42:21,222 - root - INFO: None
2026-02-09 11:42:21,222 - root - INFO: Presimplify to remove colinear, coincident points.
2026-02-09 11:42:21,222 - root - INFO: None
2026-02-09 11:42:21,223 - root - INFO: Pruning leaf reaches < 10
2026-02-09 11:42:21,223 - root - INFO: None
2026-02-09 11:42:21,223 - root - INFO: Merging internal reaches < 10
2026-02-09 11:42:21,223 - root - INFO: None
2026-02-09 11:42:21,224 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,224 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,224 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,224 - root - INFO:
2026-02-09 11:42:21,224 - root - INFO: HUC : min seg length: 223.6067977500 min geom length: 1647.2135955000
2026-02-09 11:42:21,225 - root - INFO: HUC : med seg length: 400.0000000000 med geom length: 1647.2135955000
2026-02-09 11:42:21,225 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,225 - root - INFO:
2026-02-09 11:42:21,225 - root - INFO: Snapping discrete points to make rivers and HUCs discretely consistent.
2026-02-09 11:42:21,225 - root - INFO: -- snapping HUC triple junctions to reaches
2026-02-09 11:42:21,226 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,226 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,226 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,226 - root - INFO:
2026-02-09 11:42:21,226 - root - INFO: HUC : min seg length: 223.6067977500 min geom length: 1647.2135955000
2026-02-09 11:42:21,226 - root - INFO: HUC : med seg length: 400.0000000000 med geom length: 1647.2135955000
2026-02-09 11:42:21,227 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,227 - root - INFO:
2026-02-09 11:42:21,227 - root - INFO: None
2026-02-09 11:42:21,227 - root - INFO: -- snapping reach endpoints to HUC boundaries
2026-02-09 11:42:21,229 - root - INFO: new_lss: min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,229 - root - INFO: new_lss: med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,229 - root - INFO: new_lss: max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,229 - root - INFO:
2026-02-09 11:42:21,230 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,230 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,231 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,231 - root - INFO:
2026-02-09 11:42:21,231 - root - INFO: HUC : min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,231 - root - INFO: HUC : med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,231 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,231 - root - INFO:
2026-02-09 11:42:21,232 - root - INFO: None
2026-02-09 11:42:21,232 - root - INFO: -- cutting reaches at HUC boundaries
2026-02-09 11:42:21,232 - root - INFO: intersection found
2026-02-09 11:42:21,232 - root - INFO: - cutting reach at external boundary of HUCs:
2026-02-09 11:42:21,233 - root - INFO: split HUC boundary ls into 1 pieces
2026-02-09 11:42:21,233 - root - INFO: split reach ls into 1 pieces
2026-02-09 11:42:21,235 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,235 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,235 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,235 - root - INFO:
2026-02-09 11:42:21,235 - root - INFO: HUC : min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,235 - root - INFO: HUC : med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,236 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,236 - root - INFO:
2026-02-09 11:42:21,236 - root - INFO: None
2026-02-09 11:42:21,236 - root - INFO:
2026-02-09 11:42:21,236 - root - INFO: Simplification Diagnostics
2026-02-09 11:42:21,236 - root - INFO: ------------------------------
2026-02-09 11:42:21,237 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,237 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,237 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,237 - root - INFO:
2026-02-09 11:42:21,237 - root - INFO: HUC : min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,238 - root - INFO: HUC : med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,238 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,238 - root - INFO:
2026-02-09 11:42:21,238 - root - INFO: None
2026-02-09 11:42:21,238 - root - INFO:
2026-02-09 11:42:21,238 - root - INFO: Resampling HUC and river
2026-02-09 11:42:21,239 - root - INFO: ------------------------------
2026-02-09 11:42:21,239 - root - INFO: -- resampling HUCs based on uniform target 10
2026-02-09 11:42:21,241 - root - INFO: None
2026-02-09 11:42:21,241 - root - INFO: -- resampling reaches based on uniform target 10
2026-02-09 11:42:21,243 - root - INFO: None
2026-02-09 11:42:21,243 - root - INFO:
2026-02-09 11:42:21,243 - root - INFO: Resampling Diagnostics
2026-02-09 11:42:21,243 - root - INFO: ------------------------------
2026-02-09 11:42:21,244 - root - INFO: reach: min seg length: 9.4280904158 min geom length: 141.4213562373
2026-02-09 11:42:21,244 - root - INFO: reach: med seg length: 10.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,244 - root - INFO: reach: max seg length: 10.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,244 - root - INFO:
2026-02-09 11:42:21,245 - root - INFO: HUC : min seg length: 9.7220346848 min geom length: 1647.2135955000
2026-02-09 11:42:21,245 - root - INFO: HUC : med seg length: 10.0000000000 med geom length: 1647.2135955000
2026-02-09 11:42:21,245 - root - INFO: HUC : max seg length: 10.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,245 - root - INFO:
2026-02-09 11:42:21,245 - root - INFO: None
2026-02-09 11:42:21,245 - root - INFO:
2026-02-09 11:42:21,246 - root - INFO: Clean up sharp angles, both internally and at junctions.
2026-02-09 11:42:21,246 - root - INFO: ------------------------------
2026-02-09 11:42:21,248 - root - INFO: ... cleaned up 0 sharp angles on HUCs.
2026-02-09 11:42:21,249 - root - INFO: ... cleaned up 0 sharp angles at outlets on river 0.
2026-02-09 11:42:21,252 - root - INFO: ... cleaned up 0 internal sharp angles on river 0.
2026-02-09 11:42:21,252 - root - INFO: ... cleaned up 0 sharp angles at outlets on river 1.
2026-02-09 11:42:21,252 - root - INFO: ... cleaned up 0 internal sharp angles on river 1.
2026-02-09 11:42:21,253 - root - INFO: Cleaned up 0 sharp angles.
2026-02-09 11:42:21,253 - root - INFO: reach: min seg length: 9.4280904158 min geom length: 141.4213562373
2026-02-09 11:42:21,253 - root - INFO: reach: med seg length: 10.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,253 - root - INFO: reach: max seg length: 10.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,253 - root - INFO:
2026-02-09 11:42:21,254 - root - INFO: HUC : min seg length: 9.7220346848 min geom length: 1647.2135955000
2026-02-09 11:42:21,254 - root - INFO: HUC : med seg length: 10.0000000000 med geom length: 1647.2135955000
2026-02-09 11:42:21,254 - root - INFO: HUC : max seg length: 10.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,254 - root - INFO:
2026-02-09 11:42:21,254 - root - INFO: None
# simplify the geometry using a distance-function based resampling
watershed2, rivers2 = createCopy(watershed_orig, rivers_orig)
watershed_workflow.simplify(watershed2, rivers2, 10, 50, 50, 150, keep_points=True)
2026-02-09 11:42:21,258 - root - INFO:
2026-02-09 11:42:21,258 - root - INFO: Simplifying
2026-02-09 11:42:21,258 - root - INFO: ------------------------------
2026-02-09 11:42:21,258 - root - INFO: None
2026-02-09 11:42:21,259 - root - INFO: Presimplify to remove colinear, coincident points.
2026-02-09 11:42:21,259 - root - INFO: None
2026-02-09 11:42:21,259 - root - INFO: Pruning leaf reaches < 10
2026-02-09 11:42:21,259 - root - INFO: None
2026-02-09 11:42:21,260 - root - INFO: Merging internal reaches < 10
2026-02-09 11:42:21,260 - root - INFO: None
2026-02-09 11:42:21,261 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,261 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,261 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,261 - root - INFO:
2026-02-09 11:42:21,261 - root - INFO: HUC : min seg length: 223.6067977500 min geom length: 1647.2135955000
2026-02-09 11:42:21,261 - root - INFO: HUC : med seg length: 400.0000000000 med geom length: 1647.2135955000
2026-02-09 11:42:21,262 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,262 - root - INFO:
2026-02-09 11:42:21,262 - root - INFO: Snapping discrete points to make rivers and HUCs discretely consistent.
2026-02-09 11:42:21,262 - root - INFO: -- snapping HUC triple junctions to reaches
2026-02-09 11:42:21,263 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,263 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,263 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,263 - root - INFO:
2026-02-09 11:42:21,263 - root - INFO: HUC : min seg length: 223.6067977500 min geom length: 1647.2135955000
2026-02-09 11:42:21,263 - root - INFO: HUC : med seg length: 400.0000000000 med geom length: 1647.2135955000
2026-02-09 11:42:21,263 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,264 - root - INFO:
2026-02-09 11:42:21,264 - root - INFO: None
2026-02-09 11:42:21,264 - root - INFO: -- snapping reach endpoints to HUC boundaries
2026-02-09 11:42:21,266 - root - INFO: new_lss: min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,266 - root - INFO: new_lss: med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,266 - root - INFO: new_lss: max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,267 - root - INFO:
2026-02-09 11:42:21,267 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,268 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,268 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,268 - root - INFO:
2026-02-09 11:42:21,268 - root - INFO: HUC : min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,268 - root - INFO: HUC : med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,268 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,268 - root - INFO:
2026-02-09 11:42:21,269 - root - INFO: None
2026-02-09 11:42:21,269 - root - INFO: -- cutting reaches at HUC boundaries
2026-02-09 11:42:21,269 - root - INFO: intersection found
2026-02-09 11:42:21,270 - root - INFO: - cutting reach at external boundary of HUCs:
2026-02-09 11:42:21,270 - root - INFO: split HUC boundary ls into 1 pieces
2026-02-09 11:42:21,270 - root - INFO: split reach ls into 1 pieces
2026-02-09 11:42:21,272 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,272 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,272 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,272 - root - INFO:
2026-02-09 11:42:21,272 - root - INFO: HUC : min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,272 - root - INFO: HUC : med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,272 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,273 - root - INFO:
2026-02-09 11:42:21,273 - root - INFO: None
2026-02-09 11:42:21,273 - root - INFO:
2026-02-09 11:42:21,273 - root - INFO: Simplification Diagnostics
2026-02-09 11:42:21,273 - root - INFO: ------------------------------
2026-02-09 11:42:21,274 - root - INFO: reach: min seg length: 50.0000000000 min geom length: 141.4213562373
2026-02-09 11:42:21,274 - root - INFO: reach: med seg length: 100.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,274 - root - INFO: reach: max seg length: 200.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,274 - root - INFO:
2026-02-09 11:42:21,274 - root - INFO: HUC : min seg length: 200.0000000000 min geom length: 1647.2135955000
2026-02-09 11:42:21,274 - root - INFO: HUC : med seg length: 223.6067977500 med geom length: 1647.2135955000
2026-02-09 11:42:21,274 - root - INFO: HUC : max seg length: 400.0000000000 max geom length: 1647.2135955000
2026-02-09 11:42:21,275 - root - INFO:
2026-02-09 11:42:21,275 - root - INFO: None
2026-02-09 11:42:21,275 - root - INFO:
2026-02-09 11:42:21,275 - root - INFO: Resampling HUC and river
2026-02-09 11:42:21,275 - root - INFO: ------------------------------
2026-02-09 11:42:21,275 - root - INFO: -- resampling HUCs based on distance function (50, 10, 150, 50)
2026-02-09 11:42:21,355 - root - INFO: None
2026-02-09 11:42:21,356 - root - INFO: -- resampling reaches based on uniform target 10
2026-02-09 11:42:21,357 - root - INFO: None
2026-02-09 11:42:21,358 - root - INFO:
2026-02-09 11:42:21,358 - root - INFO: Resampling Diagnostics
2026-02-09 11:42:21,358 - root - INFO: ------------------------------
2026-02-09 11:42:21,358 - root - INFO: reach: min seg length: 9.4280904158 min geom length: 141.4213562373
2026-02-09 11:42:21,359 - root - INFO: reach: med seg length: 10.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,359 - root - INFO: reach: max seg length: 10.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,359 - root - INFO:
2026-02-09 11:42:21,359 - root - INFO: HUC : min seg length: 8.7899769497 min geom length: 1647.2135955000
2026-02-09 11:42:21,359 - root - INFO: HUC : med seg length: 11.6648206334 med geom length: 1647.2135955000
2026-02-09 11:42:21,360 - root - INFO: HUC : max seg length: 49.9229443249 max geom length: 1647.2135955000
2026-02-09 11:42:21,360 - root - INFO:
2026-02-09 11:42:21,360 - root - INFO: None
2026-02-09 11:42:21,360 - root - INFO:
2026-02-09 11:42:21,360 - root - INFO: Clean up sharp angles, both internally and at junctions.
2026-02-09 11:42:21,360 - root - INFO: ------------------------------
2026-02-09 11:42:21,361 - root - INFO: ... cleaned up 0 sharp angles on HUCs.
2026-02-09 11:42:21,363 - root - INFO: ... cleaned up 0 sharp angles at outlets on river 0.
2026-02-09 11:42:21,365 - root - INFO: ... cleaned up 0 internal sharp angles on river 0.
2026-02-09 11:42:21,366 - root - INFO: ... cleaned up 0 sharp angles at outlets on river 1.
2026-02-09 11:42:21,366 - root - INFO: ... cleaned up 0 internal sharp angles on river 1.
2026-02-09 11:42:21,366 - root - INFO: Cleaned up 0 sharp angles.
2026-02-09 11:42:21,367 - root - INFO: reach: min seg length: 9.4280904158 min geom length: 141.4213562373
2026-02-09 11:42:21,367 - root - INFO: reach: med seg length: 10.0000000000 med geom length: 225.0000000000
2026-02-09 11:42:21,367 - root - INFO: reach: max seg length: 10.0000000000 max geom length: 350.0000000000
2026-02-09 11:42:21,367 - root - INFO:
2026-02-09 11:42:21,367 - root - INFO: HUC : min seg length: 8.7899769497 min geom length: 1647.2135955000
2026-02-09 11:42:21,367 - root - INFO: HUC : med seg length: 11.6648206334 med geom length: 1647.2135955000
2026-02-09 11:42:21,368 - root - INFO: HUC : max seg length: 49.9229443249 max geom length: 1647.2135955000
2026-02-09 11:42:21,368 - root - INFO:
2026-02-09 11:42:21,368 - root - INFO: None
# plot the results
def plot(ws, rivs, ax=None):
if ax is None:
fig, ax = plt.subplots(1, 1, figsize=figsize)
ws.plot(color='k', marker='+', markersize=20, ax=ax)
for river in rivs:
river.plot(marker='x', markersize=20, ax=ax)
fig, axs = plt.subplots(1,3)
plot(watershed_orig, rivers_orig, axs[0])
axs[0].set_title('original')
plot(watershed1, rivers1, axs[1])
axs[1].set_title('uniform resampling')
plot(watershed2, rivers2, axs[2])
axs[2].set_title('dist-based resampling')
Text(0.5, 1.0, 'dist-based resampling')
Surface Meshing#
## Triangulation
refine_d0 = 12
refine_d1 = 18
refine_A0 = 82
refine_A1 = 500
d0 = refine_d0; d1 = refine_d1
A0 = refine_A0; A1 = refine_A1
# Refine triangles if they get too acute
min_angle = 32 # degrees
# function that returns width for each reach. In this case, we use a constant width of 5.0 for all reaches
def widths(reach):
return 5.0
# created two versions of the discrete watershed and rivers -- let's use the second one
watershed, rivers = createCopy(watershed2, rivers2)
#info = watershed_workflow.tessalateRiverAligned(watershed, rivers, river_width=widths, diagnostics=True)
m2, areas, dists = watershed_workflow.tessalateRiverAligned(watershed, rivers, river_width=widths,
refine_min_angle=min_angle, refine_distance=[d0,A0,d1,A1],
tol=0.1, diagnostics=True)
2026-02-09 11:42:21,601 - root - INFO:
2026-02-09 11:42:21,601 - root - INFO: Stream-aligned Meshing
2026-02-09 11:42:21,601 - root - INFO: ------------------------------
2026-02-09 11:42:21,601 - root - INFO: Creating stream-aligned mesh...
2026-02-09 11:42:21,612 - root - INFO: Adjusting HUC to match reaches at outlet
2026-02-09 11:42:21,617 - root - INFO: is nonoverlapping? total_area = 4621.036555146999, summed_area = 4621.036555146998
2026-02-09 11:42:21,617 - root - INFO:
2026-02-09 11:42:21,617 - root - INFO: Triangulation
2026-02-09 11:42:21,617 - root - INFO: ------------------------------
2026-02-09 11:42:21,627 - root - INFO: Triangulating...
2026-02-09 11:42:21,628 - root - INFO: 277 points and 278 facets
2026-02-09 11:42:21,628 - root - INFO: checking graph consistency
2026-02-09 11:42:21,628 - root - INFO: tolerance is set to 0.1
2026-02-09 11:42:21,629 - root - INFO: building graph data structures
2026-02-09 11:42:21,630 - root - INFO: triangle.build...
2026-02-09 11:42:25,194 - root - INFO: ...built: 868 mesh points and 1459 triangles
2026-02-09 11:42:25,194 - root - INFO: Splitting stream triangles
2026-02-09 11:42:25,210 - root - INFO: Split 14 stream triangles
2026-02-09 11:42:25,210 - root - INFO: Plotting triangulation diagnostics
2026-02-09 11:42:25,233 - root - INFO: min area = 14.0625
2026-02-09 11:42:25,234 - root - INFO: max area = 865.1326778224702
# plot the mesh
fig, ax = plt.subplots(1,1)
plot(watershed, rivers, ax=ax)
df_mesh = gpd.GeoDataFrame(geometry=[shapely.geometry.Polygon(m2.coords[c, :]) for c in m2.conn])
df_mesh.boundary.plot(ax=ax, color='grey')
<Axes: >
mls = shapely.geometry.MultiLineString([reach.linestring for river in rivers for reach in river.preOrder()])
dist = np.array([shapely.geometry.Point(coord).distance(mls) for coord in m2.coords])
# providing elevations to each point of the mesh (real case this would come from DEMs)
points3=np.zeros((len(m2.coords),3))
points3[:,:2]=m2.coords
points3[:,2]=9+dist/200+m2.coords[:,1]/400 # grad
m2.coords = points3
# depress the river corridor to burn in streams
def conditionRiverMesh(m2, river, depress_by):
m2_copy = copy.deepcopy(m2)
for reach in river:
for i, elem in enumerate(reach['elems']):
for j in range(len(elem)):
m2.coords[elem[j]][2] = m2_copy.coords[elem[j]][2] - depress_by
for river in rivers:
conditionRiverMesh(m2, river, 1)
# plotting surface mesh with elevations
fig, ax = plt.subplots(1,1)
# Get elevation data range
elev_min = m2.coords[:,2].min()
elev_max = m2.coords[:,2].max()
# Create normalized colormap
norm = plt.Normalize(vmin=elev_min, vmax=elev_max)
mp = m2.plot(facecolors='elevation', ax=ax, linewidth=0.5, cmap='terrain', norm = norm)
cbar = fig.colorbar(mp, orientation="horizontal")
ax.set_title('surface mesh with elevations')
ax.set_aspect('equal', 'datalim')
plt.title('Surface Mesh with Elevations')
plt.show()
Mesh Extrusion#
## **** Creating a 3D mesh
# this is oversimplified case where we wille xtrude this mesh in 3D with 10 layers of prescribed thicknesses
total_thickness = 10
dzs=[0.1,0.2, 0.2, 0.25,0.25,0.25,0.25,1,1.25,1.25,2,3]
print(sum(dzs))
assert(sum(dzs)==total_thickness)
# layer extrusion
# -- data structures needed for extrusion
layer_types = []
layer_data = []
layer_ncells = []
layer_mat_ids = []
depth = 0
for dz in dzs:
depth += 0.5 * dz
layer_types.append('constant')
layer_data.append(dz)
layer_ncells.append(1)
layer_mat_ids.append(1000)
depth += 0.5 * dz
# print the summary
watershed_workflow.mesh.Mesh3D.summarizeExtrusion(layer_types, layer_data,
layer_ncells, layer_mat_ids)
2026-02-09 11:42:25,507 - root - INFO: Cell summary:
2026-02-09 11:42:25,508 - root - INFO: ------------------------------------------------------------
2026-02-09 11:42:25,508 - root - INFO: l_id | c_id |mat_id | dz | z_top
2026-02-09 11:42:25,508 - root - INFO: ------------------------------------------------------------
2026-02-09 11:42:25,508 - root - INFO: 00 | 00 | 1000 | 0.100000 | 0.000000
2026-02-09 11:42:25,509 - root - INFO: 01 | 01 | 1000 | 0.200000 | 0.100000
2026-02-09 11:42:25,509 - root - INFO: 02 | 02 | 1000 | 0.200000 | 0.300000
2026-02-09 11:42:25,509 - root - INFO: 03 | 03 | 1000 | 0.250000 | 0.500000
2026-02-09 11:42:25,509 - root - INFO: 04 | 04 | 1000 | 0.250000 | 0.750000
2026-02-09 11:42:25,509 - root - INFO: 05 | 05 | 1000 | 0.250000 | 1.000000
2026-02-09 11:42:25,509 - root - INFO: 06 | 06 | 1000 | 0.250000 | 1.250000
2026-02-09 11:42:25,509 - root - INFO: 07 | 07 | 1000 | 1.000000 | 1.500000
2026-02-09 11:42:25,509 - root - INFO: 08 | 08 | 1000 | 1.250000 | 2.500000
2026-02-09 11:42:25,510 - root - INFO: 09 | 09 | 1000 | 1.250000 | 3.750000
2026-02-09 11:42:25,510 - root - INFO: 10 | 10 | 1000 | 2.000000 | 5.000000
2026-02-09 11:42:25,510 - root - INFO: 11 | 11 | 1000 | 3.000000 | 7.000000
10.0
# extrude
m3 = watershed_workflow.mesh.Mesh3D.extruded_Mesh2D(m2,layer_types, layer_data,
layer_ncells, layer_mat_ids)
Mesh IO#
Save the mesh file as an Exodus II file.
# saving mesh as exodus file
if os.path.isfile('demo_mixed_element_mesh.exo'):
os.remove('demo_mixed_element_mesh.exo')
m3.writeExodus('demo_mixed_element_mesh.exo', 'material id')
2026-02-09 11:42:25,631 - root - INFO: adding side set: 1
2026-02-09 11:42:25,633 - root - INFO: adding side set: 2
2026-02-09 11:42:25,634 - root - INFO: adding side set: 3
You are using exodus.py v 1.21.6 (seacas-py3), a python wrapper of some of the exodus library.
Copyright (c) 2013-2023 National Technology &
Engineering Solutions of Sandia, LLC (NTESS). Under the terms of
Contract DE-NA0003525 with NTESS, the U.S. Government retains certain
rights in this software.
Opening exodus file: demo_mixed_element_mesh.exo
Closing exodus file: demo_mixed_element_mesh.exo