Example: mesh a delineated watershed
Here we mesh the Coweeta Hydrologic Laboratory as an example of how to pull data in from default locations and generate a fully functional ATS mesh.
This might be the worst example to use to learn how to use Watershed Workflows. But it is useful to demonstrate the breadth of problems this project was intended to solve.
This includes a range of datasets:
NHD Plus for river network
NRCS soils data for soil types
NLCD for land cover/transpiration/rooting depths
NED for elevation
[1]:
%matplotlib inline
[2]:
import os,sys
import numpy as np
from matplotlib import pyplot as plt
from matplotlib import cm as pcm
import shapely
import logging
import pandas
pandas.options.display.max_columns = None
import watershed_workflow
import watershed_workflow.source_list
import watershed_workflow.ui
import watershed_workflow.colors
import watershed_workflow.condition
import watershed_workflow.mesh
import watershed_workflow.split_hucs
watershed_workflow.ui.setup_logging(1,None)
figsize = (6,6)
figsize_3d = (8,6)
Sources and setup
Next we set up the source watershed and coordinate system and all data sources for our mesh. We will use the CRS that is included in the shapefile.
[3]:
# specify the input shapefile and a hint as to what HUC it is in.
coweeta_shapefile = 'Coweeta/input_data/coweeta_basin.shp'
hint = '0601' # hint: HUC 4 containing this shape.
# This is necessary to avoid downloading all HUCs to search for this shape
simplify = 30 # length scale to target average edge
logging.info("")
logging.info("Meshing shape: {}".format(coweeta_shapefile))
logging.info("="*30)
# get the shape and crs of the shape
crs, coweeta = watershed_workflow.get_split_form_shapes(coweeta_shapefile)
2021-11-19 12:20:32,538 - root - INFO:
2021-11-19 12:20:32,539 - root - INFO: Meshing shape: Coweeta/input_data/coweeta_basin.shp
2021-11-19 12:20:32,540 - root - INFO: ==============================
2021-11-19 12:20:32,541 - root - INFO:
2021-11-19 12:20:32,542 - root - INFO: Loading shapes
2021-11-19 12:20:32,542 - root - INFO: ------------------------------
2021-11-19 12:20:32,543 - root - INFO: Loading file: 'Coweeta/input_data/coweeta_basin.shp'
2021-11-19 12:20:32,560 - root - INFO: ... found 1 shapes
2021-11-19 12:20:32,561 - root - INFO: Converting to shapely
A wide range of data sources are available; here we use the defaults except for using NHD Plus for watershed boundaries and hydrography (the default is NHD, which is lower resolution and therefore smaller download sizes).
[4]:
# set up a dictionary of source objects
sources = watershed_workflow.source_list.get_default_sources()
watershed_workflow.source_list.log_sources(sources)
2021-11-19 12:20:32,576 - root - INFO: Using sources:
2021-11-19 12:20:32,577 - root - INFO: --------------
2021-11-19 12:20:32,578 - root - INFO: HUC: National Watershed Boundary Dataset (WBD)
2021-11-19 12:20:32,579 - root - INFO: hydrography: National Hydrography Dataset (NHD)
2021-11-19 12:20:32,580 - root - INFO: DEM: National Elevation Dataset (NED)
2021-11-19 12:20:32,581 - root - INFO: soil structure: National Resources Conservation Service Soil Survey (NRCS Soils)
2021-11-19 12:20:32,582 - root - INFO: geologic structure: GLHYMPS version 2.0
2021-11-19 12:20:32,582 - root - INFO: land cover: National Land Cover Database (NLCD) Layer: NLCD_2016_Land_Cover_L48
2021-11-19 12:20:32,583 - root - INFO: soil thickness: None
2021-11-19 12:20:32,584 - root - INFO: meteorology: DayMet 1km
Generate the surface mesh
First we’ll generate the flattened, 2D triangulation, which builds on hydrography data. Then we download a digital elevation map from the National Elevation Dataset, and extrude that 2D triangulation to a 3D surface mesh based on interpolation between pixels of the DEM.
[5]:
# find what HUC our shape is in
huc = watershed_workflow.find_huc(sources['HUC'], coweeta.exterior(), crs, hint, shrink_factor=0.1)
logging.info("Found Coweeta in HUC: {}".format(huc))
2021-11-19 12:20:32,597 - root - INFO:
2021-11-19 12:20:32,598 - root - INFO: Loading HUC 0601
2021-11-19 12:20:32,599 - root - INFO: ------------------------------
2021-11-19 12:20:32,599 - root - INFO:
2021-11-19 12:20:32,600 - root - INFO: Loading level 4 HUCs in 0601
2021-11-19 12:20:32,601 - root - INFO: ------------------------------
2021-11-19 12:20:32,604 - root - INFO: Using HUC file "/Users/uec/code/watershed_workflow/data-library/hydrography/WBD_06_GDB/WBD_06.gdb"
2021-11-19 12:20:32,767 - root - INFO: ... found 1 HUCs
2021-11-19 12:20:32,768 - root - INFO: -- 0601
2021-11-19 12:20:32,769 - root - INFO: Converting to out_crs
2021-11-19 12:20:32,808 - root - INFO: Converting to shapely
2021-11-19 12:20:32,827 - root - INFO: ... found 1
2021-11-19 12:20:32,834 - root - INFO:
2021-11-19 12:20:32,835 - root - INFO: Loading level 6 HUCs in 0601
2021-11-19 12:20:32,835 - root - INFO: ------------------------------
2021-11-19 12:20:32,837 - root - INFO: Using HUC file "/Users/uec/code/watershed_workflow/data-library/hydrography/WBD_06_GDB/WBD_06.gdb"
2021-11-19 12:20:32,885 - root - INFO: ... found 2 HUCs
2021-11-19 12:20:32,885 - root - INFO: -- 060102
2021-11-19 12:20:32,886 - root - INFO: -- 060101
2021-11-19 12:20:32,887 - root - INFO: Converting to out_crs
2021-11-19 12:20:32,952 - root - INFO: Converting to shapely
2021-11-19 12:20:32,986 - root - INFO:
2021-11-19 12:20:32,987 - root - INFO: Loading level 8 HUCs in 060102
2021-11-19 12:20:32,988 - root - INFO: ------------------------------
2021-11-19 12:20:32,990 - root - INFO: Using HUC file "/Users/uec/code/watershed_workflow/data-library/hydrography/WBD_06_GDB/WBD_06.gdb"
2021-11-19 12:20:33,072 - root - INFO: ... found 8 HUCs
2021-11-19 12:20:33,073 - root - INFO: -- 06010201
2021-11-19 12:20:33,074 - root - INFO: -- 06010207
2021-11-19 12:20:33,075 - root - INFO: -- 06010208
2021-11-19 12:20:33,075 - root - INFO: -- 06010205
2021-11-19 12:20:33,076 - root - INFO: -- 06010202
2021-11-19 12:20:33,077 - root - INFO: -- 06010203
2021-11-19 12:20:33,078 - root - INFO: -- 06010204
2021-11-19 12:20:33,078 - root - INFO: -- 06010206
2021-11-19 12:20:33,081 - root - INFO: Converting to out_crs
2021-11-19 12:20:33,174 - root - INFO: Converting to shapely
2021-11-19 12:20:33,207 - root - INFO:
2021-11-19 12:20:33,208 - root - INFO: Loading level 10 HUCs in 06010202
2021-11-19 12:20:33,208 - root - INFO: ------------------------------
2021-11-19 12:20:33,210 - root - INFO: Using HUC file "/Users/uec/code/watershed_workflow/data-library/hydrography/WBD_06_GDB/WBD_06.gdb"
2021-11-19 12:20:33,336 - root - INFO: ... found 5 HUCs
2021-11-19 12:20:33,337 - root - INFO: -- 0601020203
2021-11-19 12:20:33,338 - root - INFO: -- 0601020204
2021-11-19 12:20:33,338 - root - INFO: -- 0601020201
2021-11-19 12:20:33,339 - root - INFO: -- 0601020202
2021-11-19 12:20:33,339 - root - INFO: -- 0601020205
2021-11-19 12:20:33,341 - root - INFO: Converting to out_crs
2021-11-19 12:20:33,363 - root - INFO: Converting to shapely
2021-11-19 12:20:33,371 - root - INFO:
2021-11-19 12:20:33,372 - root - INFO: Loading level 12 HUCs in 0601020201
2021-11-19 12:20:33,372 - root - INFO: ------------------------------
2021-11-19 12:20:33,374 - root - INFO: Using HUC file "/Users/uec/code/watershed_workflow/data-library/hydrography/WBD_06_GDB/WBD_06.gdb"
2021-11-19 12:20:33,587 - root - INFO: ... found 6 HUCs
2021-11-19 12:20:33,588 - root - INFO: -- 060102020106
2021-11-19 12:20:33,588 - root - INFO: -- 060102020104
2021-11-19 12:20:33,589 - root - INFO: -- 060102020105
2021-11-19 12:20:33,589 - root - INFO: -- 060102020101
2021-11-19 12:20:33,589 - root - INFO: -- 060102020102
2021-11-19 12:20:33,590 - root - INFO: -- 060102020103
2021-11-19 12:20:33,592 - root - INFO: Converting to out_crs
2021-11-19 12:20:33,614 - root - INFO: Converting to shapely
2021-11-19 12:20:33,622 - root - INFO: Found Coweeta in HUC: 060102020103
[6]:
rivers = True
if rivers:
# download/collect the river network within that shape's bounds
_, reaches = watershed_workflow.get_reaches(sources['hydrography'], huc,
coweeta.exterior().bounds, crs, crs)
# simplify and prune rivers not IN the shape, constructing a tree-like data structure
# for the river network
rivers = watershed_workflow.simplify_and_prune(coweeta, reaches, filter=True, simplify=simplify, snap=True,
cut_intersections=True)
else:
rivers = list()
watershed_workflow.split_hucs.simplify(coweeta, simplify)
2021-11-19 12:20:33,630 - root - INFO:
2021-11-19 12:20:33,631 - root - INFO: Loading Hydrography
2021-11-19 12:20:33,631 - root - INFO: ------------------------------
2021-11-19 12:20:33,632 - root - INFO: Loading streams in HUC 060102020103
2021-11-19 12:20:33,633 - root - INFO: and/or bounds (273971.0911428096, 3878839.6361173145, 279140.9150949494, 3883953.7853134344)
2021-11-19 12:20:33,636 - root - INFO: Using Hydrography file "/Users/uec/code/watershed_workflow/data-library/hydrography/NHD_H_06010202_GDB/NHD_H_06010202.gdb"
2021-11-19 12:20:33,637 - root - INFO: National Hydrography Dataset (NHD): opening '/Users/uec/code/watershed_workflow/data-library/hydrography/NHD_H_06010202_GDB/NHD_H_06010202.gdb' layer 'NHDFlowline' for streams in '(273971.0911428096, 3878839.6361173145, 279140.9150949494, 3883953.7853134344)'
2021-11-19 12:20:33,855 - root - INFO: Filtering reaches not in-network
2021-11-19 12:20:33,856 - root - INFO: ... found 35 reaches
2021-11-19 12:20:33,856 - root - INFO: Converting to shapely
2021-11-19 12:20:33,869 - root - INFO: Converting to out_crs
2021-11-19 12:20:33,878 - root - INFO:
2021-11-19 12:20:33,878 - root - INFO: Simplifying and pruning
2021-11-19 12:20:33,879 - root - INFO: ------------------------------
2021-11-19 12:20:33,879 - root - INFO: Filtering rivers outside of the HUC space
2021-11-19 12:20:33,880 - root - INFO: ...filtering
2021-11-19 12:20:33,890 - root - INFO: ... filtered from 35 to 21 reaches.
2021-11-19 12:20:33,891 - root - INFO: Generating the river tree
2021-11-19 12:20:33,895 - root - INFO: Simplifying rivers
2021-11-19 12:20:33,903 - root - INFO: Simplifying HUCs
2021-11-19 12:20:33,905 - root - INFO: Snapping river and HUC (nearly) coincident nodes
2021-11-19 12:20:33,908 - root - INFO: snapping polygon segment boundaries to river endpoints
2021-11-19 12:20:33,911 - root - INFO: snapping river endpoints to the polygon
2021-11-19 12:20:33,924 - root - INFO: cutting at crossings
2021-11-19 12:20:33,934 - root - INFO: filtering rivers to HUC
2021-11-19 12:20:33,935 - root - INFO: ...filtering
2021-11-19 12:20:33,940 - root - INFO:
2021-11-19 12:20:33,941 - root - INFO: Simplification Diagnostics
2021-11-19 12:20:33,941 - root - INFO: ------------------------------
2021-11-19 12:20:33,946 - root - INFO: river min seg length: 87.37142879298393
2021-11-19 12:20:33,948 - root - INFO: river median seg length: 160.85900507901667
2021-11-19 12:20:33,950 - root - INFO: HUC min seg length: 42.195003270380795
2021-11-19 12:20:33,951 - root - INFO: HUC median seg length: 55.87002795054712
[7]:
# plot what we have so far -- an image of the HUC and its stream network
fig = plt.figure(figsize=figsize)
ax = watershed_workflow.plot.get_ax(crs, fig)
watershed_workflow.plot.hucs(coweeta, crs, ax=ax, color='k', linewidth=1)
watershed_workflow.plot.rivers(rivers, crs, ax=ax, color='red', linewidth=1)
plt.show()
/Users/Shared/ornldev/code/miniconda3/envs/watershed_workflow_DEV-2021-11-10/lib/python3.9/site-packages/pyproj/crs/crs.py:1256: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
return self._crs.to_proj4(version=version)
[8]:
# form a triangulation on the shape + river network
# triangulation refinement:
# Refine triangles if their area (in m^2) is greater than A(d), where d is the
# distance from the triangle centroid to the nearest stream.
# A(d) is a piecewise linear function -- A = A0 if d <= d0, A = A1 if d >= d1, and
# linearly interpolates between the two endpoints.
d0 = 100; d1 = 500
A0 = 1000; A1 = 5000
#A0 = 500; A1 = 2500
#A0 = 100; A1 = 500
# Refine triangles if they get too acute
min_angle = 32 # degrees
# make 2D mesh
mesh_points2, mesh_tris, areas, dists = watershed_workflow.triangulate(coweeta, rivers,
refine_distance=[d0,A0,d1,A1],
refine_min_angle=min_angle,
diagnostics=True)
#mesh_points2, mesh_tris, d = watershed_workflow.triangulate(coweeta, rivers,
# refine_max_area=100000,
# enforce_delaunay=True,
# diagnostics=True)
2021-11-19 12:20:34,201 - root - INFO:
2021-11-19 12:20:34,202 - root - INFO: Meshing
2021-11-19 12:20:34,203 - root - INFO: ------------------------------
2021-11-19 12:20:34,216 - root - INFO: Triangulating...
2021-11-19 12:20:34,217 - root - INFO: 62 points and 62 facets
2021-11-19 12:20:34,218 - root - INFO: checking graph consistency
2021-11-19 12:20:34,218 - root - INFO: tolerance is set to 1
2021-11-19 12:20:34,220 - root - INFO: building graph data structures
2021-11-19 12:20:34,222 - root - INFO: triangle.build...
2021-11-19 12:20:38,350 - root - INFO: ...built: 8086 mesh points and 15945 triangles
2021-11-19 12:20:38,351 - root - INFO: Plotting triangulation diagnostics
2021-11-19 12:20:39,481 - root - INFO: min area = 289.77447509765625
2021-11-19 12:20:39,482 - root - INFO: max area = 4979.40869140625
[9]:
# get a raster for the elevation map, based on NED
dem_profile, dem = watershed_workflow.get_raster_on_shape(sources['DEM'], coweeta.exterior(), crs)
# elevate the triangle nodes to the dem
mesh_points3 = watershed_workflow.elevate(mesh_points2, crs, dem, dem_profile)
2021-11-19 12:20:39,820 - root - INFO:
2021-11-19 12:20:39,821 - root - INFO: Loading Raster
2021-11-19 12:20:39,822 - root - INFO: ------------------------------
2021-11-19 12:20:39,823 - root - INFO: Collecting raster
2021-11-19 12:20:39,827 - root - INFO: Collecting DEMs to tile bounds: [-83.48845037186388, 35.01734099944037, -83.41165773504302, 35.08381933600275]
2021-11-19 12:20:39,829 - root - INFO: Need:
2021-11-19 12:20:39,830 - root - INFO: /Users/uec/code/watershed_workflow/data-library/dem/USGS_NED_1as_n36_w084.tif
2021-11-19 12:20:39,831 - root - INFO: source files already exist!
2021-11-19 12:20:39,853 - root - INFO: ... got raster of shape: (239, 276)
2021-11-19 12:20:39,854 - root - INFO: ... got raster bounds: (-83.48845037186388, 35.08381933600275, -83.41178370519467, 35.01743044711165)
2021-11-19 12:20:39,856 - root - INFO:
2021-11-19 12:20:39,856 - root - INFO: Elevating Triangulation to DEM
2021-11-19 12:20:39,856 - root - INFO: ------------------------------
Plotting the resulting mesh can be done in a variety of ways, including both 3D plots and mapview. We show both here, but hereafter use mapview plots as they are a bit clearer (if not so flashy)…
[10]:
# plot the resulting surface mesh
fig = plt.figure(figsize=figsize_3d)
ax = watershed_workflow.plot.get_ax('3d', fig, window=[0.0,0.2,1,0.8])
cax = fig.add_axes([0.23,0.18,0.58,0.03])
mp = ax.plot_trisurf(mesh_points3[:,0], mesh_points3[:,1], mesh_points3[:,2],
triangles=mesh_tris, cmap='viridis',
edgecolor=(0,0,0,.2), linewidth=0.5)
cb = fig.colorbar(mp, orientation="horizontal", cax=cax)
t = cax.set_title('elevation [m]')
ax.view_init(55,0)
ax.set_xticklabels(list())
ax.set_yticklabels(list())
/Users/Shared/ornldev/code/watershed_workflow/repos/setup_py/watershed_workflow/plot.py:110: MatplotlibDeprecationWarning: Axes3D(fig) adding itself to the figure is deprecated since 3.4. Pass the keyword argument auto_add_to_figure=False and use fig.add_axes(ax) to suppress this warning. The default value of auto_add_to_figure will change to False in mpl3.5 and True values will no longer work in 3.6. This is consistent with other Axes classes.
ax = Axes3D(fig, rect=window)
[10]:
[Text(3878000.0, 0, ''),
Text(3879000.0, 0, ''),
Text(3880000.0, 0, ''),
Text(3881000.0, 0, ''),
Text(3882000.0, 0, ''),
Text(3883000.0, 0, ''),
Text(3884000.0, 0, ''),
Text(3885000.0, 0, '')]
[11]:
# plot the resulting surface mesh
fig = plt.figure(figsize=figsize)
ax = watershed_workflow.plot.get_ax(crs, fig, window=[0.05,0.1,0.9,0.8])
#ax2 = watershed_workflow.plot.get_ax(crs,fig, window=[0.65,0.05,0.3,0.5])
ax2 = ax.inset_axes([0.65,0.05,0.3,0.5])
cbax = fig.add_axes([0.05,0.05,0.9,0.05])
xlim = (275900., 276400.)
ylim = (3882300., 3883000.)
mp = watershed_workflow.plot.triangulation(mesh_points3, mesh_tris, crs, ax=ax,
color='elevation', edgecolor='white', linewidth=0.2)
cbar = fig.colorbar(mp, orientation="horizontal", cax=cbax)
watershed_workflow.plot.hucs(coweeta, crs, ax=ax, color='k', linewidth=1)
watershed_workflow.plot.rivers(rivers, crs, ax=ax, color='red', linewidth=1)
ax.set_aspect('equal', 'datalim')
mp2 = watershed_workflow.plot.triangulation(mesh_points3, mesh_tris, crs, ax=ax2,
color='elevation', edgecolor='white', linewidth=0.2)
watershed_workflow.plot.hucs(coweeta, crs, ax=ax2, color='k', linewidth=1)
watershed_workflow.plot.rivers(rivers, crs, ax=ax2, color='red', linewidth=1.5)
ax2.set_xlim(xlim)
ax2.set_ylim(ylim)
ax2.set_xticks([])
ax2.set_yticks([])
ax.indicate_inset_zoom(ax2, edgecolor='k')
print(ax.get_xlim())
print(ax.get_ylim())
cbar.ax.set_title('elevation [m]')
(273391.6231597869, 279720.3829351625)
(3878583.928534328, 3884209.4927791064)
[11]:
Text(0.5, 1.0, 'elevation [m]')
[12]:
# construct the 2D mesh
m2 = watershed_workflow.mesh.Mesh2D(mesh_points3.copy(), list(mesh_tris))
[13]:
# hydrologically condition the mesh, removing pits
watershed_workflow.condition.fill_pits(m2)
# plot the change between the two meshes
diff = np.copy(mesh_points3)
diff[:,2] = m2.points[:,2] - mesh_points3[:,2]
print("max diff = ", np.abs(diff[:,2]).max())
fig, ax = watershed_workflow.plot.get_ax(crs, figsize=figsize)
watershed_workflow.plot.triangulation(diff, m2.conn, crs, color='elevation', edgecolors='gray',
linewidth=0.2, ax=ax)
ax.set_title('conditioned dz')
plt.show()
max diff = 6.589355426081738
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.
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.
[14]:
# download the NLCD raster
lc_profile, lc_raster = watershed_workflow.get_raster_on_shape(sources['land cover'],
coweeta.exterior(), crs)
# resample the raster to the triangles
lc = watershed_workflow.values_from_raster(m2.centroids(), crs, lc_raster, lc_profile)
# what land cover types did we get?
logging.info('Found land cover dtypes: {}'.format(lc.dtype))
logging.info('Found land cover types: {}'.format(set(lc)))
2021-11-19 12:20:43,447 - root - INFO:
2021-11-19 12:20:43,448 - root - INFO: Loading Raster
2021-11-19 12:20:43,449 - root - INFO: ------------------------------
2021-11-19 12:20:43,450 - root - INFO: Collecting raster
2021-11-19 12:20:43,549 - root - INFO: CRS: PROJCS["Albers_Conical_Equal_Area",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["meters",1],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
2021-11-19 12:20:43,658 - root - INFO: ... got raster of shape: (180, 173)
2021-11-19 12:20:43,746 - root - INFO: ... got raster bounds: (1129275.0, 1410015.0, 1134465.0, 1404615.0)
2021-11-19 12:20:44,122 - root - INFO: Found land cover dtypes: uint8
2021-11-19 12:20:44,124 - root - INFO: Found land cover types: {0, 41, 42, 43, 81, 52, 21, 22, 23}
[15]:
# plot the NLCD data
# -- get the NLCD colormap which uses official NLCD colors and labels
nlcd_indices, nlcd_cmap, nlcd_norm, nlcd_ticks, nlcd_labels = \
watershed_workflow.colors.generate_nlcd_colormap(lc)
# this is just hacking the label names to make them display a bit neater for a cleaner plot. Likely it
# should get put into the NLCD manager instead of here! See #8
nlcd_labels_fw = []
for label in nlcd_labels:
label_fw = label
if len(label) > 15:
if ' ' in label:
lsplit = label.split()
if len(lsplit) == 2:
label_fw = '\n'.join(lsplit)
elif len(lsplit) == 4:
label_fw = '\n'.join([' '.join(lsplit[0:2]),
' '.join(lsplit[2:])])
elif len(lsplit) == 3:
if len(lsplit[0]) > len(lsplit[-1]):
label_fw = '\n'.join([lsplit[0],
' '.join(lsplit[1:])])
else:
label_fw = '\n'.join([' '.join(lsplit[:-1]),
lsplit[-1]])
nlcd_labels_fw.append(label_fw)
# plot the image
fig = plt.figure(figsize=figsize)
ax = watershed_workflow.plot.get_ax(crs, fig)
polys = watershed_workflow.plot.mesh(m2, crs, ax=ax, color=lc, cmap=nlcd_cmap, norm=nlcd_norm, edgecolor='none',
facecolor='color', linewidth=0.5)
mp = pcm.ScalarMappable(norm=nlcd_norm, cmap=nlcd_cmap)
cb = fig.colorbar(mp)
cb.set_ticks(nlcd_ticks)
cb.set_ticklabels(nlcd_labels_fw)
ax.set_title("land cover index")
[15]:
Text(0.5, 1.0, 'land cover index')
Subsurface properties
Get soil structure from SSURGO. By soil structure, here we calculate, for each formation identified in SSURGO, a soil depth, porosity, permeability, and percent sand/silt/clay (which are then handed off to Rosetta to get a van Genuchten model).
Below this soil we also identify a geologic layer provided by GLHYMPS. This provides information about the deeper subsurface.
SSURGO Soil Properties
[16]:
# download the NRCS soils data as shapes and project it onto the mesh
# -- download the shapes
target_bounds = coweeta.exterior().bounds
logging.info('target bounds: {}'.format(target_bounds))
soil_profile, soil_survey, soil_survey_props = watershed_workflow.get_shapes(sources['soil structure'],
target_bounds, crs, crs, properties=True)
# -- determine the NRCS mukey for each soil unit; this uniquely identifies soil
# properties
soil_ids = np.array([shp.properties['mukey'] for shp in soil_survey], np.int32)
# -- color a raster by the polygons (this makes identifying a triangle's value much
# more efficient)
soil_color_profile, soil_color_raster = \
watershed_workflow.color_raster_from_shapes(target_bounds, 10, soil_survey,
soil_ids, crs)
# -- resample the raster to the triangles
soil_color = watershed_workflow.values_from_raster(m2.centroids(), crs,
soil_color_raster, soil_color_profile)
# -- select only those that appear in our color map
soil_survey_props.set_index('mukey', inplace=True, drop=False)
soil_survey_props = soil_survey_props.loc[np.unique(soil_color), :]
2021-11-19 12:20:48,483 - root - INFO: target bounds: (273971.0911428096, 3878839.6361173145, 279140.9150949494, 3883953.7853134344)
2021-11-19 12:20:48,484 - root - INFO:
2021-11-19 12:20:48,484 - root - INFO: Loading shapes
2021-11-19 12:20:48,485 - root - INFO: ------------------------------
2021-11-19 12:20:48,509 - root - INFO: Attempting to download source for target '/Users/uec/code/watershed_workflow/data-library/soil_structure/SSURGO/SSURGO_-83.4790_35.0269_-83.4208_35.0743.shp'
2021-11-19 12:20:48,528 - root - INFO: Found 490 shapes.
2021-11-19 12:20:48,530 - root - INFO: and crs: epsg:4326
2021-11-19 12:20:48,540 - root - INFO: found 43 unique MUKEYs.
2021-11-19 12:20:50,641 - root - INFO: Running Rosetta for van Genutchen parameters
2021-11-19 12:20:50,844 - root - INFO: ... done
2021-11-19 12:20:50,846 - root - INFO: requested 43 values
2021-11-19 12:20:50,847 - root - INFO: got 43 responses
2021-11-19 12:20:50,857 - root - INFO: ... found 490 shapes
2021-11-19 12:20:50,858 - root - INFO: Converting to shapely
2021-11-19 12:20:50,909 - root - INFO: Converting to requested CRS
2021-11-19 12:20:50,980 - root - INFO: Coloring shapes onto raster:
2021-11-19 12:20:50,981 - root - INFO: target_bounds = (273971.0911428096, 3878839.6361173145, 279140.9150949494, 3883953.7853134344)
2021-11-19 12:20:50,981 - root - INFO: out_bounds = [273966.0, 3878839.0, 279146.0, 3883959.0]
2021-11-19 12:20:50,981 - root - INFO: pixel_size = 10
2021-11-19 12:20:50,982 - root - INFO: width = 518, height = 512
2021-11-19 12:20:50,984 - root - INFO: and 43 independent colors of dtype int32
[17]:
soil_survey_props
[17]:
| residual saturation [-] | Rosetta porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | Rosetta permeability [m^2] | mukey | thickness [cm] | permeability [m^2] | porosity [-] | bulk density [g/cm^3] | total sand pct [%] | total silt pct [%] | total clay pct [%] | source | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| mukey | ||||||||||||||
| 545800 | 0.177165 | 0.431041 | 0.000139 | 1.470755 | 8.079687e-13 | 545800 | 203.0 | 3.429028e-15 | 0.307246 | 1.297356 | 66.356250 | 19.518750 | 14.125000 | NRCS |
| 545801 | 0.177493 | 0.432741 | 0.000139 | 1.469513 | 8.184952e-13 | 545801 | 203.0 | 3.247236e-15 | 0.303714 | 1.292308 | 66.400000 | 19.300000 | 14.300000 | NRCS |
| 545803 | 0.172412 | 0.400889 | 0.000150 | 1.491087 | 6.477202e-13 | 545803 | 203.0 | 2.800000e-12 | 0.379163 | 1.400000 | 66.799507 | 21.700493 | 11.500000 | NRCS |
| 545805 | 0.177122 | 0.388687 | 0.000083 | 1.468789 | 3.412748e-13 | 545805 | 203.0 | 2.800000e-12 | 0.384877 | 1.400000 | 46.721675 | 41.778325 | 11.500000 | NRCS |
| 545806 | 0.177122 | 0.388687 | 0.000083 | 1.468789 | 3.412748e-13 | 545806 | 203.0 | 2.800000e-12 | 0.384877 | 1.400000 | 46.721675 | 41.778325 | 11.500000 | NRCS |
| 545807 | 0.177122 | 0.388687 | 0.000083 | 1.468789 | 3.412748e-13 | 545807 | 203.0 | 2.800000e-12 | 0.384877 | 1.400000 | 46.721675 | 41.778325 | 11.500000 | NRCS |
| 545813 | 0.183468 | 0.398767 | 0.000127 | 1.445858 | 4.296896e-13 | 545813 | 203.0 | 6.219065e-14 | 0.349442 | 1.410667 | 60.007287 | 26.226047 | 13.766667 | NRCS |
| 545814 | 0.183709 | 0.398135 | 0.000126 | 1.444985 | 4.224967e-13 | 545814 | 203.0 | 5.999907e-14 | 0.344322 | 1.412931 | 59.790685 | 26.427142 | 13.782173 | NRCS |
| 545815 | 0.178116 | 0.409712 | 0.000161 | 1.496402 | 7.229891e-13 | 545815 | 203.0 | 4.813863e-14 | 0.314865 | 1.392630 | 70.125671 | 16.487364 | 13.386966 | NRCS |
| 545818 | 0.177633 | 0.449923 | 0.000153 | 1.480849 | 1.098837e-12 | 545818 | 203.0 | 2.912604e-12 | 0.305511 | 1.250023 | 70.847537 | 13.736515 | 15.415948 | NRCS |
| 545819 | 0.180064 | 0.449370 | 0.000150 | 1.470521 | 1.020176e-12 | 545819 | 203.0 | 2.906012e-12 | 0.313860 | 1.255338 | 69.872443 | 14.111475 | 16.016082 | NRCS |
| 545820 | 0.177633 | 0.449923 | 0.000153 | 1.480849 | 1.098837e-12 | 545820 | 203.0 | 2.912604e-12 | 0.305511 | 1.250023 | 70.847537 | 13.736515 | 15.415948 | NRCS |
| 545829 | 0.198664 | 0.373986 | 0.000127 | 1.404247 | 2.287242e-13 | 545829 | 203.0 | 1.867662e-12 | 0.423067 | 1.527566 | 57.045566 | 27.453895 | 15.500539 | NRCS |
| 545830 | 0.198828 | 0.374607 | 0.000129 | 1.404487 | 2.325230e-13 | 545830 | 203.0 | 1.586967e-12 | 0.423149 | 1.526750 | 57.519428 | 26.894939 | 15.585633 | NRCS |
| 545831 | 0.200323 | 0.374581 | 0.000129 | 1.401540 | 2.270522e-13 | 545831 | 203.0 | 1.568151e-12 | 0.421708 | 1.529598 | 57.567135 | 26.542643 | 15.890222 | NRCS |
| 545835 | 0.204093 | 0.420703 | 0.000124 | 1.401781 | 3.796359e-13 | 545835 | 203.0 | 9.308050e-13 | 0.378835 | 1.377864 | 58.736514 | 21.177756 | 20.085730 | NRCS |
| 545836 | 0.226801 | 0.383730 | 0.000111 | 1.358205 | 1.456048e-13 | 545836 | 203.0 | 8.624357e-13 | 0.417356 | 1.526053 | 51.125051 | 26.904748 | 21.970201 | NRCS |
| 545837 | 0.216176 | 0.382063 | 0.000117 | 1.374897 | 1.769888e-13 | 545837 | 203.0 | 1.020980e-12 | 0.371527 | 1.519777 | 53.723579 | 26.663749 | 19.612673 | NRCS |
| 545838 | 0.218648 | 0.380960 | 0.000123 | 1.368965 | 1.748097e-13 | 545838 | 203.0 | 1.032496e-12 | 0.341661 | 1.532889 | 55.303996 | 24.584165 | 20.111839 | NRCS |
| 545842 | 0.193278 | 0.412219 | 0.000142 | 1.434838 | 4.783588e-13 | 545842 | 203.0 | 9.952892e-13 | 0.386946 | 1.400000 | 64.415764 | 18.601478 | 16.982759 | NRCS |
| 545843 | 0.200565 | 0.409655 | 0.000118 | 1.409272 | 3.388668e-13 | 545843 | 203.0 | 9.952892e-13 | 0.387980 | 1.400000 | 56.982759 | 24.706897 | 18.310345 | NRCS |
| 545853 | 0.177032 | 0.396018 | 0.000122 | 1.457947 | 4.559355e-13 | 545853 | 203.0 | 2.078244e-12 | 0.466882 | 1.402853 | 58.821951 | 29.044860 | 12.133189 | NRCS |
| 545854 | 0.176966 | 0.396122 | 0.000122 | 1.458096 | 4.569603e-13 | 545854 | 203.0 | 2.103007e-12 | 0.467110 | 1.402250 | 58.811902 | 29.064959 | 12.123140 | NRCS |
| 545855 | 0.149582 | 0.384703 | 0.000247 | 1.928427 | 2.349139e-12 | 545855 | 203.0 | 6.238368e-12 | 0.235555 | 1.474297 | 84.857230 | 9.099160 | 6.043611 | NRCS |
| 545857 | 0.175471 | 0.416598 | 0.000147 | 1.484106 | 7.390878e-13 | 545857 | 203.0 | 4.981053e-15 | 0.405556 | 1.350000 | 67.400000 | 19.600000 | 13.000000 | NRCS |
| 545859 | 0.204379 | 0.380998 | 0.000120 | 1.396307 | 2.167880e-13 | 545859 | 203.0 | 1.308653e-12 | 0.268713 | 1.505902 | 55.097413 | 27.811815 | 17.090772 | NRCS |
| 545860 | 0.210140 | 0.385375 | 0.000111 | 1.388771 | 1.965020e-13 | 545860 | 203.0 | 1.079711e-12 | 0.288011 | 1.492529 | 52.434716 | 29.038207 | 18.527077 | NRCS |
| 545861 | 0.202264 | 0.393882 | 0.000117 | 1.404424 | 2.644035e-13 | 545861 | 203.0 | 1.968721e-12 | 0.309901 | 1.455172 | 55.476355 | 26.989163 | 17.534483 | NRCS |
| 545863 | 0.199673 | 0.376928 | 0.000127 | 1.403759 | 2.333276e-13 | 545863 | 203.0 | 1.515968e-12 | 0.263536 | 1.518223 | 57.066355 | 27.034187 | 15.899458 | NRCS |
| 545874 | 0.207824 | 0.402727 | 0.000134 | 1.396833 | 3.047098e-13 | 545874 | 203.0 | 1.004831e-12 | 0.366906 | 1.452913 | 60.366159 | 19.989638 | 19.644204 | NRCS |
| 545875 | 0.204926 | 0.403617 | 0.000136 | 1.404217 | 3.310600e-13 | 545875 | 203.0 | 1.116194e-12 | 0.360041 | 1.446720 | 61.381676 | 19.558623 | 19.059702 | NRCS |
| 545876 | 0.184037 | 0.452018 | 0.000141 | 1.449088 | 9.100580e-13 | 545876 | 203.0 | 2.800000e-12 | 0.324877 | 1.247752 | 67.266810 | 15.562931 | 17.170259 | NRCS |
| 545878 | 0.196924 | 0.425599 | 0.000125 | 1.416495 | 4.633338e-13 | 545878 | 185.0 | 2.282599e-12 | 0.364041 | 1.346733 | 59.965519 | 21.381982 | 18.652500 | NRCS |
| 545882 | 0.182956 | 0.364911 | 0.000102 | 1.437312 | 2.342078e-13 | 545882 | 203.0 | 2.894378e-12 | 0.280336 | 1.518895 | 49.369529 | 39.151513 | 11.478958 | NRCS |
| 545885 | 0.165660 | 0.400572 | 0.000183 | 1.574297 | 9.869368e-13 | 545885 | 203.0 | 2.800000e-12 | 0.326502 | 1.413645 | 74.559606 | 15.282759 | 10.157635 | NRCS |
[18]:
# plot the soil mukey
indices, cmap, norm, ticks, labels = watershed_workflow.colors.generate_indexed_colormap(soil_color, cmap='tab20c')
fig, ax = watershed_workflow.plot.get_ax(crs)
mp = watershed_workflow.plot.mesh(m2, crs, ax=ax, facecolor='color',
linewidth=0, color=soil_color,
cmap=cmap, norm = norm
)
watershed_workflow.colors.colorbar_index(ncolors=len(np.unique(soil_color)), cmap=cmap, labels = labels)
ax.set_title('soil type index')
ax.axis('off')
/Users/Shared/ornldev/code/miniconda3/envs/watershed_workflow_DEV-2021-11-10/lib/python3.9/site-packages/pyproj/crs/crs.py:1256: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. See: https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems
return self._crs.to_proj4(version=version)
[18]:
(273712.5998, 279399.40619999997, 3878583.92855, 3884209.4924500003)
[19]:
# Note this is not just the soil ID, but also soil properties.
print(soil_survey_props.keys())
# To demonstrate what we mean by this, plot the porosity of the soil column.
porosity_nrcs = np.empty(soil_color.shape, 'd')
porosity_rosetta = np.empty(soil_color.shape, 'd')
for mukey in soil_survey_props.index:
porosity_nrcs[soil_color == mukey] = soil_survey_props.loc[ mukey,'porosity [-]']
porosity_rosetta[soil_color == mukey] = soil_survey_props.loc[ mukey,'Rosetta porosity [-]']
pmin = min(np.nanmin(porosity_nrcs), np.nanmin(porosity_rosetta))
pmax = max(np.nanmax(porosity_nrcs), np.nanmax(porosity_rosetta))
print('min, max = ', pmin, pmax)
fig = plt.figure(figsize=(12,6))
ax1 = watershed_workflow.plot.get_ax(crs, fig, nrow=1, ncol=2, index=1)
mp = watershed_workflow.plot.triangulation(mesh_points3, mesh_tris, crs, ax=ax1,
color=porosity_nrcs, edgecolor='gray', cmap='jet', vmin=pmin, vmax=pmax)
cbar = fig.colorbar(mp)
ax1.set_title('porosity (NRCS) [-]')
ax1.axis('off')
ax2 = watershed_workflow.plot.get_ax(crs, fig, nrow=1, ncol=2, index=2)
mp = watershed_workflow.plot.triangulation(mesh_points3, mesh_tris, crs, ax=ax2,
color=porosity_rosetta, edgecolor='gray', cmap='jet', vmin=pmin, vmax=pmax)
cbar = fig.colorbar(mp)
ax2.set_title('porosity (Rosetta) [-]')
ax2.axis('off')
Index(['residual saturation [-]', 'Rosetta porosity [-]',
'van Genuchten alpha [Pa^-1]', 'van Genuchten n [-]',
'Rosetta permeability [m^2]', 'mukey', 'thickness [cm]',
'permeability [m^2]', 'porosity [-]', 'bulk density [g/cm^3]',
'total sand pct [%]', 'total silt pct [%]', 'total clay pct [%]',
'source'],
dtype='object')
min, max = 0.2355549116198203 0.4671095937710679
[19]:
(273712.5998, 279399.40619999997, 3878583.92855, 3884209.4924500003)
[20]:
# averaging permeability is a tricky beast. we average in log space, check that unit conversions make sense
fig = plt.figure(figsize=(12,6))
soil_perm_nrcs = np.empty(soil_color.shape, 'd')
soil_perm_rosetta = np.empty(soil_color.shape, 'd')
for mukey in soil_survey_props['mukey']:
soil_perm_nrcs[soil_color == mukey] = soil_survey_props.loc[soil_survey_props['mukey'] == mukey,
'permeability [m^2]']
soil_perm_rosetta[soil_color == mukey] = soil_survey_props.loc[soil_survey_props['mukey'] == mukey,
'Rosetta permeability [m^2]']
pmin = min(np.log10(soil_perm_nrcs).min(), np.log10(soil_perm_rosetta).min())
pmax = max(np.log10(soil_perm_nrcs).max(), np.log10(soil_perm_rosetta).max())
print(f'min = {pmin}, max = {pmax}')
ax1 = watershed_workflow.plot.get_ax(crs, fig, nrow=1, ncol=2, index=1)
mp = watershed_workflow.plot.triangulation(mesh_points3, mesh_tris, crs, ax=ax1,
color=np.log10(soil_perm_nrcs), edgecolor='gray', cmap='jet',
vmin=pmin, vmax=pmax)
cbar = fig.colorbar(mp)
ax1.set_title('log permeability (NRCS) [m^2]')
ax2 = watershed_workflow.plot.get_ax(crs, fig, nrow=1, ncol=2, index=2)
mp = watershed_workflow.plot.triangulation(mesh_points3, mesh_tris, crs, ax=ax2,
color=np.log10(soil_perm_rosetta), edgecolor='gray', cmap='jet',
vmin=pmin, vmax=pmax)
cbar = fig.colorbar(mp)
ax2.set_title('log permeability (Rosetta) [m^2]')
min = -14.488486163789586, max = -11.204929002586569
[20]:
Text(0.5, 1.0, 'log permeability (Rosetta) [m^2]')
GLYHMPS geologic layer
[21]:
# extract the GLYHMPS geologic structure data as shapes and project it onto the mesh
target_bounds = coweeta.exterior().bounds
logging.info('target bounds: {}'.format(target_bounds))
_, geo_survey, geo_survey_props = watershed_workflow.get_shapes(sources['geologic structure'], target_bounds,
crs, crs, properties=True)
# -- log the bounds targetted and found
logging.info('shape union bounds: {}'.format(
shapely.ops.cascaded_union(geo_survey).bounds))
# -- determine the ID for each soil unit; this uniquely identifies formation
# properties
geo_ids = np.array([shp.properties['id'] for shp in geo_survey], np.int32)
# -- color a raster by the polygons (this makes identifying a triangle's value much
# more efficient)
geo_color_profile, geo_color_raster = \
watershed_workflow.color_raster_from_shapes(target_bounds, 10, geo_survey,
geo_ids, crs)
# -- resample the raster to the triangles
geo_color = watershed_workflow.values_from_raster(m2.centroids(), crs,
geo_color_raster, geo_color_profile)
2021-11-19 12:20:58,280 - root - INFO: target bounds: (273971.0911428096, 3878839.6361173145, 279140.9150949494, 3883953.7853134344)
2021-11-19 12:20:58,281 - root - INFO:
2021-11-19 12:20:58,282 - root - INFO: Loading shapes
2021-11-19 12:20:58,283 - root - INFO: ------------------------------
2021-11-19 12:20:58,283 - root - INFO: Getting shapes of GLHYMPS on bounds: (273971.0911428096, 3878839.6361173145, 279140.9150949494, 3883953.7853134344)
2021-11-19 12:20:58,284 - root - INFO: from file: /Users/uec/code/watershed_workflow/data-library/soil_structure/GLHYMPS/GLHYMPS.shp
2021-11-19 12:20:58,425 - fiona.ogrext - INFO: Failed to auto identify EPSG: 7
2021-11-19 12:21:03,478 - root - INFO: ... found 1 shapes
2021-11-19 12:21:03,479 - root - INFO: Converting to shapely
2021-11-19 12:21:03,491 - root - INFO: Converting to requested CRS
2021-11-19 12:21:03,523 - root - INFO: shape union bounds: (159518.27011641115, 3816621.6554112737, 431027.3363569959, 4024643.4346461874)
2021-11-19 12:21:03,524 - root - INFO: Coloring shapes onto raster:
2021-11-19 12:21:03,525 - root - INFO: target_bounds = (273971.0911428096, 3878839.6361173145, 279140.9150949494, 3883953.7853134344)
2021-11-19 12:21:03,527 - root - INFO: out_bounds = [273966.0, 3878839.0, 279146.0, 3883959.0]
2021-11-19 12:21:03,528 - root - INFO: pixel_size = 10
2021-11-19 12:21:03,528 - root - INFO: width = 518, height = 512
2021-11-19 12:21:03,531 - root - INFO: and 1 independent colors of dtype int32
[22]:
# plot the geologic formation id
fig = plt.figure(figsize=figsize)
ax = watershed_workflow.plot.get_ax(crs, fig)
mp = watershed_workflow.plot.mesh(m2, crs, ax=ax, facecolor='color',
linewidth=0, color=geo_color, cmap='tab20c')
ax.set_title('soil type index')
geo_survey_props
[22]:
| id | source | permeability [m^2] | logk_stdev [-] | porosity [-] | van Genuchten alpha [Pa^-1] | van Genuchten n [-] | residual saturation [-] | |
|---|---|---|---|---|---|---|---|---|
| 0 | 1793338 | GLHYMPS | 3.019952e-11 | 1.61 | 0.01 | 0.023953 | 2.0 | 0.01 |
Mesh extrusion
Given the surface mesh and material IDs on both the surface and subsurface, we can extrude the surface mesh in the vertical to make a 3D mesh.
First, all integer IDs in Exodus files must be unique. This includes Material IDs, side sets, etc. We create the Material ID map and data frame. This is used to standardize IDs from multiple data sources. Traditionally, ATS numbers Material IDs/Side Sets as:
0-9 : reserved for boundaries, surface/bottom, etc
10-99 : Land Cover side sets, typically NLCD IDs are used
100-999 : geologic layer material IDs
1000-9999 : soil layer material IDs
[23]:
soil_survey_props['ats_id'] = range(1000, 1000+len(soil_survey_props))
soil_survey_props.set_index('ats_id', inplace=True)
geo_survey_props['ats_id'] = range(100, 100+len(geo_survey_props))
geo_survey_props.set_index('ats_id', inplace=True)
subsurface_props = pandas.concat([geo_survey_props,soil_survey_props])
[24]:
# must choose properties for geologic media. Here we choose one that has a similar porosity
subsurface_props.loc[100, ['residual saturation [-]', 'van Genuchten alpha [Pa^-1]', 'van Genuchten n [-]']] = \
subsurface_props.loc[1024, ['residual saturation [-]', 'van Genuchten alpha [Pa^-1]', 'van Genuchten n [-]']]
# save the properties to disk for use in generating input file
subsurface_props.to_csv(os.path.join('Coweeta', 'output_data', 'coweeta_subsurface_properties.csv'))
Next we extrude the DEM to create a 3D mesh.
The most difficult aspect of extrusion is creating meshes that: 1. aren’t huge numbers of cells 2. aren’t huge cell thicknesses, especially near the surface 3. follow implied interfaces, e.g. bottom of soil and bottom of geologic layer
This is an iterative process that requires some care and some art.
[25]:
# Generate a dz structure for the top 2m of soil -- it appears from above that the soil thickness is uniformly 2m
#
# 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.optimize_dzs(0.05, 0.5, 2, 10)
print(dzs)
[0.0500016 0.05030145 0.06625834 0.1110568 0.22600204 0.49640501
0.49997989 0.49999487]
[26]:
# this looks like it would work out:
dzs = [0.05, 0.05, 0.05, 0.12, 0.23, 0.5, 0.5, 0.5]
print(sum(dzs))
2.0
[27]:
# a 2m soil thickness and a 17m depth to bedrock suggests a geologic layer of 15 - 1m cells
[28]:
# layer extrusion
# -- data structures needed for extrusion
layer_types = []
layer_data = []
layer_ncells = []
layer_mat_ids = []
# -- soil layer --
for dz in dzs:
layer_types.append('constant')
layer_data.append(dz)
layer_ncells.append(1)
layer_mat_ids.append(soil_color)
# -- geologic layer --
layer_types.append('constant')
layer_data.append(15)
layer_ncells.append(15)
layer_mat_ids.append(geo_color)
# print the summary
watershed_workflow.mesh.Mesh3D.summarize_extrusion(layer_types, layer_data,
layer_ncells, layer_mat_ids)
2021-11-19 12:21:08,800 - root - INFO: Cell summary:
2021-11-19 12:21:08,802 - root - INFO: ------------------------------------------------------------
2021-11-19 12:21:08,803 - root - INFO: l_id | c_id |mat_id | dz | z_top
2021-11-19 12:21:08,803 - root - INFO: ------------------------------------------------------------
2021-11-19 12:21:08,804 - root - INFO: 00 | 00 | 545882 | 0.050000 | 0.000000
2021-11-19 12:21:08,805 - root - INFO: 01 | 01 | 545882 | 0.050000 | 0.050000
2021-11-19 12:21:08,805 - root - INFO: 02 | 02 | 545882 | 0.050000 | 0.100000
2021-11-19 12:21:08,806 - root - INFO: 03 | 03 | 545882 | 0.120000 | 0.150000
2021-11-19 12:21:08,806 - root - INFO: 04 | 04 | 545882 | 0.230000 | 0.270000
2021-11-19 12:21:08,807 - root - INFO: 05 | 05 | 545882 | 0.500000 | 0.500000
2021-11-19 12:21:08,807 - root - INFO: 06 | 06 | 545882 | 0.500000 | 1.000000
2021-11-19 12:21:08,808 - root - INFO: 07 | 07 | 545882 | 0.500000 | 1.500000
2021-11-19 12:21:08,809 - root - INFO: 08 | 08 | 1793338 | 1.000000 | 2.000000
2021-11-19 12:21:08,809 - root - INFO: 08 | 09 | 1793338 | 1.000000 | 3.000000
2021-11-19 12:21:08,810 - root - INFO: 08 | 10 | 1793338 | 1.000000 | 4.000000
2021-11-19 12:21:08,810 - root - INFO: 08 | 11 | 1793338 | 1.000000 | 5.000000
2021-11-19 12:21:08,811 - root - INFO: 08 | 12 | 1793338 | 1.000000 | 6.000000
2021-11-19 12:21:08,812 - root - INFO: 08 | 13 | 1793338 | 1.000000 | 7.000000
2021-11-19 12:21:08,812 - root - INFO: 08 | 14 | 1793338 | 1.000000 | 8.000000
2021-11-19 12:21:08,814 - root - INFO: 08 | 15 | 1793338 | 1.000000 | 9.000000
2021-11-19 12:21:08,814 - root - INFO: 08 | 16 | 1793338 | 1.000000 | 10.000000
2021-11-19 12:21:08,816 - root - INFO: 08 | 17 | 1793338 | 1.000000 | 11.000000
2021-11-19 12:21:08,816 - root - INFO: 08 | 18 | 1793338 | 1.000000 | 12.000000
2021-11-19 12:21:08,817 - root - INFO: 08 | 19 | 1793338 | 1.000000 | 13.000000
2021-11-19 12:21:08,818 - root - INFO: 08 | 20 | 1793338 | 1.000000 | 14.000000
2021-11-19 12:21:08,819 - root - INFO: 08 | 21 | 1793338 | 1.000000 | 15.000000
2021-11-19 12:21:08,819 - root - INFO: 08 | 22 | 1793338 | 1.000000 | 16.000000
[29]:
# extrude
m3 = watershed_workflow.mesh.Mesh3D.extruded_Mesh2D(m2, layer_types, layer_data,
layer_ncells, layer_mat_ids)
[30]:
# add back on land cover side sets
surf_ss = m3.side_sets[1]
for index, name in zip(nlcd_indices, nlcd_labels):
where = np.where(lc == index)[0]
ss = watershed_workflow.mesh.SideSet(name, int(index),
[surf_ss.elem_list[w] for w in where],
[surf_ss.side_list[w] for w in where])
m3.side_sets.append(ss)
[31]:
# save to disk
try:
os.remove(os.path.join('Coweeta', 'output_data', 'coweeta_basin.exo'))
except FileNotFoundError:
pass
m3.write_exodus(os.path.join('Coweeta', 'output_data', 'coweeta_basin.exo'))
You are using exodus.py v 1.20.10 (seacas-py3), a python wrapper of some of the exodus library.
Copyright (c) 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021 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: Coweeta/output_data/coweeta_basin.exo
2021-11-19 12:21:19,092 - root - INFO: adding side set: 1
2021-11-19 12:21:19,144 - root - INFO: adding side set: 2
2021-11-19 12:21:19,200 - root - INFO: adding side set: 3
2021-11-19 12:21:19,240 - root - INFO: adding side set: 0
2021-11-19 12:21:19,275 - root - INFO: adding side set: 21
2021-11-19 12:21:19,311 - root - INFO: adding side set: 22
2021-11-19 12:21:19,347 - root - INFO: adding side set: 23
2021-11-19 12:21:19,381 - root - INFO: adding side set: 41
2021-11-19 12:21:19,428 - root - INFO: adding side set: 42
2021-11-19 12:21:19,463 - root - INFO: adding side set: 43
2021-11-19 12:21:19,507 - root - INFO: adding side set: 52
2021-11-19 12:21:19,542 - root - INFO: adding side set: 81
Closing exodus file: Coweeta/output_data/coweeta_basin.exo
[ ]: