Mixed-Element Mesh for the Coweeta Watershed#
This workflow provides a complete working example to develop an streamaligned mixed-element mesh for Coweeta watershed. Long quad elements with pentagons at junctions are placed along NHDPlus flowlines to represent rivers/streams. Rest of the domain is meshed with standard TIN.
It uses the following datasets:
NHD Plusfor the watershed boundary and hydrography.NEDfor elevationNLCDfor land cover/transpiration/rooting depthsGLYHMPSgeology data for structural formationsSoilGrids 2017for depth to bedrock and soil texture informationSSURGOfor soil data, where available, in the top 2m.
This workflow creates the following files:
Mesh file:
Coweeta.exo, includes all labeled sets
# these can be turned on for development work
%load_ext autoreload
%autoreload 2
%matplotlib ipympl
## FIX ME -- why is this broken without importing netcdf first?
import netCDF4
# setting up logging first or else it gets preempted by another package
import watershed_workflow.ui
watershed_workflow.ui.setup_logging(1)
import os,sys
import logging
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm as pcm
import shapely
import pandas as pd
import geopandas as gpd
import cftime, datetime
pd.options.display.max_columns = None
import watershed_workflow
import watershed_workflow.config
import watershed_workflow.sources
import watershed_workflow.utils
import watershed_workflow.plot
import watershed_workflow.mesh
import watershed_workflow.regions
import watershed_workflow.meteorology
import watershed_workflow.land_cover_properties
import watershed_workflow.resampling
import watershed_workflow.condition
import watershed_workflow.io
import watershed_workflow.sources.standard_names as names
# set the default figure size for notebooks
plt.rcParams["figure.figsize"] = (8, 6)
Input: Parameters and other source data#
Note, this section will need to be modified for other runs of this workflow in other regions.
# Force Watershed Workflow to pull data from this directory rather than a shared data directory.
# This picks up the Coweeta-specific datasets set up here to avoid large file downloads for
# demonstration purposes.
#
def splitPathFull(path):
"""
Splits an absolute path into a list of components such that
os.path.join(*splitPathFull(path)) == path
"""
parts = []
while True:
head, tail = os.path.split(path)
if head == path: # root on Unix or drive letter with backslash on Windows (e.g., C:\)
parts.insert(0, head)
break
elif tail == path: # just a single file or directory
parts.insert(0, tail)
break
else:
parts.insert(0, tail)
path = head
return parts
cwd = splitPathFull(os.getcwd())
# REMOVE THIS PORTION OF THE CELL for general use outside of Coweeta -- this is just locating
# the working directory within the WW directory structure
if cwd[-1] == 'Coweeta':
pass
elif cwd[-1] == 'examples':
cwd.append('Coweeta')
else:
cwd.extend(['examples','Coweeta'])
# END REMOVE THIS PORTION
# Note, this directory is where downloaded data will be put as well
data_dir = os.path.join(*(cwd + ['input_data',]))
def toInput(filename):
return os.path.join(data_dir, filename)
output_dir = os.path.join(*(cwd + ['output_data',]))
def toOutput(filename):
return os.path.join(output_dir, filename)
work_dir = os.path.join(*cwd)
def toWorkingDir(filename):
return os.path.join(work_dir, filename)
# Set the data directory to the local space to get the locally downloaded files
# REMOVE THIS CELL for general use outside fo Coweeta
watershed_workflow.config.setDataDirectory(data_dir)
## Parameters cell -- this provides all parameters that can be changed via pipelining to generate a new watershed.
name = 'Coweeta'
coweeta_shapefile = './Coweeta/input_data/coweeta_basin.shp'
# Geometric parameters
# -- parameters to clean and reduce the river network prior to meshing
simplify = 60 # length scale to target average edge
ignore_small_rivers = 2 # remove rivers with fewer than this number of reaches -- important for NHDPlus HR
prune_by_area_fraction = 0.01 # prune any reaches whose contributing area is less than this fraction of the domain
# -- mesh triangle refinement control
refine_d0 = 200
refine_d1 = 600
#refine_L0 = 75
#refine_L1 = 200
refine_L0 = 125
refine_L1 = 300
refine_A0 = refine_L0**2 / 2
refine_A1 = refine_L1**2 / 2
min_angle = 32 # degrees
# Simulation control
# - note that we use the NoLeap calendar, same as DayMet. Simulations are typically run over the "water year"
# which starts August 1.
start = cftime.DatetimeNoLeap(2010,8,1)
end = cftime.DatetimeNoLeap(2011,8,1)
# Global Soil Properties
min_porosity = 0.05 # minimum porosity considered "too small"
max_permeability = 1.e-10 # max value considered "too permeable"
max_vg_alpha = 1.e-3 # max value of van Genuchten's alpha -- our correlation is not valid for some soils
# a dictionary of output_filenames -- will include all filenames generated
output_filenames = {}
# Note that, by default, we tend to work in the DayMet CRS because this allows us to avoid
# reprojecting meteorological forcing datasets.
crs = watershed_workflow.crs.daymet_crs
# get the shape and crs of the shape
coweeta_source = watershed_workflow.sources.ManagerShapefile(coweeta_shapefile, id_name='BASIN_CODE')
coweeta = coweeta_source.getShapes(out_crs=crs)
2026-03-09 16:30:42,990 - root - INFO: fixing column: geometry
# set up a dictionary of source objects
#
# Data sources, also called managers, deal with downloading and parsing data files from a variety of online APIs.
sources = watershed_workflow.sources.getDefaultSources()
sources['hydrography'] = watershed_workflow.sources.hydrography_sources['NHDPlus HR']
#
# This demo uses a few datasets that have been clipped out of larger, national
# datasets and are distributed with the code. This is simply to save download
# time for this simple problem and to lower the barrier for trying out
# Watershed Workflow. A more typical workflow would delete these lines (as
# these files would not exist for other watersheds).
#
# The default versions of these download large raster and shapefile files that
# are defined over a very large region (globally or the entire US).
#
# DELETE THIS SECTION for non-Coweeta runs
dtb_file = os.path.join(data_dir, 'soil_structure', 'DTB', 'DTB.tif')
geo_file = os.path.join(data_dir, 'soil_structure', 'GLHYMPS', 'GLHYMPS.shp')
# GLHYMPs is a several-GB download, so we have sliced it and included the slice here
sources['geologic structure'] = watershed_workflow.sources.ManagerGLHYMPS(geo_file)
# The Pelletier DTB map is not particularly accurate at Coweeta -- the SoilGrids map seems to be better.
# Here we will use a clipped version of that map.
sources['depth to bedrock'] = watershed_workflow.sources.ManagerRaster(dtb_file)
# END DELETE THIS SECTION
# log the sources that will be used here
watershed_workflow.sources.logSources(sources)
2026-03-09 16:30:43,009 - root - INFO: Using sources:
2026-03-09 16:30:43,009 - root - INFO: --------------
2026-03-09 16:30:43,009 - root - INFO: HUC: WBD
2026-03-09 16:30:43,009 - root - INFO: hydrography: NHDPlus HR
2026-03-09 16:30:43,010 - root - INFO: DEM: 3DEP 60m
2026-03-09 16:30:43,010 - root - INFO: soil structure: National Resources Conservation Service Soil Survey (NRCS Soils)
2026-03-09 16:30:43,010 - root - INFO: geologic structure: shapefile: "GLHYMPS.shp"
2026-03-09 16:30:43,010 - root - INFO: land cover: NLCD 2021 L48
2026-03-09 16:30:43,010 - root - INFO: LAI: MODIS
2026-03-09 16:30:43,011 - root - INFO: depth to bedrock: raster: "DTB.tif"
2026-03-09 16:30:43,011 - root - INFO: meteorology: AORC v1.1
Basin Geometry#
In this section, we choose the basin, the streams to be included in the stream-aligned mesh, and make sure that all are resolved discretely at appropriate length scales for this work.
the Watershed#
coweeta
| AREA | PERIMETER | CWTBASINNA | CWTBASIN_1 | BASIN_CODE | SPOT | LABEL | geometry | ID | name | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1.626020e+07 | 17521.768 | 2 | 1 | 1 | -9999 | Coweeta Hydrologic Lab | POLYGON ((1443453.518 -645937.256, 1443489.888... | 1 | 1 |
# Construct and plot the WW object used for storing watersheds
watershed = watershed_workflow.split_hucs.SplitHUCs(coweeta)
watershed.plot()
2026-03-09 16:30:43,051 - root - INFO: Removing holes on 1 polygons
2026-03-09 16:30:43,052 - root - INFO: -- removed interior
2026-03-09 16:30:43,052 - root - INFO: -- union
2026-03-09 16:30:43,052 - root - INFO: Parsing 1 components for holes
2026-03-09 16:30:43,053 - root - INFO: -- complete
<Axes: >
the Rivers#
# download/collect the river network within that shape's bounds
reaches = sources['hydrography'].getShapesByGeometry(watershed.exterior, crs, out_crs=crs)
rivers = watershed_workflow.river_tree.createRivers(reaches, method='hydroseq')
watershed_orig, rivers_orig = watershed, rivers
sources['hydrography'].name, sources['hydrography'].source
2026-03-09 16:30:43,121 - root - INFO: Shapes manager: buffered+snapped query box = (-83.482, 35.024, -83.418, 35.077)
2026-03-09 16:30:43,208 - root - INFO: fixing column: geometry
2026-03-09 16:30:43,212 - root - INFO: fixing column: catchment
('NHDPlus HR', 'HyRiver.NHDPlusHR')
# plot the rivers and watershed
def plot(ws, rivs, ax=None):
if ax is None:
fig, ax = plt.subplots(1, 1)
ws.plot(color='k', marker='+', markersize=10, ax=ax)
for river in rivs:
river.plot(marker='x', markersize=10, ax=ax)
plot(watershed, rivers)
# keeping the originals for plotting comparisons
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]
watershed, rivers = createCopy(watershed_orig, rivers_orig)
# simplifying -- this sets the discrete length scale of both the watershed boundary and the rivers
watershed_workflow.simplify(watershed, rivers, refine_L0, refine_L1, refine_d0, refine_d1)
# simplify may remove reaches from the rivers object
# -- this call removes any reaches from the dataframe as well, signaling we are all done removing reaches
#
# ETC: NOTE -- can this be moved into the simplify call?
for river in rivers:
river.resetDataFrame()
# Now that the river network is set, find the watershed boundary outlets
for river in rivers:
watershed_workflow.hydrography.findOutletsByCrossings(watershed, river)
2026-03-09 16:30:43,339 - root - INFO:
2026-03-09 16:30:43,339 - root - INFO: Simplifying
2026-03-09 16:30:43,339 - root - INFO: ------------------------------
2026-03-09 16:30:43,339 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,339 - root - INFO: Presimplify to remove colinear, coincident points.
2026-03-09 16:30:43,346 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,346 - root - INFO: Pruning leaf reaches < 125
2026-03-09 16:30:43,347 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,347 - root - INFO: Merging internal reaches < 125
2026-03-09 16:30:43,348 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,350 - root - INFO: reach: min seg length: 4.5823534935 min geom length: 160.2183496376
2026-03-09 16:30:43,351 - root - INFO: reach: med seg length: 5.5205081646 med geom length: 957.0878211542
2026-03-09 16:30:43,352 - root - INFO: reach: max seg length: 27.5794801882 max geom length: 3922.9535292603
2026-03-09 16:30:43,354 - root - INFO:
2026-03-09 16:30:43,355 - root - INFO: HUC : min seg length: 5.8226468933 min geom length: 16872.4086261245
2026-03-09 16:30:43,355 - root - INFO: HUC : med seg length: 19.4631658810 med geom length: 16872.4086261245
2026-03-09 16:30:43,356 - root - INFO: HUC : max seg length: 405.2969428374 max geom length: 16872.4086261245
2026-03-09 16:30:43,356 - root - INFO:
2026-03-09 16:30:43,357 - root - INFO: Snapping discrete points to make rivers and HUCs discretely consistent.
2026-03-09 16:30:43,357 - root - INFO: -- snapping HUC triple junctions to reaches
2026-03-09 16:30:43,360 - root - INFO: reach: min seg length: 4.5823534935 min geom length: 160.2183496376
2026-03-09 16:30:43,360 - root - INFO: reach: med seg length: 5.5205081646 med geom length: 957.0878211542
2026-03-09 16:30:43,360 - root - INFO: reach: max seg length: 27.5794801882 max geom length: 3922.9535292603
2026-03-09 16:30:43,361 - root - INFO:
2026-03-09 16:30:43,361 - root - INFO: HUC : min seg length: 5.8226468933 min geom length: 16872.4086261245
2026-03-09 16:30:43,364 - root - INFO: HUC : med seg length: 19.4631658810 med geom length: 16872.4086261245
2026-03-09 16:30:43,364 - root - INFO: HUC : max seg length: 405.2969428374 max geom length: 16872.4086261245
2026-03-09 16:30:43,365 - root - INFO:
2026-03-09 16:30:43,365 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,365 - root - INFO: -- snapping reach endpoints to HUC boundaries
2026-03-09 16:30:43,376 - root - INFO: reach: min seg length: 4.5823534935 min geom length: 160.2183496376
2026-03-09 16:30:43,376 - root - INFO: reach: med seg length: 5.5205081646 med geom length: 957.0878211542
2026-03-09 16:30:43,376 - root - INFO: reach: max seg length: 27.5794801882 max geom length: 3922.9535292603
2026-03-09 16:30:43,376 - root - INFO:
2026-03-09 16:30:43,377 - root - INFO: HUC : min seg length: 5.8226468933 min geom length: 16872.4086261245
2026-03-09 16:30:43,377 - root - INFO: HUC : med seg length: 19.4631658810 med geom length: 16872.4086261245
2026-03-09 16:30:43,377 - root - INFO: HUC : max seg length: 405.2969428374 max geom length: 16872.4086261245
2026-03-09 16:30:43,377 - root - INFO:
2026-03-09 16:30:43,378 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,378 - root - INFO: -- cutting reaches at HUC boundaries
2026-03-09 16:30:43,378 - root - INFO: intersection found
/Users/Shared/ornldev/code/watershed_workflow/repos/master/watershed_workflow/utils.py:832: RuntimeWarning: invalid value encountered in divide
ds = ds / np.linalg.norm(ds)
2026-03-09 16:30:43,380 - root - INFO: - cutting reach at external boundary of HUCs:
2026-03-09 16:30:43,380 - root - INFO: split HUC boundary ls into 2 pieces
2026-03-09 16:30:43,380 - root - INFO: split reach ls into 2 pieces
2026-03-09 16:30:43,389 - root - INFO: reach: min seg length: 1.9676019217 min geom length: 160.2183496376
2026-03-09 16:30:43,389 - root - INFO: reach: med seg length: 5.5225100483 med geom length: 904.8489049176
2026-03-09 16:30:43,390 - root - INFO: reach: max seg length: 27.5794801882 max geom length: 3922.9535292603
2026-03-09 16:30:43,390 - root - INFO:
2026-03-09 16:30:43,390 - root - INFO: HUC : min seg length: 5.8226468933 min geom length: 4467.6054418386
2026-03-09 16:30:43,390 - root - INFO: HUC : med seg length: 19.4623135569 med geom length: 8436.2043130623
2026-03-09 16:30:43,390 - root - INFO: HUC : max seg length: 405.2969428374 max geom length: 12404.8031842860
2026-03-09 16:30:43,391 - root - INFO:
2026-03-09 16:30:43,391 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,391 - root - INFO:
2026-03-09 16:30:43,391 - root - INFO: Simplification Diagnostics
2026-03-09 16:30:43,391 - root - INFO: ------------------------------
2026-03-09 16:30:43,393 - root - INFO: reach: min seg length: 1.9676019217 min geom length: 160.2183496376
2026-03-09 16:30:43,393 - root - INFO: reach: med seg length: 5.5225100483 med geom length: 904.8489049176
2026-03-09 16:30:43,393 - root - INFO: reach: max seg length: 27.5794801882 max geom length: 3922.9535292603
2026-03-09 16:30:43,394 - root - INFO:
2026-03-09 16:30:43,394 - root - INFO: HUC : min seg length: 5.8226468933 min geom length: 4467.6054418386
2026-03-09 16:30:43,394 - root - INFO: HUC : med seg length: 19.4623135569 med geom length: 8436.2043130623
2026-03-09 16:30:43,394 - root - INFO: HUC : max seg length: 405.2969428374 max geom length: 12404.8031842860
2026-03-09 16:30:43,395 - root - INFO:
2026-03-09 16:30:43,395 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,395 - root - INFO:
2026-03-09 16:30:43,395 - root - INFO: Resampling HUC and river
2026-03-09 16:30:43,395 - root - INFO: ------------------------------
2026-03-09 16:30:43,396 - root - INFO: -- resampling HUCs based on distance function (200, 125, 600, 300)
2026-03-09 16:30:43,443 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,443 - root - INFO: -- resampling reaches based on uniform target 125
2026-03-09 16:30:43,449 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,449 - root - INFO:
2026-03-09 16:30:43,449 - root - INFO: Resampling Diagnostics
2026-03-09 16:30:43,449 - root - INFO: ------------------------------
2026-03-09 16:30:43,450 - root - INFO: reach: min seg length: 78.3369140284 min geom length: 157.8784571387
2026-03-09 16:30:43,450 - root - INFO: reach: med seg length: 117.7932453013 med geom length: 893.9653921062
2026-03-09 16:30:43,451 - root - INFO: reach: max seg length: 124.2541375207 max geom length: 3870.3355820777
2026-03-09 16:30:43,451 - root - INFO:
2026-03-09 16:30:43,451 - root - INFO: HUC : min seg length: 120.9289139295 min geom length: 4224.3783428488
2026-03-09 16:30:43,451 - root - INFO: HUC : med seg length: 212.8777438811 med geom length: 7933.7936132856
2026-03-09 16:30:43,452 - root - INFO: HUC : max seg length: 294.9177323115 max geom length: 11643.2088837225
2026-03-09 16:30:43,452 - root - INFO:
2026-03-09 16:30:43,452 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,452 - root - INFO:
2026-03-09 16:30:43,452 - root - INFO: Clean up sharp angles, both internally and at junctions.
2026-03-09 16:30:43,453 - root - INFO: ------------------------------
2026-03-09 16:30:43,454 - root - INFO: ... cleaned up 0 sharp angles on HUCs.
2026-03-09 16:30:43,456 - root - INFO: ... cleaned up 0 sharp angles at outlets on river 0.
2026-03-09 16:30:43,477 - root - INFO: ... cleaned up 0 internal sharp angles on river 0.
2026-03-09 16:30:43,477 - root - INFO: Cleaned up 0 sharp angles.
2026-03-09 16:30:43,478 - root - INFO: reach: min seg length: 78.3369140284 min geom length: 157.8784571387
2026-03-09 16:30:43,478 - root - INFO: reach: med seg length: 117.7932453013 med geom length: 893.9653921062
2026-03-09 16:30:43,478 - root - INFO: reach: max seg length: 124.2541375207 max geom length: 3870.3355820777
2026-03-09 16:30:43,478 - root - INFO:
2026-03-09 16:30:43,479 - root - INFO: HUC : min seg length: 120.9289139295 min geom length: 4224.3783428488
2026-03-09 16:30:43,479 - root - INFO: HUC : med seg length: 212.8777438811 med geom length: 7933.7936132856
2026-03-09 16:30:43,479 - root - INFO: HUC : max seg length: 294.9177323115 max geom length: 11643.2088837225
2026-03-09 16:30:43,479 - root - INFO:
2026-03-09 16:30:43,479 - root - INFO: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs +type=crs
2026-03-09 16:30:43,519 - root - INFO: Crossings by Polygon:
2026-03-09 16:30:43,519 - root - INFO: Polygon 0
2026-03-09 16:30:43,519 - root - INFO: crossing: [1446635.61839213 -646208.54462176]
2026-03-09 16:30:43,520 - root - INFO: Constructing outlet list
2026-03-09 16:30:43,520 - root - INFO: last outlet is 0 in polygon 0 at {crossings_clusters_centroids[last_outlet]}
plot(watershed, rivers)
watershed.df.columns
Index(['AREA', 'PERIMETER', 'CWTBASINNA', 'CWTBASIN_1', 'BASIN_CODE', 'SPOT',
'LABEL', 'geometry', 'ID', 'name', 'outlet'],
dtype='str')
# this generates a zoomable map, showing different reaches and watersheds,
# with discrete points. Problem areas are clickable to get IDs for manual
# modifications.
m = watershed.explore(marker=False)
for river in rivers_orig:
m = river.explore(m=m, column=None, color='black', name=river['name']+' raw', marker=False)
for river in rivers:
m = river.explore(m=m)
m = watershed_workflow.makeMap(m)
m
KWARGS:
{'tooltip': False, 'popup': ['ID', 'name', 'AREA', 'PERIMETER', 'CWTBASINNA', 'CWTBASIN_1', 'BASIN_CODE', 'SPOT'], 'legend': True, 'style_kwds': {'weight': 5, 'fillOpacity': 0.2}, 'cmap': <matplotlib.colors.ListedColormap object at 0x32152f680>, 'vmin': np.int32(1), 'vmax': np.int32(1), 'highlight_kwds': {'fillOpacity': 0.4}}
Make this Notebook Trusted to load map: File -> Trust Notebook
Mesh Geometry#
Discretely create the stream-aligned mesh. Download elevation data, and condition the mesh discretely to make for better topography.
# Refine triangles if they get too acute
# width of reach by stream order (order:width)
def widths(reach):
mapping = {1: 8, 2: 12, 3: 16}
order = reach.properties['stream_order']
return mapping.get(order, 8) # default width if order not found
# create the mesh
m2, areas, dists = watershed_workflow.tessalateRiverAligned(watershed, rivers,
river_width=widths,
refine_min_angle=min_angle,
refine_distance=[refine_d0, refine_A0, refine_d1, refine_A1],
diagnostics=True)
2026-03-09 16:30:44,516 - root - INFO:
2026-03-09 16:30:44,516 - root - INFO: Stream-aligned Meshing
2026-03-09 16:30:44,517 - root - INFO: ------------------------------
2026-03-09 16:30:44,517 - root - INFO: Creating stream-aligned mesh...
2026-03-09 16:30:44,545 - root - INFO: ... created a mesh with 184 elements for river 0
2026-03-09 16:30:44,546 - root - INFO: Adjusting HUC to match reaches at outlet
2026-03-09 16:30:44,555 - root - INFO: is nonoverlapping? total_area = 198327.62368762706, summed_area = 198327.6236876241
2026-03-09 16:30:44,555 - root - INFO: Building the remaining triangular mesh...
2026-03-09 16:30:44,556 - root - INFO:
2026-03-09 16:30:44,556 - root - INFO: Triangulation
2026-03-09 16:30:44,556 - root - INFO: ------------------------------
2026-03-09 16:30:44,570 - root - INFO: Triangulating...
2026-03-09 16:30:44,571 - root - INFO: 443 points and 444 facets
2026-03-09 16:30:44,571 - root - INFO: checking graph consistency
2026-03-09 16:30:44,571 - root - INFO: tolerance is set to 1.0
2026-03-09 16:30:44,572 - root - INFO: building graph data structures
2026-03-09 16:30:44,573 - root - INFO: triangle.build...
2026-03-09 16:30:51,071 - root - INFO: ...built: 1427 mesh points and 2409 triangles
2026-03-09 16:30:51,072 - root - INFO: Plotting triangulation diagnostics
2026-03-09 16:30:51,111 - root - INFO: min area = 2570.255859375
2026-03-09 16:30:51,112 - root - INFO: max area = 39949.32360839844
2026-03-09 16:30:51,129 - root - INFO: Fixing corridor-spanning triangles...
2026-03-09 16:30:51,261 - root - INFO: ... split 9 spanning or junction triangles
# get a raster for the elevation map, based on 3DEP
dem = sources['DEM'].getDataset(watershed.exterior, watershed.crs)['dem']
# provide surface mesh elevations
watershed_workflow.elevate(m2, dem, method='linear')
# also elevate the river network linestrings
watershed_workflow.condition.setProfileByDEM(rivers, dem)
2026-03-09 16:30:51,325 - root - INFO: Incoming shape area = 0.0016041250506451986
2026-03-09 16:30:51,325 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.00324
2026-03-09 16:30:51,325 - root - INFO: ... buffered shape area = 0.0021711505755106406
2026-03-09 16:30:51,325 - root - INFO: ... snapped bounding box = (-83.48184, 35.0244, -83.41812, 35.07732)
2026-03-09 16:30:51,325 - root - INFO: Getting DEM with map of area = 0.0021711505755106406
# Plot the DEM raster
fig, ax = plt.subplots(1,1)
# Plot the DEM data
im = dem.plot(ax=ax, cmap='terrain', add_colorbar=False)
# Add colorbar
cbar = plt.colorbar(im, ax=ax, shrink=0.8)
cbar.set_label('Elevation (m)', rotation=270, labelpad=15)
# Add title and labels
ax.set_title('Digital Elevation Model (DEM)', fontsize=14, fontweight='bold')
ax.set_xlabel('X Coordinate')
ax.set_ylabel('Y Coordinate')
# Set equal aspect ratio
ax.set_aspect('equal')
plt.tight_layout()
plt.show()
There are a range of options to condition river corridor mesh. We hydrologically condition the river mesh, ensuring unimpeded water flow in river corridors by globally adjusting flowlines to rectify artificial obstructions from inconsistent DEM elevations or misalignments. Please read the documentation for more information
In the pit-filling algorithm, we want to make sure that river corridor is not filled up. Hence we exclude river corridor cells from the pit-filling algorithm.
# now condition the river to fix mis-hits, where the corridor centroids do not fall in the DEM's idea of the river, enforcing monotonicity of the river network
def computeBurnInDepth(da_sq_miles):
"""burn-in depth as a function of drainage area"""
depth_in_feet = 1.22 * da_sq_miles**0.317
return 0.3048 * depth_in_feet # ft --> meters
def computeBurnInDepthFromReach(reach):
depth = computeBurnInDepth(reach['drainage_area_sqkm'] * 0.386102)
logging.debug(f"reach of DA {reach['drainage_area_sqkm']} has depth {depth}")
return depth
watershed_workflow.condition.conditionRiverMeshes(m2,
rivers,
network_burn_in_depth=computeBurnInDepthFromReach)
# hydrologically condition the non-corridor portion of the mesh, removing pits
outlet_edge = watershed_workflow.mesh.Edge(rivers[0]['elems'][-1][0], rivers[0]['elems'][-1][-1])
preserved_pits = [c for (c,conn) in enumerate(m2.conn) if len(conn) > 3]
m2r, res = watershed_workflow.condition.conditionMesh(m2,
preserved_pits=preserved_pits,
forced_outlet_edges=[outlet_edge,],
epsilon = 0.01
)
2026-03-09 16:30:52,676 - root - INFO:
2026-03-09 16:30:52,677 - root - INFO: Running marching_iterative: 5 initial pits
2026-03-09 16:30:52,677 - root - INFO: ==============================================================================
2026-03-09 16:30:52,832 - root - INFO: ... iteration 0 of marching: 0 pits
2026-03-09 16:30:52,832 - root - INFO: ... done iterating in 1 iterations, 0 pits.
2026-03-09 16:30:52,864 - root - INFO: completed: 0 final pits, 5 removed
2026-03-09 16:30:52,864 - root - INFO: RMSE of dz: 14.27023393620527
2026-03-09 16:30:52,865 - root - INFO: MAE of dz: 2.1890246643363964
2026-03-09 16:30:52,865 - root - INFO: MAX of dz: 217.65307542325263
# plotting surface mesh with elevations
fig, ax = plt.subplots(1,1)
ax2 = ax.inset_axes([0.65,0.05,0.3,0.5])
mp = m2.plot(facecolors='elevation', edgecolors=None, ax=ax, linewidth=0.5, colorbar=False)
cbar = fig.colorbar(mp, orientation="horizontal")
ax.set_title('surface mesh with elevations')
ax.set_aspect('equal', 'datalim')
mp2 = m2.plot(facecolors='elevation', edgecolors='white', ax=ax2, colorbar=False)
ax2.set_aspect('equal', 'datalim')
xlim = (1.4433e6, 1.4438e6)
ylim = (-647000, -647500)
ax2.set_xlim(xlim)
ax2.set_ylim(ylim)
ax2.set_xticks([])
ax2.set_yticks([])
ax.indicate_inset_zoom(ax2, edgecolor='k')
cbar.ax.set_title('elevation [m]')
plt.show()
# [1445055.35066667 -646519.51766667]
2026-03-09 16:30:52,957 - matplotlib.axes._base - WARNING: Ignoring fixed x limits to fulfill fixed data aspect with adjustable data limits.
# add labeled sets for subcatchments and outlets
watershed_workflow.regions.addWatershedAndOutletRegions(m2, watershed, outlet_width=250, exterior_outlet=True)
# add labeled sets for river corridor cells
watershed_workflow.regions.addRiverCorridorRegions(m2, rivers)
# add labeled sets for river corridor cells by order
watershed_workflow.regions.addStreamOrderRegions(m2, rivers)
2026-03-09 16:30:53,057 - root - INFO: Adding regions for 1 polygons
2026-03-09 16:30:53,087 - root - INFO: Exterior outlet point (from attribute): POINT (1446635.6183921252 -646208.5446217591)
for ls in m2.labeled_sets:
print(f'{ls.setid} : {ls.entity} : {len(ls.ent_ids)} : "{ls.name}"')
10000 : CELL : 2611 : "1"
10001 : CELL : 2611 : "1 surface"
10002 : FACE : 76 : "1 boundary"
10003 : FACE : 5 : "1 outlet"
10004 : FACE : 5 : "surface domain outlet"
10005 : CELL : 184 : "river corridor 0 surface"
10006 : CELL : 16 : "stream order 3"
10007 : CELL : 46 : "stream order 2"
10008 : CELL : 122 : "stream order 1"
Surface properties#
Meshes interact with data to provide forcing, parameters, and more in the actual simulation. Specifically, we need vegetation type on the surface to provide information about transpiration and subsurface structure to provide information about water retention curves, etc.
NLCD for LULC#
We’ll start by downloading and collecting land cover from the NLCD dataset, and generate sets for each land cover type that cover the surface. Likely these will be some combination of grass, deciduous forest, coniferous forest, and mixed.
# download the NLCD raster
nlcd = sources['land cover'].getDataset(watershed.exterior.buffer(100), watershed.crs)['cover']
# what land cover types did we get?
logging.info('Found land cover dtypes: {}'.format(nlcd.dtype))
logging.info('Found land cover types: {}'.format(set(list(nlcd.values.ravel()))))
2026-03-09 16:30:53,234 - root - INFO: Incoming shape area = 0.001776229723905591
2026-03-09 16:30:53,235 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.00081
2026-03-09 16:30:53,235 - root - INFO: ... buffered shape area = 0.0019172956875437115
2026-03-09 16:30:53,235 - root - INFO: ... snapped bounding box = (-83.48022, 35.02629, -83.42028, 35.0757)
2026-03-09 16:30:53,279 - root - INFO: Found land cover dtypes: uint8
2026-03-09 16:30:53,280 - root - INFO: Found land cover types: {np.uint8(71), np.uint8(41), np.uint8(42), np.uint8(43), np.uint8(81), np.uint8(52), np.uint8(21), np.uint8(22), np.uint8(23), np.uint8(127)}
# create a colormap for the data
nlcd_indices, nlcd_cmap, nlcd_norm, nlcd_ticks, nlcd_labels = \
watershed_workflow.colors.createNLCDColormap(np.unique(nlcd))
nlcd_cmap
from_list
under
bad
over
fig, ax = plt.subplots(1,1)
nlcd.plot.imshow(ax=ax, cmap=nlcd_cmap, norm=nlcd_norm, add_colorbar=False)
watershed_workflow.colors.createIndexedColorbar(ncolors=len(nlcd_indices),
cmap=nlcd_cmap, labels=nlcd_labels, ax=ax)
ax.set_title('Land Cover')
plt.show()
# map nlcd onto the mesh
m2_nlcd = watershed_workflow.getDatasetOnMesh(m2, nlcd, method='nearest')
m2.cell_data['land_cover'] = m2_nlcd
# double-check that nan not in the values
assert 127 not in m2_nlcd
# create a new set of labels and indices with only those that actually appear on the mesh
nlcd_indices, nlcd_cmap, nlcd_norm, nlcd_ticks, nlcd_labels = \
watershed_workflow.colors.createNLCDColormap(np.unique(m2_nlcd))
mp = m2.plot(facecolors=m2_nlcd, cmap=nlcd_cmap, norm=nlcd_norm, edgecolors=None, colorbar=False)
watershed_workflow.colors.createIndexedColorbar(ncolors=len(nlcd_indices),
cmap=nlcd_cmap, labels=nlcd_labels, ax=plt.gca())
plt.show()
# add labeled sets to the mesh for NLCD
nlcd_labels_dict = dict(zip(nlcd_indices, nlcd_labels))
watershed_workflow.regions.addSurfaceRegions(m2, names=nlcd_labels_dict)
nlcd_labels_dict
{np.uint8(21): 'Developed, Open Space',
np.uint8(22): 'Developed, Low Intensity',
np.uint8(23): 'Developed, Medium Intensity',
np.uint8(41): 'Deciduous Forest',
np.uint8(42): 'Evergreen Forest',
np.uint8(43): 'Mixed Forest',
np.uint8(81): 'Pasture/Hay'}
for ls in m2.labeled_sets:
print(f'{ls.setid} : {ls.entity} : {len(ls.ent_ids)} : "{ls.name}"')
10000 : CELL : 2611 : "1"
10001 : CELL : 2611 : "1 surface"
10002 : FACE : 76 : "1 boundary"
10003 : FACE : 5 : "1 outlet"
10004 : FACE : 5 : "surface domain outlet"
10005 : CELL : 184 : "river corridor 0 surface"
10006 : CELL : 16 : "stream order 3"
10007 : CELL : 46 : "stream order 2"
10008 : CELL : 122 : "stream order 1"
21 : CELL : 105 : "Developed, Open Space"
22 : CELL : 3 : "Developed, Low Intensity"
23 : CELL : 1 : "Developed, Medium Intensity"
41 : CELL : 1350 : "Deciduous Forest"
42 : CELL : 42 : "Evergreen Forest"
43 : CELL : 1104 : "Mixed Forest"
81 : CELL : 6 : "Pasture/Hay"
MODIS LAI#
Leaf area index is needed on each land cover type – this is used in the Evapotranspiration calculation.
# download LAI and corresponding LULC datasets -- these are actually already downloaded,
# as the MODIS AppEEARS API is quite slow
#
# Note that MODIS does NOT work with the noleap calendar, so we have to convert to actual dates first
start_leap = cftime.DatetimeGregorian(start.year, start.month, start.day)
end_leap = cftime.DatetimeGregorian(end.year, end.month, end.day)
modis_data = sources['LAI'].getDataset(watershed.exterior, crs, start_leap, end_leap)
2026-03-09 16:30:53,721 - root - INFO: Incoming shape area = 0.0016041250506451986
2026-03-09 16:30:53,722 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.013500000000000002
2026-03-09 16:30:53,722 - root - INFO: ... buffered shape area = 0.004353081553256633
2026-03-09 16:30:53,722 - root - INFO: ... snapped bounding box = (-83.49300000000001, 35.014500000000005, -83.40750000000001, 35.091)
2026-03-09 16:30:53,722 - root - INFO: Building MODIS request for bounds: (-83.49300000000001, 35.014500000000005, -83.40750000000001, 35.091), years 2010-2011
2026-03-09 16:30:53,723 - root - INFO: Using superset cache for LAI: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/input_data/land_cover/MODIS/MODIS_LAI_2010-2011_-83.4930_35.0145_-83.4075_35.0910.nc
2026-03-09 16:30:53,723 - root - INFO: Using superset cache for LULC: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/input_data/land_cover/MODIS/MODIS_LULC_2010-2011_-83.4930_35.0145_-83.4075_35.0910.nc
2026-03-09 16:30:53,723 - root - INFO: Cache filenames:
2026-03-09 16:30:53,723 - root - INFO: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/input_data/land_cover/MODIS/MODIS_LAI_2010-2011_-83.4930_35.0145_-83.4075_35.0910.nc
2026-03-09 16:30:53,723 - root - INFO: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/input_data/land_cover/MODIS/MODIS_LULC_2010-2011_-83.4930_35.0145_-83.4075_35.0910.nc
2026-03-09 16:30:53,723 - root - INFO: All files exist locally.
assert modis_data['LAI'].rio.crs is not None
print(modis_data['LULC'].rio.crs, modis_data['LULC'].dtype)
EPSG:4269 float64
# MODIS data comes with time-dependent LAI AND time-dependent LULC -- just take the mode to find the most common LULC
modis_data['LULC'] = watershed_workflow.data.computeMode(modis_data['LULC'], 'time_LULC')
# now it is safe to have only one time
modis_data = modis_data.rename({'time_LAI':'time'})
# remove leap day (366th day of any leap year) to match our Noleap Calendar
modis_data = watershed_workflow.data.filterLeapDay(modis_data)
# plot the MODIS data -- note the entire domain is covered with one type for Coweeta (it is small!)
fig, axs = plt.subplots(1,2, figsize=(12,6))
modis_data['LULC'].plot.imshow(ax=axs[0])
modis_data['LAI'][0].plot.imshow(ax=axs[1])
<matplotlib.image.AxesImage at 0x375e67980>
# compute the transient time series
modis_lai = watershed_workflow.land_cover_properties.computeTimeSeries(modis_data['LAI'], modis_data['LULC'],
polygon=watershed.exterior, polygon_crs=watershed.crs)
modis_lai
| Deciduous Broadleaf Forests LAI [-] | time | |
|---|---|---|
| 0 | 3.243011 | 2010-08-01 00:00:00 |
| 1 | 4.847312 | 2010-08-05 00:00:00 |
| 2 | 3.476344 | 2010-08-09 00:00:00 |
| 3 | 4.193548 | 2010-08-13 00:00:00 |
| 4 | 3.186022 | 2010-08-17 00:00:00 |
| ... | ... | ... |
| 84 | 4.968817 | 2011-06-30 00:00:00 |
| 85 | 3.388172 | 2011-07-04 00:00:00 |
| 86 | 4.444086 | 2011-07-08 00:00:00 |
| 87 | 2.693548 | 2011-07-12 00:00:00 |
| 88 | 6.363441 | 2011-07-16 00:00:00 |
89 rows × 2 columns
# smooth the data in time
modis_lai_smoothed = watershed_workflow.data.smoothTimeSeries(modis_lai, 'time')
# save the MODIS time series to disk
output_filenames['modis_lai_transient'] = toOutput(f'{name}_LAI_MODIS_transient.h5')
watershed_workflow.io.writeTimeseriesToHDF5(output_filenames['modis_lai_transient'], modis_lai_smoothed)
watershed_workflow.land_cover_properties.plotLAI(modis_lai_smoothed, indices='MODIS')
2026-03-09 16:30:54,109 - root - INFO: Writing HDF5 file: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/output_data/Coweeta_LAI_MODIS_transient.h5
# compute a typical year
td = datetime.timedelta(days=365 * 10)
modis_lai_typical = watershed_workflow.data.computeAverageYear(modis_lai_smoothed,
start_date=(start - td),
output_nyears=10)
output_filenames['modis_lai_cyclic_steadystate'] = toOutput(f'{name}_LAI_MODIS_CyclicSteadystate.h5')
watershed_workflow.io.writeTimeseriesToHDF5(output_filenames['modis_lai_cyclic_steadystate'], modis_lai_typical)
watershed_workflow.land_cover_properties.plotLAI(modis_lai_typical, indices='MODIS')
2026-03-09 16:30:54,175 - root - INFO: Writing HDF5 file: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/output_data/Coweeta_LAI_MODIS_CyclicSteadystate.h5
Crosswalk of LAI to NLCD LC#
crosswalk = watershed_workflow.land_cover_properties.computeCrosswalk(modis_data['LULC'], nlcd, method='fractional area')
2026-03-09 16:30:54,312 - root - INFO: Compute the crosswalk between MODIS and NLCD:
2026-03-09 16:30:54,312 - root - INFO: unique MODIS: [np.float64(4.0)]
2026-03-09 16:30:54,312 - root - INFO: unique NLCD: [np.uint8(21), np.uint8(22), np.uint8(23), np.uint8(41), np.uint8(42), np.uint8(43), np.uint8(52), np.uint8(71), np.uint8(81)]
# Compute the NLCD-based time series
nlcd_lai_cyclic_steadystate = watershed_workflow.land_cover_properties.applyCrosswalk(crosswalk, modis_lai_typical)
nlcd_lai_transient = watershed_workflow.land_cover_properties.applyCrosswalk(crosswalk, modis_lai_smoothed)
watershed_workflow.land_cover_properties.removeNullLAI(nlcd_lai_cyclic_steadystate)
watershed_workflow.land_cover_properties.removeNullLAI(nlcd_lai_transient)
nlcd_lai_transient
None LAI [-] False
Open Water LAI [-] False
Perrenial Ice/Snow LAI [-] False
Developed, Medium Intensity LAI [-] True
Developed, High Intensity LAI [-] False
Barren Land LAI [-] False
None LAI [-] False
Open Water LAI [-] False
Perrenial Ice/Snow LAI [-] False
Developed, Medium Intensity LAI [-] True
Developed, High Intensity LAI [-] False
Barren Land LAI [-] False
| time | Developed, Open Space LAI [-] | Developed, Low Intensity LAI [-] | Developed, Medium Intensity LAI [-] | Deciduous Forest LAI [-] | Evergreen Forest LAI [-] | Mixed Forest LAI [-] | Shrub/Scrub LAI [-] | Grassland/Herbaceous LAI [-] | Pasture/Hay LAI [-] | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2010-08-01 00:00:00 | 3.441091 | 3.441091 | 0.0 | 3.441091 | 3.441091 | 3.441091 | 3.441091 | 3.441091 | 3.441091 |
| 1 | 2010-08-05 00:00:00 | 4.188735 | 4.188735 | 0.0 | 4.188735 | 4.188735 | 4.188735 | 4.188735 | 4.188735 | 4.188735 |
| 2 | 2010-08-09 00:00:00 | 4.219150 | 4.219150 | 0.0 | 4.219150 | 4.219150 | 4.219150 | 4.219150 | 4.219150 | 4.219150 |
| 3 | 2010-08-13 00:00:00 | 3.846493 | 3.846493 | 0.0 | 3.846493 | 3.846493 | 3.846493 | 3.846493 | 3.846493 | 3.846493 |
| 4 | 2010-08-17 00:00:00 | 3.173528 | 3.173528 | 0.0 | 3.173528 | 3.173528 | 3.173528 | 3.173528 | 3.173528 | 3.173528 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 84 | 2011-06-30 00:00:00 | 4.754583 | 4.754583 | 0.0 | 4.754583 | 4.754583 | 4.754583 | 4.754583 | 4.754583 | 4.754583 |
| 85 | 2011-07-04 00:00:00 | 3.759242 | 3.759242 | 0.0 | 3.759242 | 3.759242 | 3.759242 | 3.759242 | 3.759242 | 3.759242 |
| 86 | 2011-07-08 00:00:00 | 3.311572 | 3.311572 | 0.0 | 3.311572 | 3.311572 | 3.311572 | 3.311572 | 3.311572 | 3.311572 |
| 87 | 2011-07-12 00:00:00 | 3.859421 | 3.859421 | 0.0 | 3.859421 | 3.859421 | 3.859421 | 3.859421 | 3.859421 | 3.859421 |
| 88 | 2011-07-16 00:00:00 | 5.988454 | 5.988454 | 0.0 | 5.988454 | 5.988454 | 5.988454 | 5.988454 | 5.988454 | 5.988454 |
89 rows × 10 columns
# write the NLCD-based time series to disk
output_filenames['nlcd_lai_cyclic_steadystate'] = toOutput(f'{name}_LAI_NLCD_CyclicSteadystate.h5')
watershed_workflow.io.writeTimeseriesToHDF5(output_filenames['nlcd_lai_cyclic_steadystate'], nlcd_lai_cyclic_steadystate)
output_filenames['nlcd_lai_transient'] = toOutput(f'{name}_LAI_NLCD_{start.year}_{end.year}.h5')
watershed_workflow.io.writeTimeseriesToHDF5(output_filenames['nlcd_lai_transient'], nlcd_lai_transient)
2026-03-09 16:30:54,455 - root - INFO: Writing HDF5 file: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/output_data/Coweeta_LAI_NLCD_CyclicSteadystate.h5
2026-03-09 16:30:54,458 - root - INFO: Writing HDF5 file: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/output_data/Coweeta_LAI_NLCD_2010_2011.h5
Subsurface Soil, Geologic Structure#
NRCS Soils#
# get NRCS shapes, on a reasonable crs
nrcs = sources['soil structure'].getShapesByGeometry(watershed.exterior, watershed.crs, out_crs=crs)
2026-03-09 16:30:54,496 - root - INFO: Shapes manager: buffered+snapped query box = (-83.482, 35.025, -83.419, 35.077)
2026-03-09 16:30:54,497 - root - INFO: Using superset cache: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/input_data/soil_structure/SSURGO/SSURGO_-83.4820_35.0250_-83.4190_35.0770.gpkg
2026-03-09 16:30:54,504 - root - INFO: fixing column: geometry
nrcs
| mukey | residual saturation [-] | Rosetta porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | Rosetta permeability [m^2] | thickness [m] | permeability [m^2] | porosity [-] | bulk density [g/cm^3] | total sand pct [%] | total silt pct [%] | total clay pct [%] | source | geometry | ID | name | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 545800 | 0.177165 | 0.431041 | 0.000139 | 1.470755 | 8.079687e-13 | 2.03 | 3.429028e-15 | 0.307246 | 1.297356 | 66.356250 | 19.518750 | 14.125000 | NRCS | MULTIPOLYGON (((1444403.769 -650469.761, 14443... | 545800 | NRCS-545800 |
| 1 | 545801 | 0.177493 | 0.432741 | 0.000139 | 1.469513 | 8.184952e-13 | 2.03 | 3.247236e-15 | 0.303714 | 1.292308 | 66.400000 | 19.300000 | 14.300000 | NRCS | MULTIPOLYGON (((1442214.588 -649448.055, 14422... | 545801 | NRCS-545801 |
| 2 | 545803 | 0.172412 | 0.400889 | 0.000150 | 1.491087 | 6.477202e-13 | 2.03 | 2.800000e-12 | 0.379163 | 1.400000 | 66.799507 | 21.700493 | 11.500000 | NRCS | MULTIPOLYGON (((1443525.158 -646896.263, 14435... | 545803 | NRCS-545803 |
| 3 | 545805 | 0.177122 | 0.388687 | 0.000083 | 1.468789 | 3.412748e-13 | 2.03 | 2.800000e-12 | 0.384877 | 1.400000 | 46.721675 | 41.778325 | 11.500000 | NRCS | MULTIPOLYGON (((1447216.231 -649752.647, 14472... | 545805 | NRCS-545805 |
| 4 | 545806 | 0.177122 | 0.388687 | 0.000083 | 1.468789 | 3.412748e-13 | 2.03 | 2.800000e-12 | 0.384877 | 1.400000 | 46.721675 | 41.778325 | 11.500000 | NRCS | MULTIPOLYGON (((1447110.457 -649608.807, 14471... | 545806 | NRCS-545806 |
| 5 | 545807 | 0.177122 | 0.388687 | 0.000083 | 1.468789 | 3.412748e-13 | 2.03 | 2.800000e-12 | 0.384877 | 1.400000 | 46.721675 | 41.778325 | 11.500000 | NRCS | MULTIPOLYGON (((1444160.857 -649303.56, 144415... | 545807 | NRCS-545807 |
| 6 | 545811 | 0.185732 | 0.387543 | 0.000162 | 1.466606 | 4.631920e-13 | 2.03 | 1.196738e-12 | 0.330802 | 1.484203 | 68.011736 | 18.128229 | 13.860034 | NRCS | MULTIPOLYGON (((1442781.585 -650867.273, 14427... | 545811 | NRCS-545811 |
| 8 | 545813 | 0.183468 | 0.398767 | 0.000127 | 1.445858 | 4.296896e-13 | 2.03 | 6.219065e-14 | 0.349442 | 1.410667 | 60.007287 | 26.226047 | 13.766667 | NRCS | MULTIPOLYGON (((1443175.399 -650635.458, 14431... | 545813 | NRCS-545813 |
| 9 | 545814 | 0.183709 | 0.398135 | 0.000126 | 1.444985 | 4.224967e-13 | 2.03 | 5.999907e-14 | 0.344322 | 1.412931 | 59.790685 | 26.427142 | 13.782173 | NRCS | MULTIPOLYGON (((1444319.018 -650468.397, 14443... | 545814 | NRCS-545814 |
| 10 | 545815 | 0.178116 | 0.409712 | 0.000161 | 1.496402 | 7.229891e-13 | 2.03 | 4.813863e-14 | 0.314865 | 1.392630 | 70.125671 | 16.487364 | 13.386966 | NRCS | MULTIPOLYGON (((1443295.561 -650613.903, 14432... | 545815 | NRCS-545815 |
| 13 | 545818 | 0.177633 | 0.449923 | 0.000153 | 1.480849 | 1.098837e-12 | 2.03 | 2.912604e-12 | 0.305511 | 1.250023 | 70.847537 | 13.736515 | 15.415948 | NRCS | MULTIPOLYGON (((1442501.014 -650852.369, 14425... | 545818 | NRCS-545818 |
| 14 | 545819 | 0.180064 | 0.449370 | 0.000150 | 1.470521 | 1.020176e-12 | 2.03 | 2.906012e-12 | 0.313860 | 1.255338 | 69.872443 | 14.111475 | 16.016082 | NRCS | MULTIPOLYGON (((1443778.887 -649765.75, 144377... | 545819 | NRCS-545819 |
| 15 | 545820 | 0.177633 | 0.449923 | 0.000153 | 1.480849 | 1.098837e-12 | 2.03 | 2.912604e-12 | 0.305511 | 1.250023 | 70.847537 | 13.736515 | 15.415948 | NRCS | MULTIPOLYGON (((1443488.297 -648842.039, 14434... | 545820 | NRCS-545820 |
| 17 | 545829 | 0.198664 | 0.373986 | 0.000127 | 1.404247 | 2.287242e-13 | 2.03 | 1.867662e-12 | 0.423067 | 1.527566 | 57.045566 | 27.453895 | 15.500539 | NRCS | MULTIPOLYGON (((1444974.354 -650311.087, 14449... | 545829 | NRCS-545829 |
| 18 | 545830 | 0.198828 | 0.374607 | 0.000129 | 1.404487 | 2.325230e-13 | 2.03 | 1.586967e-12 | 0.423149 | 1.526750 | 57.519428 | 26.894939 | 15.585633 | NRCS | MULTIPOLYGON (((1443444.381 -650734.223, 14434... | 545830 | NRCS-545830 |
| 19 | 545831 | 0.200323 | 0.374581 | 0.000129 | 1.401540 | 2.270522e-13 | 2.03 | 1.568151e-12 | 0.421708 | 1.529598 | 57.567135 | 26.542643 | 15.890222 | NRCS | MULTIPOLYGON (((1442467.988 -650266.306, 14424... | 545831 | NRCS-545831 |
| 20 | 545835 | 0.204093 | 0.420703 | 0.000124 | 1.401781 | 3.796359e-13 | 2.03 | 9.308050e-13 | 0.378835 | 1.377864 | 58.736514 | 21.177756 | 20.085730 | NRCS | MULTIPOLYGON (((1443974.275 -648217.413, 14439... | 545835 | NRCS-545835 |
| 21 | 545836 | 0.226801 | 0.383730 | 0.000111 | 1.358205 | 1.456048e-13 | 2.03 | 8.624357e-13 | 0.417356 | 1.526053 | 51.125051 | 26.904748 | 21.970201 | NRCS | MULTIPOLYGON (((1445257.26 -648565.095, 144528... | 545836 | NRCS-545836 |
| 22 | 545837 | 0.216176 | 0.382063 | 0.000117 | 1.374897 | 1.769888e-13 | 2.03 | 1.020980e-12 | 0.371527 | 1.519777 | 53.723579 | 26.663749 | 19.612673 | NRCS | MULTIPOLYGON (((1445383.222 -649594.905, 14454... | 545837 | NRCS-545837 |
| 23 | 545838 | 0.218648 | 0.380960 | 0.000123 | 1.368965 | 1.748097e-13 | 2.03 | 1.032496e-12 | 0.341661 | 1.532889 | 55.303996 | 24.584165 | 20.111839 | NRCS | MULTIPOLYGON (((1445485.734 -648533.032, 14454... | 545838 | NRCS-545838 |
| 24 | 545842 | 0.193278 | 0.412219 | 0.000142 | 1.434838 | 4.783588e-13 | 2.03 | 9.952892e-13 | 0.386946 | 1.400000 | 64.415764 | 18.601478 | 16.982759 | NRCS | MULTIPOLYGON (((1444928.723 -646359.362, 14449... | 545842 | NRCS-545842 |
| 25 | 545843 | 0.200565 | 0.409655 | 0.000118 | 1.409272 | 3.388668e-13 | 2.03 | 9.952892e-13 | 0.387980 | 1.400000 | 56.982759 | 24.706897 | 18.310345 | NRCS | MULTIPOLYGON (((1443481.962 -646966.156, 14434... | 545843 | NRCS-545843 |
| 28 | 545853 | 0.177032 | 0.396018 | 0.000122 | 1.457947 | 4.559355e-13 | 2.03 | 2.078244e-12 | 0.466882 | 1.402853 | 58.821951 | 29.044860 | 12.133189 | NRCS | MULTIPOLYGON (((1443043.461 -650547.797, 14430... | 545853 | NRCS-545853 |
| 29 | 545854 | 0.176966 | 0.396122 | 0.000122 | 1.458096 | 4.569603e-13 | 2.03 | 2.103007e-12 | 0.467110 | 1.402250 | 58.811902 | 29.064959 | 12.123140 | NRCS | MULTIPOLYGON (((1444241.168 -650400.199, 14442... | 545854 | NRCS-545854 |
| 30 | 545855 | 0.149582 | 0.384703 | 0.000247 | 1.928427 | 2.349139e-12 | 2.03 | 6.238368e-12 | 0.235555 | 1.474297 | 84.857230 | 9.099160 | 6.043611 | NRCS | MULTIPOLYGON (((1446223.595 -646598.551, 14462... | 545855 | NRCS-545855 |
| 31 | 545857 | 0.175471 | 0.416598 | 0.000147 | 1.484106 | 7.390878e-13 | 2.03 | 4.981053e-15 | 0.405556 | 1.350000 | 67.400000 | 19.600000 | 13.000000 | NRCS | MULTIPOLYGON (((1445114.169 -650391.022, 14451... | 545857 | NRCS-545857 |
| 32 | 545859 | 0.204379 | 0.380998 | 0.000120 | 1.396307 | 2.167880e-13 | 2.03 | 1.308653e-12 | 0.268713 | 1.505902 | 55.097413 | 27.811815 | 17.090772 | NRCS | MULTIPOLYGON (((1446199.661 -646550.119, 14461... | 545859 | NRCS-545859 |
| 33 | 545860 | 0.210140 | 0.385375 | 0.000111 | 1.388771 | 1.965020e-13 | 2.03 | 1.079711e-12 | 0.288011 | 1.492529 | 52.434716 | 29.038207 | 18.527077 | NRCS | MULTIPOLYGON (((1445904.904 -648371.627, 14459... | 545860 | NRCS-545860 |
| 34 | 545861 | 0.202264 | 0.393882 | 0.000117 | 1.404424 | 2.644035e-13 | 2.03 | 1.968721e-12 | 0.309901 | 1.455172 | 55.476355 | 26.989163 | 17.534483 | NRCS | MULTIPOLYGON (((1445925.111 -646332.769, 14459... | 545861 | NRCS-545861 |
| 36 | 545863 | 0.199673 | 0.376928 | 0.000127 | 1.403759 | 2.333276e-13 | 2.03 | 1.515968e-12 | 0.263536 | 1.518223 | 57.066355 | 27.034187 | 15.899458 | NRCS | MULTIPOLYGON (((1446608.795 -646100.816, 14466... | 545863 | NRCS-545863 |
| 37 | 545874 | 0.207824 | 0.402727 | 0.000134 | 1.396833 | 3.047098e-13 | 2.03 | 1.004831e-12 | 0.366906 | 1.452913 | 60.366159 | 19.989638 | 19.644204 | NRCS | MULTIPOLYGON (((1446310.135 -649449.476, 14463... | 545874 | NRCS-545874 |
| 38 | 545875 | 0.204926 | 0.403617 | 0.000136 | 1.404217 | 3.310600e-13 | 2.03 | 1.116194e-12 | 0.360041 | 1.446720 | 61.381676 | 19.558623 | 19.059702 | NRCS | MULTIPOLYGON (((1445646.413 -648834.921, 14456... | 545875 | NRCS-545875 |
| 39 | 545876 | 0.184037 | 0.452018 | 0.000141 | 1.449088 | 9.100580e-13 | 2.03 | 2.800000e-12 | 0.324877 | 1.247752 | 67.266810 | 15.562931 | 17.170259 | NRCS | MULTIPOLYGON (((1447374.926 -647819.074, 14473... | 545876 | NRCS-545876 |
| 40 | 545878 | 0.196924 | 0.425599 | 0.000125 | 1.416495 | 4.633338e-13 | 1.85 | 2.282599e-12 | 0.364041 | 1.346733 | 59.965519 | 21.381982 | 18.652500 | NRCS | MULTIPOLYGON (((1442931.528 -650568.692, 14429... | 545878 | NRCS-545878 |
| 41 | 545882 | 0.182956 | 0.364911 | 0.000102 | 1.437312 | 2.342078e-13 | 2.03 | 2.894378e-12 | 0.280336 | 1.518895 | 49.369529 | 39.151513 | 11.478958 | NRCS | POLYGON ((1446277.275 -646404.037, 1446291.853... | 545882 | NRCS-545882 |
| 42 | 545885 | 0.165660 | 0.400572 | 0.000183 | 1.574297 | 9.869368e-13 | 2.03 | 2.800000e-12 | 0.326502 | 1.413645 | 74.559606 | 15.282759 | 10.157635 | NRCS | MULTIPOLYGON (((1443029.67 -648603.973, 144299... | 545885 | NRCS-545885 |
| 43 | 545886 | 0.167329 | 0.400668 | 0.000178 | 1.555671 | 9.088520e-13 | 2.03 | 1.531754e-12 | 0.311297 | 1.412906 | 73.402973 | 16.083642 | 10.513386 | NRCS | MULTIPOLYGON (((1441954.621 -648148.049, 14419... | 545886 | NRCS-545886 |
| 44 | 545887 | 0.165406 | 0.401866 | 0.000178 | 1.562874 | 9.604167e-13 | 2.03 | 1.754223e-12 | 0.313298 | 1.404491 | 73.621601 | 16.263744 | 10.114655 | NRCS | POLYGON ((1442719.417 -648589.931, 1442738.136... | 545887 | NRCS-545887 |
# create a clean dataframe with just the data we will need for ATS
def replace_column_nans(df, col_nan, col_replacement):
"""In a df, replace col_nan entries by col_replacement if is nan. In Place!"""
row_indexer = df[col_nan].isna()
df.loc[row_indexer, col_nan] = df.loc[row_indexer, col_replacement]
return
# where poro or perm is nan, put Rosetta poro
replace_column_nans(nrcs, 'porosity [-]', 'Rosetta porosity [-]')
replace_column_nans(nrcs, 'permeability [m^2]', 'Rosetta permeability [m^2]')
# drop unnecessary columns
for col in ['Rosetta porosity [-]', 'Rosetta permeability [m^2]', 'bulk density [g/cm^3]', 'total sand pct [%]',
'total silt pct [%]', 'total clay pct [%]']:
nrcs.pop(col)
# drop nans
nan_mask = nrcs.isna().any(axis=1)
dropped_mukeys = nrcs.index[nan_mask]
# Drop those rows
nrcs = nrcs[~nan_mask]
assert nrcs['porosity [-]'][:].min() >= min_porosity
assert nrcs['permeability [m^2]'][:].max() <= max_permeability
nrcs
# check for nans
nrcs.isna().any()
mukey False
residual saturation [-] False
van Genuchten alpha [Pa^-1] False
van Genuchten n [-] False
thickness [m] False
permeability [m^2] False
porosity [-] False
source False
geometry False
ID False
name False
dtype: bool
# Compute the soil color of each cell of the mesh
# Note, we use mukey here because it is an int, while ID is a string
soil_color_mukey = watershed_workflow.getShapePropertiesOnMesh(m2, nrcs, 'mukey',
resolution=50, nodata=-999)
nrcs.set_index('mukey', drop=False, inplace=True)
unique_soil_colors = list(np.unique(soil_color_mukey))
if -999 in unique_soil_colors:
unique_soil_colors.remove(-999)
# retain only the unique values of soil_color
nrcs = nrcs.loc[unique_soil_colors]
# renumber the ones we know will appear with an ATS ID using ATS conventions
nrcs['ATS ID'] = range(1000, 1000+len(unique_soil_colors))
nrcs.set_index('ATS ID', drop=True, inplace=True)
# create a new soil color and a soil thickness map using the ATS IDs
soil_color = -np.ones_like(soil_color_mukey)
soil_thickness = np.nan * np.ones(soil_color.shape, 'd')
for ats_ID, ID, thickness in zip(nrcs.index, nrcs.mukey, nrcs['thickness [m]']):
mask = np.where(soil_color_mukey == ID)
soil_thickness[mask] = thickness
soil_color[mask] = ats_ID
m2.cell_data['soil_color'] = soil_color
m2.cell_data['soil thickness'] = soil_thickness
# plot the soil color
# -- get a cmap for soil color
sc_indices, sc_cmap, sc_norm, sc_ticks, sc_labels = \
watershed_workflow.colors.createIndexedColormap(nrcs.index)
mp = m2.plot(facecolors=m2.cell_data['soil_color'], cmap=sc_cmap, norm=sc_norm, edgecolors=None, colorbar=False)
watershed_workflow.colors.createIndexedColorbar(ncolors=len(nrcs),
cmap=sc_cmap, labels=sc_labels, ax=plt.gca())
plt.show()
Depth to Bedrock from SoilGrids#
dtb = sources['depth to bedrock'].getDataset(watershed.exterior, watershed.crs)['band_1']
# the SoilGrids dataset is in cm --> convert to meters
dtb.values = dtb.values/100.
2026-03-09 16:30:54,868 - root - INFO: Incoming shape area = 0.0016041250506451986
2026-03-09 16:30:54,868 - root - INFO: ... buffering incoming shape by 3x native resolution = 0.006249999000004891
2026-03-09 16:30:54,868 - root - INFO: ... buffered shape area = 0.0027499628230988967
2026-03-09 16:30:54,868 - root - INFO: ... snapped bounding box = (-83.48540330906533, 35.020827730027406, -83.41456998706528, 35.08124438702745)
# map to the mesh
m2.cell_data['dtb'] = watershed_workflow.getDatasetOnMesh(m2, dtb, method='linear')
gons = m2.plot(facecolors=m2.cell_data['dtb'], cmap='RdBu', edgecolors=None)
plt.show()
GLHYMPs Geology#
glhymps = sources['geologic structure'].getShapesByGeometry(watershed.exterior.buffer(1000), watershed.crs, out_crs=crs)
glhymps = watershed_workflow.soil_properties.mangleGLHYMPSProperties(glhymps,
min_porosity=min_porosity,
max_permeability=max_permeability,
max_vg_alpha=max_vg_alpha)
# intersect with the buffered geometry -- don't keep extras
glhymps = glhymps[glhymps.intersects(watershed.exterior.buffer(10))]
glhymps
2026-03-09 16:30:55,060 - root - INFO: Shapes manager: buffered+snapped query box = (-9294644.860408794, 3640385.2153247977, -9284578.959260395, 3647327.2161167976)
2026-03-09 16:30:55,069 - root - INFO: fixing column: geometry
| ID | name | source | permeability [m^2] | logk_stdev [-] | porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | residual saturation [-] | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 1793338 | GLHYMPS-1793338 | GLHYMPS | 3.019952e-11 | 1.61 | 0.4 | 0.001 | 1.5 | 0.01 | MULTIPOLYGON (((1377015.684 -683020.035, 13769... |
# quality check -- make sure glymps shapes cover the watershed
print(glhymps.union_all().contains(watershed.exterior))
glhymps
True
| ID | name | source | permeability [m^2] | logk_stdev [-] | porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | residual saturation [-] | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 1793338 | GLHYMPS-1793338 | GLHYMPS | 3.019952e-11 | 1.61 | 0.4 | 0.001 | 1.5 | 0.01 | MULTIPOLYGON (((1377015.684 -683020.035, 13769... |
# clean the data
glhymps.pop('logk_stdev [-]')
assert glhymps['porosity [-]'][:].min() >= min_porosity
assert glhymps['permeability [m^2]'][:].max() <= max_permeability
assert glhymps['van Genuchten alpha [Pa^-1]'][:].max() <= max_vg_alpha
glhymps.isna().any()
ID False
name False
source False
permeability [m^2] False
porosity [-] False
van Genuchten alpha [Pa^-1] False
van Genuchten n [-] False
residual saturation [-] False
geometry False
dtype: bool
# note that for larger areas there are often common regions -- two labels with the same properties -- no need to duplicate those with identical values.
def reindex_remove_duplicates(df, index):
"""Removes duplicates, creating a new index and saving the old index as tuples of duplicate values. In place!"""
if index is not None:
if index in df:
df.set_index(index, drop=True, inplace=True)
index_name = df.index.name
# identify duplicate rows
duplicates = list(df.groupby(list(df)).apply(lambda x: tuple(x.index)))
# order is preserved
df.drop_duplicates(inplace=True)
df.reset_index(inplace=True)
df[index_name] = duplicates
return
reindex_remove_duplicates(glhymps, 'ID')
glhymps
| ID | name | source | permeability [m^2] | porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | residual saturation [-] | geometry | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | (1793338,) | GLHYMPS-1793338 | GLHYMPS | 3.019952e-11 | 0.4 | 0.001 | 1.5 | 0.01 | MULTIPOLYGON (((1377015.684 -683020.035, 13769... |
# Compute the geo color of each cell of the mesh
geology_color_glhymps = watershed_workflow.getShapePropertiesOnMesh(m2, glhymps, 'index',
resolution=50, nodata=-999)
# retain only the unique values of geology that actually appear in our cell mesh
unique_geology_colors = list(np.unique(geology_color_glhymps))
if -999 in unique_geology_colors:
unique_geology_colors.remove(-999)
# retain only the unique values of geology_color
glhymps = glhymps.loc[unique_geology_colors]
# renumber the ones we know will appear with an ATS ID using ATS conventions
glhymps['ATS ID'] = range(100, 100+len(unique_geology_colors))
glhymps['TMP_ID'] = glhymps.index
glhymps.reset_index(drop=True, inplace=True)
glhymps.set_index('ATS ID', drop=True, inplace=True)
# create a new geology color using the ATS IDs
geology_color = -np.ones_like(geology_color_glhymps)
for ats_ID, tmp_ID in zip(glhymps.index, glhymps.TMP_ID):
geology_color[np.where(geology_color_glhymps == tmp_ID)] = ats_ID
glhymps.pop('TMP_ID')
m2.cell_data['geology_color'] = geology_color
geology_color_glhymps.min()
<xarray.DataArray 'index' ()> Size: 8B
array(0)
Attributes:
resolution: 50
source_column: index
nodata: -999
crs: +proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_...Combine to form a complete subsurface dataset#
bedrock = watershed_workflow.soil_properties.getDefaultBedrockProperties()
# merge the properties databases
subsurface_props = pd.concat([glhymps, nrcs, bedrock])
# save the properties to disk for use in generating input file
output_filenames['subsurface_properties'] = toOutput(f'{name}_subsurface_properties.csv')
subsurface_props.to_csv(output_filenames['subsurface_properties'])
subsurface_props
| ID | name | source | permeability [m^2] | porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | residual saturation [-] | geometry | mukey | thickness [m] | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 100 | (1793338,) | GLHYMPS-1793338 | GLHYMPS | 3.019952e-11 | 0.400000 | 0.001000 | 1.500000 | 0.010000 | MULTIPOLYGON (((1377015.684 -683020.035, 13769... | NaN | NaN |
| 1000 | 545800 | NRCS-545800 | NRCS | 3.429028e-15 | 0.307246 | 0.000139 | 1.470755 | 0.177165 | MULTIPOLYGON (((1444403.769 -650469.761, 14443... | 545800.0 | 2.03 |
| 1001 | 545801 | NRCS-545801 | NRCS | 3.247236e-15 | 0.303714 | 0.000139 | 1.469513 | 0.177493 | MULTIPOLYGON (((1442214.588 -649448.055, 14422... | 545801.0 | 2.03 |
| 1002 | 545803 | NRCS-545803 | NRCS | 2.800000e-12 | 0.379163 | 0.000150 | 1.491087 | 0.172412 | MULTIPOLYGON (((1443525.158 -646896.263, 14435... | 545803.0 | 2.03 |
| 1003 | 545805 | NRCS-545805 | NRCS | 2.800000e-12 | 0.384877 | 0.000083 | 1.468789 | 0.177122 | MULTIPOLYGON (((1447216.231 -649752.647, 14472... | 545805.0 | 2.03 |
| 1004 | 545806 | NRCS-545806 | NRCS | 2.800000e-12 | 0.384877 | 0.000083 | 1.468789 | 0.177122 | MULTIPOLYGON (((1447110.457 -649608.807, 14471... | 545806.0 | 2.03 |
| 1005 | 545807 | NRCS-545807 | NRCS | 2.800000e-12 | 0.384877 | 0.000083 | 1.468789 | 0.177122 | MULTIPOLYGON (((1444160.857 -649303.56, 144415... | 545807.0 | 2.03 |
| 1006 | 545813 | NRCS-545813 | NRCS | 6.219065e-14 | 0.349442 | 0.000127 | 1.445858 | 0.183468 | MULTIPOLYGON (((1443175.399 -650635.458, 14431... | 545813.0 | 2.03 |
| 1007 | 545814 | NRCS-545814 | NRCS | 5.999907e-14 | 0.344322 | 0.000126 | 1.444985 | 0.183709 | MULTIPOLYGON (((1444319.018 -650468.397, 14443... | 545814.0 | 2.03 |
| 1008 | 545815 | NRCS-545815 | NRCS | 4.813863e-14 | 0.314865 | 0.000161 | 1.496402 | 0.178116 | MULTIPOLYGON (((1443295.561 -650613.903, 14432... | 545815.0 | 2.03 |
| 1009 | 545818 | NRCS-545818 | NRCS | 2.912604e-12 | 0.305511 | 0.000153 | 1.480849 | 0.177633 | MULTIPOLYGON (((1442501.014 -650852.369, 14425... | 545818.0 | 2.03 |
| 1010 | 545819 | NRCS-545819 | NRCS | 2.906012e-12 | 0.313860 | 0.000150 | 1.470521 | 0.180064 | MULTIPOLYGON (((1443778.887 -649765.75, 144377... | 545819.0 | 2.03 |
| 1011 | 545820 | NRCS-545820 | NRCS | 2.912604e-12 | 0.305511 | 0.000153 | 1.480849 | 0.177633 | MULTIPOLYGON (((1443488.297 -648842.039, 14434... | 545820.0 | 2.03 |
| 1012 | 545829 | NRCS-545829 | NRCS | 1.867662e-12 | 0.423067 | 0.000127 | 1.404247 | 0.198664 | MULTIPOLYGON (((1444974.354 -650311.087, 14449... | 545829.0 | 2.03 |
| 1013 | 545830 | NRCS-545830 | NRCS | 1.586967e-12 | 0.423149 | 0.000129 | 1.404487 | 0.198828 | MULTIPOLYGON (((1443444.381 -650734.223, 14434... | 545830.0 | 2.03 |
| 1014 | 545831 | NRCS-545831 | NRCS | 1.568151e-12 | 0.421708 | 0.000129 | 1.401540 | 0.200323 | MULTIPOLYGON (((1442467.988 -650266.306, 14424... | 545831.0 | 2.03 |
| 1015 | 545835 | NRCS-545835 | NRCS | 9.308050e-13 | 0.378835 | 0.000124 | 1.401781 | 0.204093 | MULTIPOLYGON (((1443974.275 -648217.413, 14439... | 545835.0 | 2.03 |
| 1016 | 545836 | NRCS-545836 | NRCS | 8.624357e-13 | 0.417356 | 0.000111 | 1.358205 | 0.226801 | MULTIPOLYGON (((1445257.26 -648565.095, 144528... | 545836.0 | 2.03 |
| 1017 | 545837 | NRCS-545837 | NRCS | 1.020980e-12 | 0.371527 | 0.000117 | 1.374897 | 0.216176 | MULTIPOLYGON (((1445383.222 -649594.905, 14454... | 545837.0 | 2.03 |
| 1018 | 545838 | NRCS-545838 | NRCS | 1.032496e-12 | 0.341661 | 0.000123 | 1.368965 | 0.218648 | MULTIPOLYGON (((1445485.734 -648533.032, 14454... | 545838.0 | 2.03 |
| 1019 | 545842 | NRCS-545842 | NRCS | 9.952892e-13 | 0.386946 | 0.000142 | 1.434838 | 0.193278 | MULTIPOLYGON (((1444928.723 -646359.362, 14449... | 545842.0 | 2.03 |
| 1020 | 545843 | NRCS-545843 | NRCS | 9.952892e-13 | 0.387980 | 0.000118 | 1.409272 | 0.200565 | MULTIPOLYGON (((1443481.962 -646966.156, 14434... | 545843.0 | 2.03 |
| 1021 | 545853 | NRCS-545853 | NRCS | 2.078244e-12 | 0.466882 | 0.000122 | 1.457947 | 0.177032 | MULTIPOLYGON (((1443043.461 -650547.797, 14430... | 545853.0 | 2.03 |
| 1022 | 545854 | NRCS-545854 | NRCS | 2.103007e-12 | 0.467110 | 0.000122 | 1.458096 | 0.176966 | MULTIPOLYGON (((1444241.168 -650400.199, 14442... | 545854.0 | 2.03 |
| 1023 | 545855 | NRCS-545855 | NRCS | 6.238368e-12 | 0.235555 | 0.000247 | 1.928427 | 0.149582 | MULTIPOLYGON (((1446223.595 -646598.551, 14462... | 545855.0 | 2.03 |
| 1024 | 545857 | NRCS-545857 | NRCS | 4.981053e-15 | 0.405556 | 0.000147 | 1.484106 | 0.175471 | MULTIPOLYGON (((1445114.169 -650391.022, 14451... | 545857.0 | 2.03 |
| 1025 | 545859 | NRCS-545859 | NRCS | 1.308653e-12 | 0.268713 | 0.000120 | 1.396307 | 0.204379 | MULTIPOLYGON (((1446199.661 -646550.119, 14461... | 545859.0 | 2.03 |
| 1026 | 545860 | NRCS-545860 | NRCS | 1.079711e-12 | 0.288011 | 0.000111 | 1.388771 | 0.210140 | MULTIPOLYGON (((1445904.904 -648371.627, 14459... | 545860.0 | 2.03 |
| 1027 | 545861 | NRCS-545861 | NRCS | 1.968721e-12 | 0.309901 | 0.000117 | 1.404424 | 0.202264 | MULTIPOLYGON (((1445925.111 -646332.769, 14459... | 545861.0 | 2.03 |
| 1028 | 545874 | NRCS-545874 | NRCS | 1.004831e-12 | 0.366906 | 0.000134 | 1.396833 | 0.207824 | MULTIPOLYGON (((1446310.135 -649449.476, 14463... | 545874.0 | 2.03 |
| 1029 | 545875 | NRCS-545875 | NRCS | 1.116194e-12 | 0.360041 | 0.000136 | 1.404217 | 0.204926 | MULTIPOLYGON (((1445646.413 -648834.921, 14456... | 545875.0 | 2.03 |
| 1030 | 545876 | NRCS-545876 | NRCS | 2.800000e-12 | 0.324877 | 0.000141 | 1.449088 | 0.184037 | MULTIPOLYGON (((1447374.926 -647819.074, 14473... | 545876.0 | 2.03 |
| 1031 | 545878 | NRCS-545878 | NRCS | 2.282599e-12 | 0.364041 | 0.000125 | 1.416495 | 0.196924 | MULTIPOLYGON (((1442931.528 -650568.692, 14429... | 545878.0 | 1.85 |
| 1032 | 545882 | NRCS-545882 | NRCS | 2.894378e-12 | 0.280336 | 0.000102 | 1.437312 | 0.182956 | POLYGON ((1446277.275 -646404.037, 1446291.853... | 545882.0 | 2.03 |
| 1033 | 545885 | NRCS-545885 | NRCS | 2.800000e-12 | 0.326502 | 0.000183 | 1.574297 | 0.165660 | MULTIPOLYGON (((1443029.67 -648603.973, 144299... | 545885.0 | 2.03 |
| 999 | 999 | bedrock | n/a | 1.000000e-16 | 0.050000 | 0.000019 | 1.500000 | 0.010000 | None | NaN | NaN |
Extrude the 2D Mesh to make a 3D mesh#
# set the floor of the domain as max DTB
dtb_max = np.nanmax(m2.cell_data['dtb'].values)
m2.cell_data['dtb'] = m2.cell_data['dtb'].fillna(dtb_max)
print(f'total thickness: {dtb_max} m')
total_thickness = 50.
total thickness: 19.65074688662983 m
# Generate a dz structure for the top 2m of soil
#
# here we try for 10 cells, starting at 5cm at the top and going to 50cm at the bottom of the 2m thick soil
dzs, res = watershed_workflow.mesh.optimizeDzs(0.05, 0.5, 2, 10)
print(dzs)
print(sum(dzs))
[0.05252036 0.06479527 0.09181838 0.14694819 0.25946146 0.39251928
0.49439017 0.49754689]
2.0
# this looks like it would work out, with rounder numbers:
dzs_soil = [0.05, 0.05, 0.05, 0.12, 0.23, 0.5, 0.5, 0.5]
print(sum(dzs_soil))
2.0
# 50m total thickness, minus 2m soil thickness, leaves us with 48 meters to make up.
# optimize again...
dzs2, res2 = watershed_workflow.mesh.optimizeDzs(1, 10, 48, 8)
print(dzs2)
print(sum(dzs2))
# how about...
dzs_geo = [1.0, 2.0, 4.0, 8.0, 11, 11, 11]
print(dzs_geo)
print(sum(dzs_geo))
[ 2.71470626 5.85798457 9.42730917 10. 10. 10. ]
48.0
[1.0, 2.0, 4.0, 8.0, 11, 11, 11]
48.0
# layer extrusion
DTB = m2.cell_data['dtb'].values
soil_color = m2.cell_data['soil_color'].values
geo_color = m2.cell_data['geology_color'].values
soil_thickness = m2.cell_data['soil thickness'].values
# -- data structures needed for extrusion
layer_types = []
layer_data = []
layer_ncells = []
layer_mat_ids = []
# -- soil layer --
depth = 0
for dz in dzs_soil:
depth += 0.5 * dz
layer_types.append('constant')
layer_data.append(dz)
layer_ncells.append(1)
# use glhymps params
br_or_geo = np.where(depth < DTB, geo_color, 999)
soil_or_br_or_geo = np.where(np.bitwise_and(soil_color > 0, depth < soil_thickness),
soil_color,
br_or_geo)
layer_mat_ids.append(soil_or_br_or_geo)
depth += 0.5 * dz
# -- geologic layer --
for dz in dzs_geo:
depth += 0.5 * dz
layer_types.append('constant')
layer_data.append(dz)
layer_ncells.append(1)
geo_or_br = np.where(depth < DTB, geo_color, 999)
layer_mat_ids.append(geo_or_br)
depth += 0.5 * dz
# print the summary
watershed_workflow.mesh.Mesh3D.summarizeExtrusion(layer_types, layer_data,
layer_ncells, layer_mat_ids)
# downselect subsurface properties to only those that are used
layer_mat_id_used = list(np.unique(np.array(layer_mat_ids)))
subsurface_props_used = subsurface_props.loc[layer_mat_id_used]
subsurface_props_used
2026-03-09 16:30:55,576 - root - INFO: Cell summary:
2026-03-09 16:30:55,577 - root - INFO: ------------------------------------------------------------
2026-03-09 16:30:55,577 - root - INFO: l_id | c_id |mat_id | dz | z_top
2026-03-09 16:30:55,577 - root - INFO: ------------------------------------------------------------
2026-03-09 16:30:55,577 - root - INFO: 00 | 00 | 1008 | 0.050000 | 0.000000
2026-03-09 16:30:55,577 - root - INFO: 01 | 01 | 1008 | 0.050000 | 0.050000
2026-03-09 16:30:55,577 - root - INFO: 02 | 02 | 1008 | 0.050000 | 0.100000
2026-03-09 16:30:55,578 - root - INFO: 03 | 03 | 1008 | 0.120000 | 0.150000
2026-03-09 16:30:55,578 - root - INFO: 04 | 04 | 1008 | 0.230000 | 0.270000
2026-03-09 16:30:55,578 - root - INFO: 05 | 05 | 1008 | 0.500000 | 0.500000
2026-03-09 16:30:55,578 - root - INFO: 06 | 06 | 1008 | 0.500000 | 1.000000
2026-03-09 16:30:55,578 - root - INFO: 07 | 07 | 1008 | 0.500000 | 1.500000
2026-03-09 16:30:55,578 - root - INFO: 08 | 08 | 100 | 1.000000 | 2.000000
2026-03-09 16:30:55,578 - root - INFO: 09 | 09 | 100 | 2.000000 | 3.000000
2026-03-09 16:30:55,578 - root - INFO: 10 | 10 | 100 | 4.000000 | 5.000000
2026-03-09 16:30:55,579 - root - INFO: 11 | 11 | 100 | 8.000000 | 9.000000
2026-03-09 16:30:55,579 - root - INFO: 12 | 12 | 999 | 11.000000 | 17.000000
2026-03-09 16:30:55,579 - root - INFO: 13 | 13 | 999 | 11.000000 | 28.000000
2026-03-09 16:30:55,579 - root - INFO: 14 | 14 | 999 | 11.000000 | 39.000000
| ID | name | source | permeability [m^2] | porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | residual saturation [-] | geometry | mukey | thickness [m] | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 100 | (1793338,) | GLHYMPS-1793338 | GLHYMPS | 3.019952e-11 | 0.400000 | 0.001000 | 1.500000 | 0.010000 | MULTIPOLYGON (((1377015.684 -683020.035, 13769... | NaN | NaN |
| 999 | 999 | bedrock | n/a | 1.000000e-16 | 0.050000 | 0.000019 | 1.500000 | 0.010000 | None | NaN | NaN |
| 1000 | 545800 | NRCS-545800 | NRCS | 3.429028e-15 | 0.307246 | 0.000139 | 1.470755 | 0.177165 | MULTIPOLYGON (((1444403.769 -650469.761, 14443... | 545800.0 | 2.03 |
| 1001 | 545801 | NRCS-545801 | NRCS | 3.247236e-15 | 0.303714 | 0.000139 | 1.469513 | 0.177493 | MULTIPOLYGON (((1442214.588 -649448.055, 14422... | 545801.0 | 2.03 |
| 1002 | 545803 | NRCS-545803 | NRCS | 2.800000e-12 | 0.379163 | 0.000150 | 1.491087 | 0.172412 | MULTIPOLYGON (((1443525.158 -646896.263, 14435... | 545803.0 | 2.03 |
| 1003 | 545805 | NRCS-545805 | NRCS | 2.800000e-12 | 0.384877 | 0.000083 | 1.468789 | 0.177122 | MULTIPOLYGON (((1447216.231 -649752.647, 14472... | 545805.0 | 2.03 |
| 1004 | 545806 | NRCS-545806 | NRCS | 2.800000e-12 | 0.384877 | 0.000083 | 1.468789 | 0.177122 | MULTIPOLYGON (((1447110.457 -649608.807, 14471... | 545806.0 | 2.03 |
| 1005 | 545807 | NRCS-545807 | NRCS | 2.800000e-12 | 0.384877 | 0.000083 | 1.468789 | 0.177122 | MULTIPOLYGON (((1444160.857 -649303.56, 144415... | 545807.0 | 2.03 |
| 1006 | 545813 | NRCS-545813 | NRCS | 6.219065e-14 | 0.349442 | 0.000127 | 1.445858 | 0.183468 | MULTIPOLYGON (((1443175.399 -650635.458, 14431... | 545813.0 | 2.03 |
| 1007 | 545814 | NRCS-545814 | NRCS | 5.999907e-14 | 0.344322 | 0.000126 | 1.444985 | 0.183709 | MULTIPOLYGON (((1444319.018 -650468.397, 14443... | 545814.0 | 2.03 |
| 1008 | 545815 | NRCS-545815 | NRCS | 4.813863e-14 | 0.314865 | 0.000161 | 1.496402 | 0.178116 | MULTIPOLYGON (((1443295.561 -650613.903, 14432... | 545815.0 | 2.03 |
| 1009 | 545818 | NRCS-545818 | NRCS | 2.912604e-12 | 0.305511 | 0.000153 | 1.480849 | 0.177633 | MULTIPOLYGON (((1442501.014 -650852.369, 14425... | 545818.0 | 2.03 |
| 1010 | 545819 | NRCS-545819 | NRCS | 2.906012e-12 | 0.313860 | 0.000150 | 1.470521 | 0.180064 | MULTIPOLYGON (((1443778.887 -649765.75, 144377... | 545819.0 | 2.03 |
| 1011 | 545820 | NRCS-545820 | NRCS | 2.912604e-12 | 0.305511 | 0.000153 | 1.480849 | 0.177633 | MULTIPOLYGON (((1443488.297 -648842.039, 14434... | 545820.0 | 2.03 |
| 1012 | 545829 | NRCS-545829 | NRCS | 1.867662e-12 | 0.423067 | 0.000127 | 1.404247 | 0.198664 | MULTIPOLYGON (((1444974.354 -650311.087, 14449... | 545829.0 | 2.03 |
| 1013 | 545830 | NRCS-545830 | NRCS | 1.586967e-12 | 0.423149 | 0.000129 | 1.404487 | 0.198828 | MULTIPOLYGON (((1443444.381 -650734.223, 14434... | 545830.0 | 2.03 |
| 1014 | 545831 | NRCS-545831 | NRCS | 1.568151e-12 | 0.421708 | 0.000129 | 1.401540 | 0.200323 | MULTIPOLYGON (((1442467.988 -650266.306, 14424... | 545831.0 | 2.03 |
| 1015 | 545835 | NRCS-545835 | NRCS | 9.308050e-13 | 0.378835 | 0.000124 | 1.401781 | 0.204093 | MULTIPOLYGON (((1443974.275 -648217.413, 14439... | 545835.0 | 2.03 |
| 1016 | 545836 | NRCS-545836 | NRCS | 8.624357e-13 | 0.417356 | 0.000111 | 1.358205 | 0.226801 | MULTIPOLYGON (((1445257.26 -648565.095, 144528... | 545836.0 | 2.03 |
| 1017 | 545837 | NRCS-545837 | NRCS | 1.020980e-12 | 0.371527 | 0.000117 | 1.374897 | 0.216176 | MULTIPOLYGON (((1445383.222 -649594.905, 14454... | 545837.0 | 2.03 |
| 1018 | 545838 | NRCS-545838 | NRCS | 1.032496e-12 | 0.341661 | 0.000123 | 1.368965 | 0.218648 | MULTIPOLYGON (((1445485.734 -648533.032, 14454... | 545838.0 | 2.03 |
| 1019 | 545842 | NRCS-545842 | NRCS | 9.952892e-13 | 0.386946 | 0.000142 | 1.434838 | 0.193278 | MULTIPOLYGON (((1444928.723 -646359.362, 14449... | 545842.0 | 2.03 |
| 1020 | 545843 | NRCS-545843 | NRCS | 9.952892e-13 | 0.387980 | 0.000118 | 1.409272 | 0.200565 | MULTIPOLYGON (((1443481.962 -646966.156, 14434... | 545843.0 | 2.03 |
| 1021 | 545853 | NRCS-545853 | NRCS | 2.078244e-12 | 0.466882 | 0.000122 | 1.457947 | 0.177032 | MULTIPOLYGON (((1443043.461 -650547.797, 14430... | 545853.0 | 2.03 |
| 1022 | 545854 | NRCS-545854 | NRCS | 2.103007e-12 | 0.467110 | 0.000122 | 1.458096 | 0.176966 | MULTIPOLYGON (((1444241.168 -650400.199, 14442... | 545854.0 | 2.03 |
| 1023 | 545855 | NRCS-545855 | NRCS | 6.238368e-12 | 0.235555 | 0.000247 | 1.928427 | 0.149582 | MULTIPOLYGON (((1446223.595 -646598.551, 14462... | 545855.0 | 2.03 |
| 1024 | 545857 | NRCS-545857 | NRCS | 4.981053e-15 | 0.405556 | 0.000147 | 1.484106 | 0.175471 | MULTIPOLYGON (((1445114.169 -650391.022, 14451... | 545857.0 | 2.03 |
| 1025 | 545859 | NRCS-545859 | NRCS | 1.308653e-12 | 0.268713 | 0.000120 | 1.396307 | 0.204379 | MULTIPOLYGON (((1446199.661 -646550.119, 14461... | 545859.0 | 2.03 |
| 1026 | 545860 | NRCS-545860 | NRCS | 1.079711e-12 | 0.288011 | 0.000111 | 1.388771 | 0.210140 | MULTIPOLYGON (((1445904.904 -648371.627, 14459... | 545860.0 | 2.03 |
| 1027 | 545861 | NRCS-545861 | NRCS | 1.968721e-12 | 0.309901 | 0.000117 | 1.404424 | 0.202264 | MULTIPOLYGON (((1445925.111 -646332.769, 14459... | 545861.0 | 2.03 |
| 1028 | 545874 | NRCS-545874 | NRCS | 1.004831e-12 | 0.366906 | 0.000134 | 1.396833 | 0.207824 | MULTIPOLYGON (((1446310.135 -649449.476, 14463... | 545874.0 | 2.03 |
| 1029 | 545875 | NRCS-545875 | NRCS | 1.116194e-12 | 0.360041 | 0.000136 | 1.404217 | 0.204926 | MULTIPOLYGON (((1445646.413 -648834.921, 14456... | 545875.0 | 2.03 |
| 1030 | 545876 | NRCS-545876 | NRCS | 2.800000e-12 | 0.324877 | 0.000141 | 1.449088 | 0.184037 | MULTIPOLYGON (((1447374.926 -647819.074, 14473... | 545876.0 | 2.03 |
| 1031 | 545878 | NRCS-545878 | NRCS | 2.282599e-12 | 0.364041 | 0.000125 | 1.416495 | 0.196924 | MULTIPOLYGON (((1442931.528 -650568.692, 14429... | 545878.0 | 1.85 |
| 1032 | 545882 | NRCS-545882 | NRCS | 2.894378e-12 | 0.280336 | 0.000102 | 1.437312 | 0.182956 | POLYGON ((1446277.275 -646404.037, 1446291.853... | 545882.0 | 2.03 |
| 1033 | 545885 | NRCS-545885 | NRCS | 2.800000e-12 | 0.326502 | 0.000183 | 1.574297 | 0.165660 | MULTIPOLYGON (((1443029.67 -648603.973, 144299... | 545885.0 | 2.03 |
# extrude
m3 = watershed_workflow.mesh.Mesh3D.extruded_Mesh2D(m2, layer_types, layer_data,
layer_ncells, layer_mat_ids)
print('2D labeled sets')
print('---------------')
for ls in m2.labeled_sets:
print(f'{ls.setid} : {ls.entity} : {len(ls.ent_ids)} : "{ls.name}"')
print('')
print('Extruded 3D labeled sets')
print('------------------------')
for ls in m3.labeled_sets:
print(f'{ls.setid} : {ls.entity} : {len(ls.ent_ids)} : "{ls.name}"')
print('')
print('Extruded 3D side sets')
print('---------------------')
for ls in m3.side_sets:
print(f'{ls.setid} : FACE : {len(ls.cell_list)} : "{ls.name}"')
2D labeled sets
---------------
10000 : CELL : 2611 : "1"
10001 : CELL : 2611 : "1 surface"
10002 : FACE : 76 : "1 boundary"
10003 : FACE : 5 : "1 outlet"
10004 : FACE : 5 : "surface domain outlet"
10005 : CELL : 184 : "river corridor 0 surface"
10006 : CELL : 16 : "stream order 3"
10007 : CELL : 46 : "stream order 2"
10008 : CELL : 122 : "stream order 1"
21 : CELL : 105 : "Developed, Open Space"
22 : CELL : 3 : "Developed, Low Intensity"
23 : CELL : 1 : "Developed, Medium Intensity"
41 : CELL : 1350 : "Deciduous Forest"
42 : CELL : 42 : "Evergreen Forest"
43 : CELL : 1104 : "Mixed Forest"
81 : CELL : 6 : "Pasture/Hay"
Extruded 3D labeled sets
------------------------
10000 : CELL : 39165 : "1"
Extruded 3D side sets
---------------------
1 : FACE : 2611 : "bottom"
2 : FACE : 2611 : "surface"
3 : FACE : 1140 : "external sides"
10001 : FACE : 2611 : "1 surface"
10002 : FACE : 1140 : "1 boundary"
10003 : FACE : 75 : "1 outlet"
10004 : FACE : 75 : "surface domain outlet"
10005 : FACE : 184 : "river corridor 0 surface"
10006 : FACE : 16 : "stream order 3"
10007 : FACE : 46 : "stream order 2"
10008 : FACE : 122 : "stream order 1"
21 : FACE : 105 : "Developed, Open Space"
22 : FACE : 3 : "Developed, Low Intensity"
23 : FACE : 1 : "Developed, Medium Intensity"
41 : FACE : 1350 : "Deciduous Forest"
42 : FACE : 42 : "Evergreen Forest"
43 : FACE : 1104 : "Mixed Forest"
81 : FACE : 6 : "Pasture/Hay"
# save the mesh to disk
output_filenames['mesh'] = toOutput(f'{name}.exo')
try:
os.remove(output_filenames['mesh'])
except FileNotFoundError:
pass
m3.writeExodus(output_filenames['mesh'], 'material id')
2026-03-09 16:30:55,971 - root - INFO: adding side set: 1
2026-03-09 16:30:55,974 - root - INFO: adding side set: 2
2026-03-09 16:30:55,977 - root - INFO: adding side set: 3
2026-03-09 16:30:55,979 - root - INFO: adding side set: 10001
2026-03-09 16:30:55,981 - root - INFO: adding side set: 10002
2026-03-09 16:30:55,982 - root - INFO: adding side set: 10003
2026-03-09 16:30:55,983 - root - INFO: adding side set: 10004
2026-03-09 16:30:55,985 - root - INFO: adding side set: 10005
2026-03-09 16:30:55,986 - root - INFO: adding side set: 10006
2026-03-09 16:30:55,987 - root - INFO: adding side set: 10007
2026-03-09 16:30:55,988 - root - INFO: adding side set: 10008
2026-03-09 16:30:55,989 - root - INFO: adding side set: 21
2026-03-09 16:30:55,990 - root - INFO: adding side set: 22
2026-03-09 16:30:55,991 - root - INFO: adding side set: 23
2026-03-09 16:30:55,991 - root - INFO: adding side set: 41
2026-03-09 16:30:55,993 - root - INFO: adding side set: 42
2026-03-09 16:30:55,994 - root - INFO: adding side set: 43
2026-03-09 16:30:55,996 - root - INFO: adding side set: 81
2026-03-09 16:30:55,997 - root - INFO: adding elem set: 10000
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: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/output_data/Coweeta.exo
Closing exodus file: /Users/Shared/ornldev/code/watershed_workflow/repos/master/examples/Coweeta/output_data/Coweeta.exo