20260108
21244500
1.0
FALSE
WettedArea_11000cfs
002
6248347.470000
6341755.310000
2118756.960000
2166597.630000
1
Projected
GCS_WGS_1984
Linear Unit: Meter (1.000000)
WGS_1984_Web_Mercator_Auxiliary_Sphere
<ProjectedCoordinateSystem xsi:type='typens:ProjectedCoordinateSystem' xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' xmlns:xs='http://www.w3.org/2001/XMLSchema' xmlns:typens='http://www.esri.com/schemas/ArcGIS/3.3.0'><WKT>PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0],AUTHORITY["EPSG",3857]]</WKT><XOrigin>-20037700</XOrigin><YOrigin>-30241100</YOrigin><XYScale>10000</XYScale><ZOrigin>-100000</ZOrigin><ZScale>10000</ZScale><MOrigin>-100000</MOrigin><MScale>10000</MScale><XYTolerance>0.001</XYTolerance><ZTolerance>0.001</ZTolerance><MTolerance>0.001</MTolerance><HighPrecision>true</HighPrecision><WKID>102100</WKID><LatestWKID>3857</LatestWKID></ProjectedCoordinateSystem>
20250501
17371000
20250501
17371000
150000000
5000
ISO19139
WettedArea_11000cfs
Feature Class
0
OBJECTID
OID
OID
4
0
0
Internal feature number.
Esri
Sequential unique whole numbers that are automatically generated.
Shape
Shape
Geometry
0
0
0
Feature geometry.
Esri
Coordinates defining the features.
InPoly_FID
InPoly_FID
Integer
4
0
0
Shape_Length
Shape_Length
Double
8
0
0
Length of feature in internal units.
Esri
Positive real numbers that are automatically generated.
Shape_Area
Shape_Area
Double
8
0
0
Area of feature in internal units squared.
Esri
Positive real numbers that are automatically generated.
dataset
Eric Peterson
Trinity River Restoration Program, U.S. Bureau of Reclamation
Science Coordinator
530-623-1800
Weaverville
CA
96093
US
ebpeterson@usbr.gov
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
20250501
ArcGIS Metadata
1.0
File Geodatabase Feature Class
ArcGIS 10.5
Eric Peterson
Trinity River Restoration Program, U.S. Bureau of Reclamation
Data Steward
530-623-1800
Weaverville
CA
96093
US
ebpeterson@usbr.gov
ZIP file
http://www.trrp.net/library/data?id=104
Wetted area model for 11000 cfs on the 2022 topography model
2018-02-28
2018-03-01
Eric Peterson
Trinity River Restoration Program, U.S. Bureau of Reclamation
Science Coordinator
530-623-1800
Weaverville
CA
96093
US
ebpeterson@usbr.gov
TRRP (Trinity River Restoration Program). 2023. Wetted area models on the 2022 topography model. Data Package. TRRP, Weaverville, California.
<DIV STYLE="text-align:Left;"><DIV><DIV><P><SPAN>The SRH2D hydrologic model software was used in conjunction with the 2022 topography model of the Trinity River between Lewiston Dam and the North Fork Trinity River to estimate flows and fish habitat at a variety of flow rates. The model operates within a two-dimensional mesh of variable cell sizes. Wetted area output directly from the model therefore follows model mesh cells, yielding a blocky configuration to the wetted area extent. For the data presented here, the water surface elevation estimates for mesh cells were used to develop a water surface (ArcGIS Terrain format), which was intersected with the 2022 topography model (also ArcGIS Terrain format) for a more accurate representation of the actual aerial extent. Preparation of mass point files for inclusion in terrains included removal of a few modeling artifact polygon centroids where major tributaries approach the river and at the tight river bend in Junction City. Post intersection processing included the removal of polygons covering less than 5 ft2, followed by smoothing of the angular results that come from Terrain intersection, as well as manual removal of additional model artifacts. See the associated document " Visual Explanation of Manual Wetted Area Polygon Edits.pdf" or the lineage section of these metadata for details.</SPAN></P></DIV></DIV></DIV>
Visualization and approximation of the wetted extent of the Trinity River at 11000 cfs given the 2022 topography model.
Processing to wetted area: Jeanne McSloy, Trinity River Restoration Program, U.S. Bureau of Reclamation.
Underlying hydrologic modeling: Dr. Nathan Bradley, U.S. Bureau of Reclamation, Technical Services center.
Eric Peterson
Trinity River Restoration Program, U.S. Bureau of Reclamation
Science Coordinator
530-623-1800
Weaverville
CA
96093
US
ebpeterson@usbr.gov
Trinity River
Trinity County
California
Bank Lines
Inundation Extent
Edge of Water
Trinity River
Trinity County
California
Bank Lines
Inundation Extent
Edge of Water
<DIV STYLE="text-align:Left;"><DIV><DIV><P><SPAN>Public domain; data provided without guarantee or warrantee.</SPAN></P></DIV></DIV></DIV>
Esri ArcGIS Pro 3.3 Patch 2
true
-123.13
-122.79
40.64
40.78
true
-123.1268249
-122.7935261
40.6424812
40.7703269
1
-123.131340
-122.792521
40.774983
40.640852
Eric Peterson
Trinity River Restoration Program, U.S. Bureau of Reclamation
Data Steward
530-623-1800
Weaverville
CA
96093
US
ebpeterson@usbr.gov
Public domain; data provided without guarantee or warrantee.
###### DATA PREPARATION FOR TERRAIN INPUT ######
# PREPARATION: All GDBs to be used be were placed in the same folder
# A GDB was created for output, containing a Feature Dataset with
# x, y and z resolution= 0.01 US ft and tolerance = 0.02 US ft
# consistent with settings used in previous terrain creation
###### STEP 1 - Select wetted mesh cell polygons (WaterDepth>0) and export to a new feature class ######
import arcpy
import os
import sys
import traceback
def setup_environment(workspace_path):
"""
Sets up the environment for processing.
Parameters:
workspace_path (str): Path to the workspace geodatabase.
"""
arcpy.CheckOutExtension("3D") # Check out 3D Analyst extension
arcpy.env.overwriteOutput = True
arcpy.env.workspace = workspace_path
print(f"Environment set. Workspace: {workspace_path}")
def select_and_export_polygons(input_gdb, output_gdb, field_name, value_threshold, suffix):
"""
Selects polygons where a field value exceeds a threshold and exports them to new polygon feature classes.
Parameters:
input_gdb (str): Path to the input geodatabase containing polygon feature classes.
output_gdb (str): Path to the output geodatabase for storing new polygon feature classes.
field_name (str): The field to query (e.g., "WaterDepth").
value_threshold (float): The threshold value for the query.
suffix (str): Suffix to append to the exported feature class name.
"""
arcpy.env.workspace = input_gdb
feature_classes = arcpy.ListFeatureClasses(feature_type="POLYGON")
if not feature_classes:
print("No polygon feature classes found in the input geodatabase.")
return
print(f"Found {len(feature_classes)} polygon feature classes to process.")
for fc in feature_classes:
try:
# Construct the output feature class name and path
output_fc_name = f"{os.path.basename(fc)}_{suffix}"
output_fc_path = os.path.join(output_gdb, output_fc_name)
# Construct the query to filter the data
where_clause = f'"{field_name}" > {value_threshold}'
# Select the polygons that meet the query criteria
selected_layer = arcpy.management.MakeFeatureLayer(fc, "selected_layer", where_clause)
# Skip if no records meet the criteria
if int(arcpy.management.GetCount(selected_layer)[0]) == 0:
print(f"No records found in {fc} with {field_name} > {value_threshold}. Skipping...")
arcpy.management.Delete(selected_layer)
continue
# Export the selected polygons to a new feature class
arcpy.management.CopyFeatures(selected_layer, output_fc_path)
print(f"Exported selected polygons from {fc} to {output_fc_path}")
# Clean up the temporary layer
arcpy.management.Delete(selected_layer)
except Exception as e:
print(f"Error processing {fc}: {e}")
if __name__ == "__main__":
try:
# Define input and output geodatabases
input_gdb = r"Hydraulics"
output_gdb = r"WettedAreaModels_2022.gdb"
print("Starting polygon selection and export script...")
# Set up the environment
setup_environment(input_gdb)
# Step 1: Select polygons and export them
print("\n--- Step 1: Select and Export Polygons ---")
select_and_export_polygons(
input_gdb=input_gdb,
output_gdb=output_gdb,
field_name="WaterDepth",
value_threshold=0,
suffix="selected"
)
print("\n--> Finished script successfully.")
except Exception as e:
print(f"Error running the script: {e}")
traceback.print_exc()
###### STEP 2: Export wetted polygons using field mapping to retain only relevant fields ######
import arcpy
import os
def export_features_with_field_mapping(input_gdb, field_to_retain, input_suffix, output_suffix):
"""
Exports feature classes with a specific suffix, retaining only one specified field,
and appends a new suffix to the output feature classes.
Parameters:
input_gdb (str): Path to the geodatabase containing the feature classes.
field_to_retain (str): Name of the field to retain in the output.
input_suffix (str): Suffix of the input feature classes to process (e.g., "_selected").
output_suffix (str): Suffix to append to the output feature classes (e.g., "_WSE").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Construct the output feature class name
output_fc_name = fc.replace(input_suffix, output_suffix)
output_fc_path = os.path.join(input_gdb, output_fc_name)
# Create the field mapping parameter as a string
field_mapping = f'{field_to_retain} "{field_to_retain}" true true false 8 Double 0 0,First,#,{fc},{field_to_retain},-1,-1'
# Use Export Features tool
arcpy.conversion.ExportFeatures(
in_features=fc,
out_features=output_fc_path,
field_mapping=field_mapping
)
print(f"Exported {fc} to {output_fc_name}, retaining only the field: {field_to_retain}")
except Exception as e:
print(f"Error processing {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase (output from the first script)
input_gdb = r"WettedAreaModels_2022.gdb"
field_to_retain = "WSE" # Field to retain in the output
input_suffix = "_selected" # Suffix of input feature classes
output_suffix = "_WSE" # Suffix for the output feature classes
print("Starting Export Features script...")
# Run the export process
export_features_with_field_mapping(input_gdb, field_to_retain, input_suffix, output_suffix)
print("\n--> Finished script successfully.")
###### STEP 2a: Delete Intermediate Data from GDB ######
import arcpy
import os
def delete_selected_feature_classes(input_gdb, input_suffix):
"""
Deletes feature classes in the geodatabase that end with a specific suffix.
Parameters:
input_gdb (str): Path to the geodatabase containing feature classes to delete.
input_suffix (str): Suffix of the feature classes to delete (e.g., "_selected").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Delete the feature class
arcpy.management.Delete(fc)
print(f"Deleted feature class: {fc}")
except Exception as e:
print(f"Error deleting {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase
input_gdb = r"WettedAreaModels_2022.gdb"
input_suffix = "_selected" # Suffix of feature classes to delete
print("Starting deletion of intermediate files...")
# Run the deletion process
delete_selected_feature_classes(input_gdb, input_suffix)
print("\n--> Finished script successfully.")
###### STEP 3: Export wetted mesh cell centroids as point feature classes in a Feature Dataset (FDS) ######
import arcpy
import os
def feature_to_point(input_gdb, feature_dataset, input_suffix, output_suffix):
"""
Converts polygon feature classes with a specific suffix to point feature classes,
saving the output to a specified feature dataset.
Parameters:
input_gdb (str): Path to the geodatabase containing polygon feature classes.
feature_dataset (str): Name of the feature dataset for saving the output point feature classes.
input_suffix (str): Suffix of the input polygon feature classes to process (e.g., "_WSE").
output_suffix (str): Suffix for the output point feature classes (e.g., "_wetpts").
"""
print("=== Starting Debugging Process ===")
# Validate input geodatabase
if arcpy.Exists(input_gdb):
print(f"Geodatabase exists: {input_gdb}")
else:
print(f"Geodatabase does NOT exist: {input_gdb}")
return
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Validate feature dataset
feature_dataset_path = os.path.join(input_gdb, feature_dataset)
if arcpy.Exists(feature_dataset_path):
print(f"Feature dataset exists: {feature_dataset_path}")
else:
print(f"Feature dataset does NOT exist: {feature_dataset_path}")
return
# List all polygon feature classes in the geodatabase (top-level only)
feature_classes = arcpy.ListFeatureClasses(feature_type="Polygon")
if not feature_classes:
print("No polygon feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} polygon feature classes in the geodatabase: {feature_classes}")
# Filter for feature classes matching the input suffix
filtered_classes = [fc for fc in feature_classes if fc.endswith(input_suffix)]
if not filtered_classes:
print(f"No feature classes found with the suffix '{input_suffix}'.")
return
print(f"Feature classes matching suffix '{input_suffix}': {filtered_classes}")
# Process each matching feature class
for fc in filtered_classes:
try:
# Construct output feature class name and path
output_fc_name = fc.replace(input_suffix, output_suffix)
output_fc_path = os.path.join(feature_dataset_path, output_fc_name)
# Debugging output
print(f"Processing input feature class: {fc}")
print(f"Output path: {output_fc_path}")
# Test Feature to Point operation
arcpy.management.FeatureToPoint(fc, output_fc_path, "INSIDE")
print(f"Converted {fc} to points: {output_fc_name}")
except arcpy.ExecuteError:
print(f"ArcPy error for {fc}: {arcpy.GetMessages(2)}")
except Exception as e:
print(f"General error for {fc}: {e}")
print("\n=== Debugging Process Complete ===")
if __name__ == "__main__":
# Define paths and parameters
input_gdb = r"WettedAreaModels_2022.gdb"
feature_dataset = "AllFLows_FDS" # Feature dataset within the geodatabase
input_suffix = "_WSE" # Input polygon feature class suffix
output_suffix = "_wetpts" # Output point feature class suffix
print("Starting Feature to Point script...")
# Run the Feature to Point process
feature_to_point(input_gdb, feature_dataset, input_suffix, output_suffix)
print("\n--> Finished processing.")
###### STEP 3a Delete Intermediate Data from GDB ######
import arcpy
import os
def feature_to_point(input_gdb, feature_dataset, input_suffix, output_suffix):
"""
Converts polygon feature classes with a specific suffix to point feature classes,
saving the output to a specified feature dataset.
Parameters:
input_gdb (str): Path to the geodatabase containing polygon feature classes.
feature_dataset (str): Name of the feature dataset for saving the output point feature classes.
input_suffix (str): Suffix of the input polygon feature classes to process (e.g., "_WSE").
output_suffix (str): Suffix for the output point feature classes (e.g., "_wetpts").
"""
print("=== Starting Debugging Process ===")
# Validate input geodatabase
if arcpy.Exists(input_gdb):
print(f"Geodatabase exists: {input_gdb}")
else:
print(f"Geodatabase does NOT exist: {input_gdb}")
return
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Validate feature dataset
feature_dataset_path = os.path.join(input_gdb, feature_dataset)
if arcpy.Exists(feature_dataset_path):
print(f"Feature dataset exists: {feature_dataset_path}")
else:
print(f"Feature dataset does NOT exist: {feature_dataset_path}")
return
# List all polygon feature classes in the geodatabase (top-level only)
feature_classes = arcpy.ListFeatureClasses(feature_type="Polygon")
if not feature_classes:
print("No polygon feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} polygon feature classes in the geodatabase: {feature_classes}")
# Filter for feature classes matching the input suffix
filtered_classes = [fc for fc in feature_classes if fc.endswith(input_suffix)]
if not filtered_classes:
print(f"No feature classes found with the suffix '{input_suffix}'.")
return
print(f"Feature classes matching suffix '{input_suffix}': {filtered_classes}")
# Process each matching feature class
for fc in filtered_classes:
try:
# Construct output feature class name and path
output_fc_name = fc.replace(input_suffix, output_suffix)
output_fc_path = os.path.join(feature_dataset_path, output_fc_name)
# Debugging output
print(f"Processing input feature class: {fc}")
print(f"Output path: {output_fc_path}")
# Test Feature to Point operation
arcpy.management.FeatureToPoint(fc, output_fc_path, "INSIDE")
print(f"Converted {fc} to points: {output_fc_name}")
except arcpy.ExecuteError:
print(f"ArcPy error for {fc}: {arcpy.GetMessages(2)}")
except Exception as e:
print(f"General error for {fc}: {e}")
print("\n=== Debugging Process Complete ===")
if __name__ == "__main__":
# Define paths and parameters
input_gdb = r"WettedAreaModels_2022.gdb"
feature_dataset = "AllFLows_FDS" # Feature dataset within the geodatabase
input_suffix = "_WSE" # Input polygon feature class suffix
output_suffix = "_wetpts" # Output point feature class suffix
print("Starting Feature to Point script...")
# Run the Feature to Point process
feature_to_point(input_gdb, feature_dataset, input_suffix, output_suffix)
print("\n--> Finished processing.")
###### STEP 4: Convert points to 3D and assign WSE value as Z for all point feature classes in FDS using Feature to 3D By Attribute ######
import arcpy
import os
def feature_to_3d(input_gdb, feature_dataset, input_suffix, output_suffix, height_field):
"""
Converts point feature classes to 3D by attribute without deleting input feature classes.
Parameters:
input_gdb (str): Path to the geodatabase containing the feature dataset.
feature_dataset (str): Name of the feature dataset containing input point feature classes.
input_suffix (str): Suffix of the input point feature classes to process (e.g., "_wetpts").
output_suffix (str): Suffix for the output 3D point feature classes (e.g., "_WSEz").
height_field (str): Name of the attribute field used for elevation in the 3D conversion.
"""
arcpy.env.workspace = os.path.join(input_gdb, feature_dataset)
arcpy.env.overwriteOutput = True
# Validate feature dataset
feature_dataset_path = os.path.join(input_gdb, feature_dataset)
if not arcpy.Exists(feature_dataset_path):
print(f"Feature dataset does NOT exist: {feature_dataset_path}")
return
print(f"Processing feature dataset: {feature_dataset_path}")
# List all point feature classes in the feature dataset
feature_classes = arcpy.ListFeatureClasses(feature_type="Point")
if not feature_classes:
print("No point feature classes found in the feature dataset.")
return
print(f"Found {len(feature_classes)} point feature classes in the feature dataset.")
# Filter for feature classes matching the input suffix
filtered_classes = [fc for fc in feature_classes if fc.endswith(input_suffix)]
if not filtered_classes:
print(f"No feature classes found with the suffix '{input_suffix}'.")
return
print(f"Feature classes matching suffix '{input_suffix}': {filtered_classes}")
# Process each matching feature class
for fc in filtered_classes:
try:
# Construct the output feature class name and path
output_fc_name = fc.replace(input_suffix, output_suffix)
output_fc_path = os.path.join(feature_dataset_path, output_fc_name)
# Debugging output
print(f"Processing input feature class: {fc}")
print(f"Saving output 3D feature class: {output_fc_path}")
# Run Feature to 3D By Attribute
arcpy.ddd.FeatureTo3DByAttribute(
in_features=os.path.join(feature_dataset_path, fc),
out_feature_class=output_fc_path,
height_field=height_field
)
print(f"Converted {fc} to 3D: {output_fc_name}")
except arcpy.ExecuteError:
print(f"ArcPy error for {fc}: {arcpy.GetMessages(2)}")
except Exception as e:
print(f"General error for {fc}: {e}")
if __name__ == "__main__":
# Define paths and parameters
input_gdb = r"WettedAreaModels_2022.gdb"
feature_dataset = "AllFLows_FDS" # Feature dataset within the geodatabase
input_suffix = "_wetpts" # Input point feature class suffix
output_suffix = "_WSEz" # Output 3D point feature class suffix
height_field = "WSE" # Field used for 3D conversion
print("Starting Feature to 3D By Attribute script...")
# Run the process
feature_to_3d(input_gdb, feature_dataset, input_suffix, output_suffix, height_field)
print("\n--> Finished processing.")
###### STEP 4a - Delete Intermediate Data in Feature Data Set ######
import arcpy
import os
def delete_selected_feature_classes(input_gdb, input_suffix):
"""
Deletes feature classes in the geodatabase that end with a specific suffix.
Parameters:
input_gdb (str): Path to the geodatabase containing feature classes to delete.
input_suffix (str): Suffix of the feature classes to delete (e.g., "_selected").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Delete the feature class
arcpy.management.Delete(fc)
print(f"Deleted feature class: {fc}")
except Exception as e:
print(f"Error deleting {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase
input_gdb = r"AllFlows_FDS"
input_suffix = "_wetpts" # Suffix of feature classes to delete
print("Starting deletion of intermediate files...")
# Run the deletion process
delete_selected_feature_classes(input_gdb, input_suffix)
print("\n--> Finished script successfully.")
###### Step 5 Remove Extraneous Points from non-TA feature classes ######
import arcpy
# Set the workspace to the geodatabase
arcpy.env.workspace = r"WettedAreaModels_2022_A.gdb"
feature_dataset
= "AllFlows_FDS"
# List of ORIG_FID values to delete
orig_fids_to_delete = [101245, 101246, 203705, 255637, 255638, 255639, 443203]
# Debugging: Print workspace and feature dataset info
print("Workspace set to:", arcpy.env.workspace)
desc = arcpy.Describe(feature_dataset)
print(f"Feature Dataset: {desc.name}, Type: {desc.dataType}")
# List all feature classes in the feature dataset
feature_classes = arcpy.ListFeatureClasses(feature_dataset=feature_dataset)
print("Feature Classes in Feature Dataset:", feature_classes)
if feature_classes:
for fc in feature_classes:
# Skip feature classes with "TA" in their names
if "TA" not in fc:
print(f"Processing feature class: {fc}")
try:
# Build SQL query using ORIG_FID
oid_query = f"ORIG_FID IN ({','.join(map(str, orig_fids_to_delete))})"
print("OID Query:", oid_query)
# Create a temporary layer for deletion
arcpy.MakeFeatureLayer_management(fc, "temp_layer", where_clause=oid_query)
# Count matching features
count = arcpy.GetCount_management("temp_layer")[0]
print(f"Matching Features Found: {count}")
if int(count) > 0:
# Delete selected features
arcpy.DeleteFeatures_management("temp_layer")
print(f"Deleted {count} points from {fc}")
else:
print(f"No matching points found in {fc}")
except Exception as e:
print(f"Error processing {fc}: {e}")
finally:
# Always delete the temporary layer to release locks
if arcpy.Exists("temp_layer"):
arcpy.Delete_management("temp_layer")
print("Temporary layer deleted.")
else:
print(f"Skipped feature class: {fc}")
# Compact the geodatabase after processing to ensure locks are released
print("Compacting the geodatabase to remove any lingering locks...")
arcpy.Compact_management(arcpy.env.workspace)
print("All feature classes have been processed.")
else:
print(f"No feature classes found in the feature dataset: {feature_dataset}")
###### Step 6 - MANUAL - Add Hard Line and Soft Clip feature classes to Feature Dataset containing Mass Points ######
# Before proceeding this step was performed manually: Import
# Hard Line (DtDEM_Beyond_OG) and
# Soft Clip (T2022_Terrain_Extent_polygon)
# feature classes for use in terrain creation into the FDS
# (AllFlows_FDS) that contains the set of mass point files.
###### Step 7 - Create Terrains ######
import arcpy
# Enable overwrite for outputs, if necessary
arcpy.env.overwriteOutput = True
# Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")
# Set the path to the feature dataset
feature_dataset = r"AllFlows_FDS"
arcpy.env.workspace
= feature_dataset
# List all feature classes in the feature dataset
feature_classes = arcpy.ListFeatureClasses()
# Ensure feature classes are found
if feature_classes is None or len(feature_classes) == 0:
print(f"No feature classes found in the feature dataset: {feature_dataset}")
else:
# Loop through each feature class and create a terrain with the adjusted naming convention
for fc in feature_classes:
if fc.endswith("WSEz"):
# Remove the last 4 characters ("WSEz") and append "ter"
terrain_name = f"{fc[:-4]}ter"
else:
# Default fallback if the naming convention doesn't match
terrain_name = f"{fc}_ter"
# Create the terrain
arcpy.ddd.CreateTerrain(
in_feature_dataset=feature_dataset,
out_terrain_name=terrain_name,
average_point_spacing=50,
pyramid_type="ZTOLERANCE"
)
print(f"Created terrain: {terrain_name}")
# Release the 3D Analyst license
arcpy.CheckInExtension("3D")
print("Terrain creation process completed.")
###### Step 8 - Add Terrain Pyramid Levels ######
import arcpy
# Enable overwrite for outputs, if necessary
arcpy.env.overwriteOutput = False
# Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")
# Set the feature dataset path
feature_dataset = r"AllFlows_FDS"
arcpy.env.workspace
= feature_dataset
# List all terrains in the feature dataset
terrains = arcpy.ListDatasets("*", "Terrain")
# Ensure terrains are found
if terrains is None or len(terrains) == 0:
print(f"No terrains found in the feature dataset: {feature_dataset}")
else:
# Loop through each terrain and add pyramid levels
for terrain in terrains:
arcpy.ddd.AddTerrainPyramidLevel(
in_terrain=terrain, pyramid_level_definition=["1 1000"]
)
print(f"Added pyramid levels to terrain: {terrain}")
# Release the 3D Analyst license
arcpy.CheckInExtension("3D")
print("Pyramid levels have been added to all terrains.")
###### Step 9 - Add Feature Classes to Terrains ######
import arcpy
def add_feature_classes_to_terrain(output_fds):
"""
Adds corresponding point feature classes as mass points, along with the existing hard line and soft clip feature classes,
to terrains in the feature dataset.
Parameters:
output_fds (str): Path to the output feature dataset containing the terrains and feature classes.
"""
# Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")
try:
# Set workspace to the feature dataset
arcpy.env.workspace = output_fds
# List all point feature classes and terrains
feature_classes = arcpy.ListFeatureClasses(feature_type="Point")
terrains = arcpy.ListDatasets("*", "Terrain")
# Locate the hard line and soft clip feature classes in the FDS
hard_line = "DtDEM_beyond_OG"
soft_clip = "T2022_Terrain_Extent_polygon"
if not feature_classes:
print("No point feature classes found.")
return
if not terrains:
print("No terrains found.")
return
# Process each terrain
for terrain in terrains:
terrain_name = terrain.split("_ter")[0] # Extract base name of terrain
corresponding_fc = f"{terrain_name}" # Assume feature class matches terrain name (without "_ter")
# Ensure the corresponding point feature class exists
if corresponding_fc in feature_classes:
terrain_path = f"{output_fds}{terrain}"
print(f"Processing terrain: {terrain}")
# Add Hard Line feature class
arcpy.ddd.AddFeatureClassToTerrain(
in_terrain=terrain_path,
in_features=[
[hard_line, "SHAPE", "Hard_Line", 3, 0, 1, True, False, "DtDEM_beyond_OG_embed", "<None>", False]
]
)
print(f" - Added Hard Line: {hard_line}")
# Add Soft Clip feature class
arcpy.ddd.AddFeatureClassToTerrain(
in_terrain=terrain_path,
in_features=[
[soft_clip, "<None>", "Soft_Clip", 4, 0, 1, True, False, "T2022_Terrain_Extent_polygon_embed", "<None>", False]
]
)
print(f" - Added Soft Clip: {soft_clip}")
# Add Mass Points feature class
arcpy.ddd.AddFeatureClassToTerrain(
in_terrain=terrain_path,
in_features=[
[corresponding_fc, "Shape", "Mass_Points", 5, 0, 1, True, False, f"{terrain}_embed", "<None>", False]
]
)
print(f" - Added Mass Points: {corresponding_fc}")
else:
print(f"Warning: No corresponding point feature class found for terrain {terrain}")
finally:
# Ensure 3D Analyst extension is checked in
arcpy.CheckInExtension("3D")
print("Feature class addition process completed.")
if __name__ == "__main__":
# Define the output feature dataset
output_fds = r"AllFlows_FDS"
# Call the function
add_feature_classes_to_terrain(output_fds)
###### Step 10 - Build Terrains ######
import arcpy
import os
import time
def build_all_terrains(feature_dataset_path):
# Enable overwriting of outputs if needed
arcpy.env.overwriteOutput = False
# Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")
# Set the workspace
arcpy.env.workspace = feature_dataset_path
# List all terrains in the feature dataset
terrains = arcpy.ListDatasets(feature_type="ALL") # Adjusting to "ALL" as "Terrain" may not be recognized
# Check if terrains exist
if terrains:
print(f"Found {len(terrains)} terrain(s) in the feature dataset:")
for terrain in terrains:
terrain_path = os.path.join(feature_dataset_path, terrain)
print(f"Building terrain: {terrain}")
# Start timing
start_time = time.time()
# Build the terrain (using the 3D Analyst Build Terrain tool)
try:
arcpy.ddd.BuildTerrain(terrain_path)
# Stop timing
elapsed_time = time.time() - start_time
minutes, seconds = divmod(elapsed_time, 60)
hours, minutes = divmod(minutes, 60)
print(f"{terrain} is complete. That took ({int(hours)} hours, {int(minutes)} minutes).")
except arcpy.ExecuteError:
print(f"Failed to build terrain: {terrain}. Error: {arcpy.GetMessages(2)}")
else:
print("No terrains found in the feature dataset.")
if __name__ == "__main__":
# Path to the feature dataset
feature_dataset_path = r"AllFlows_FDS"
# Global Environment settings
with arcpy.EnvManager(scratchWorkspace="P:WettedAreaModels_2022.gdb",
workspace="P:WettedAreaModels_2022.gdb"):
build_all_terrains(feature_dataset_path)
###### Step 10 (alternate) - Build Terrains from List ######
# If processing is interrupted, prompts for terrain path and builds remaining terrains from dynamically generated list
import arcpy
import time
def build_selected_terrains():
# Enable overwriting of outputs if needed
arcpy.env.overwriteOutput = False
# Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")
# Prompt the user to enter terrain paths
print("Enter the full paths to the terrains you want to build, one at a time.")
print("Type 'done' (without quotes) when you are finished entering terrains.")
terrain_paths = []
while True:
terrain_path = input("Enter terrain path: ").strip()
if terrain_path.lower() == 'done':
break
elif arcpy.Exists(terrain_path):
terrain_paths.append(terrain_path)
else:
print(f"Invalid path or terrain does not exist: {terrain_path}")
# Check if any terrains were entered
if not terrain_paths:
print("No valid terrains entered. Exiting the script.")
return
print(f"Found {len(terrain_paths)} terrain(s) to process.")
for terrain_path in terrain_paths:
terrain_name = arcpy.Describe(terrain_path).name # Get terrain name
print(f"Building terrain: {terrain_name}")
# Start timing
start_time = time.time()
try:
# Build the terrain
arcpy.ddd.BuildTerrain(terrain_path)
# Stop timing
elapsed_time = time.time() - start_time
minutes, seconds = divmod(elapsed_time, 60)
hours, minutes = divmod(minutes, 60)
print(f"{terrain_name} is complete. That took ({int(hours)} hours, {int(minutes)} minutes).")
except arcpy.ExecuteError:
print(f"Failed to build terrain: {terrain_name}. Error: {arcpy.GetMessages(2)}")
if __name__ == "__main__":
# Run the function
build_selected_terrains()
###### Step 11 - Surface Difference ######
# Reference Terrain = original DTM Terrain
# Output is stored at GDB level
import arcpy
import os
import time # Import the time module for timing functionality
# Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")
# Path to the feature dataset
feature_dataset = r"AllFlows_FDS"
arcpy.env.workspace
= feature_dataset
# Reference terrain
reference_terrain = r"Terrain_2022"
#
Define the target spatial reference
target_sr = arcpy.SpatialReference(2225, 6360) # WKID 2225 horizontal, 6360 vertical
# List all terrains in the feature dataset
terrains = arcpy.ListDatasets(feature_type="ALL") # Adjust if specific dataset types are needed
# Process each terrain in the feature dataset
if terrains:
print(f"Found {len(terrains)} terrain(s) in the feature dataset.")
for terrain in terrains:
try:
# Print a message before starting the process
print(f"Starting Surface Difference for: {terrain}")
# Start timing
start_time = time.time()
terrain_path = os.path.join(feature_dataset, terrain)
# Extract the name for the output feature class
name = terrain.split("")[-1]
output_fc = f"{name}_SD"
output_path = f"P:{output_fc}"
# Run the Surface Difference tool
arcpy.SurfaceDifference_3d(terrain_path, reference_terrain, output_path)
# Set the coordinate system of the output to the target spatial reference
arcpy.DefineProjection_management(output_path, target_sr)
# End timing
end_time = time.time()
duration = end_time - start_time
# Print a message after completing the process
print(f"Surface Difference completed for: {terrain} in {duration:.2f} seconds.")
except Exception as e:
print(f"Error processing {terrain}: {e}")
else:
print("No terrains found in the feature dataset.")
# Check in the 3D Analyst extension
arcpy.CheckInExtension("3D")
print("Script completed.")
###### Step 11a (alternative) - Surface Difference from dynamic list ######
# If processing is interrupted, prompts for terrain path and creates surface difference for remaining terrains from dynamically generated list
import arcpy
import time # Import the time module for timing functionality
# Check out the 3D Analyst extension
arcpy.CheckOutExtension("3D")
# Prompt user for terrain paths
print("Enter the full paths to the terrains you want to process, one at a time.")
print("Type 'done' (without quotes) when you are finished entering terrains.")
valid_terrains = []
while True:
terrain_path = input("Enter terrain path: ").strip()
if terrain_path.lower() == 'done':
break
elif arcpy.Exists(terrain_path):
valid_terrains.append(terrain_path)
else:
print(f"Invalid path or terrain does not exist: {terrain_path}")
# Reference terrain
reference_terrain = "P:Terrain_2022"
#
Define the target spatial reference
target_sr = arcpy.SpatialReference(2225, 6360) # WKID 2225 horizontal, 6360 vertical
# Process each terrain
if valid_terrains:
for terrain in valid_terrains:
try:
# Print a message before starting the process
print(f"Starting Surface Difference for: {terrain}")
# Start timing
start_time = time.time()
# Extract the name for the output feature class
name = terrain.split("")[-1]
output_fc = f"{name}_SD"
output_path = f"P:{output_fc}"
# Run the Surface Difference tool
arcpy.SurfaceDifference_3d(terrain, reference_terrain, output_path)
# Set the coordinate system of the output to the target spatial reference
arcpy.DefineProjection_management(output_path, target_sr)
# End timing
end_time = time.time()
duration = end_time - start_time
# Print a message after completing the process
print(f"Surface Difference completed for: {terrain} in {duration:.2f} seconds.")
except Exception as e:
print(f"Error processing {terrain}: {e}")
else:
print("No valid terrains provided. Exiting script.")
# Check in the 3D Analyst extension
arcpy.CheckInExtension("3D")
###### Step 12 - Select Wetted Area and Micro-Islands and Export to new Feature Class ######
import arcpy
def select_and_export(input_layer, output_path, where_clause):
"""
Selects features from a feature class based on an attribute query and exports them to a new feature class.
Parameters:
input_layer (str): The input feature class.
output_path (str): The path to the output feature class.
where_clause (str): The SQL query for selecting features.
"""
try:
# Use the Select tool to directly create a new feature class
arcpy.analysis.Select(in_features=input_layer, out_feature_class=output_path, where_clause=where_clause)
print(f"Selection completed for {input_layer}. Exported to {output_path}")
except Exception as e:
print(f"Error processing {input_layer}: {e}")
if __name__ == "__main__":
# Path to the geodatabase
gdb_path = r"WettedAreaModels_2022.gdb"
arcpy.env.workspace = gdb_path
# SQL query for selecting features
where_clause = "Code = 1 Or Code = -1 And Shape_Area <= 5"
# List all feature classes directly in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
print("Feature Classes in Geodatabase:", feature_classes)
if feature_classes:
for fc in feature_classes:
try:
# Define output feature class name
output_fc_name = f"{fc}_selMicIslnds"
output_path = f"{gdb_path}{output_fc_name}"
print(f"Processing feature class: {fc}")
print(f"Output will be saved as: {output_fc_name}")
# Call the function to select and export features
select_and_export(fc, output_path, where_clause)
except Exception as e:
print(f"Error processing feature class {fc}: {e}")
else:
print(f"No feature classes found directly in the geodatabase: {gdb_path}")
print("Script completed successfully!")
###### Step 12a Delete Intermediate Data from GDB ######
import arcpy
import os
def delete_selected_feature_classes(input_gdb, input_suffix):
"""
Deletes feature classes in the geodatabase that end with a specific suffix.
Parameters:
input_gdb (str): Path to the geodatabase containing feature classes to delete.
input_suffix (str): Suffix of the feature classes to delete (e.g., "_selected").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Delete the feature class
arcpy.management.Delete(fc)
print(f"Deleted feature class: {fc}")
except Exception as e:
print(f"Error deleting {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase
input_gdb = r"WettedAreaModels_2022.gdb"
input_suffix = "_SD" # Suffix of feature classes to delete
print("Starting deletion of intermediate files...")
# Run the deletion process
delete_selected_feature_classes(input_gdb, input_suffix)
print("\n--> Finished script successfully.")
###### Step 13 - Dissolve Lines and Absorb MicroIslands ######
import arcpy
# Path to the geodatabase
gdb_path = r"WettedAreaModels_2022.gdb"
def
main():
# Validate the geodatabase
if not arcpy.Exists(gdb_path):
print(f"Error: Geodatabase does not exist: {gdb_path}")
return
# Set the workspace to the geodatabase
arcpy.env.workspace = gdb_path
# Allow overwriting of output
arcpy.env.overwriteOutput = True
# Get a list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print(f"No feature classes found in the geodatabase: {gdb_path}")
return
# Loop through each feature class
for fc in feature_classes:
try:
# Replace the "selMicIslnds" suffix with "Diss" in the output name
if fc.endswith("selMicIslnds"):
output_fc_name = fc.replace("selMicIslnds", "Diss")
else:
output_fc_name = f"{fc}_Diss" # Fallback if the suffix doesn't exist
# Define the full path to the output feature class
output_fc = f"{gdb_path}{output_fc_name}"
# Run the Dissolve tool
arcpy.management.Dissolve(
in_features=fc,
out_feature_class=output_fc,
dissolve_field="", # No fields are specified for dissolving
multi_part="SINGLE_PART" # Ensure multipart features are not created
)
print(f"Successfully dissolved {fc} to {output_fc}")
except Exception as e:
print(f"Error processing {fc}: {e}")
# Print completion message
print("Dissolve operation completed for all feature classes.")
if __name__ == "__main__":
main()
###### Step 13a - Delete Intermediate Data from GDB ######
import arcpy
import os
def delete_selected_feature_classes(input_gdb, input_suffix):
"""
Deletes feature classes in the geodatabase that end with a specific suffix.
Parameters:
input_gdb (str): Path to the geodatabase containing feature classes to delete.
input_suffix (str): Suffix of the feature classes to delete (e.g., "_selected").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Delete the feature class
arcpy.management.Delete(fc)
print(f"Deleted feature class: {fc}")
except Exception as e:
print(f"Error deleting {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase
input_gdb = r"WettedAreaModels_2022.gdb"
input_suffix = "_SD_selMicIslnds" # Suffix of feature classes to delete
print("Starting deletion of intermediate files...")
# Run the deletion process
delete_selected_feature_classes(input_gdb, input_suffix)
print("\n--> Finished script successfully.")
###### Step 14 - Select and Export to Remove Micro-Ponds ######
import arcpy
# Set the workspace
gdb_path = r"WettedAreaModels_2022.gdb"
arcpy.env.workspace
= gdb_path
# List all feature classes in the workspace
feature_classes = arcpy.ListFeatureClasses()
# SQL query for selection
query = "Shape_Area > 5"
# Loop through each feature class
for fc in feature_classes:
try:
# Print message before starting
print(f"Starting selection for feature class: {fc}")
# Define the output feature class name
output_fc = f"{fc}_no_MicPonds"
# Execute the Select tool
arcpy.analysis.Select(in_features=fc, out_feature_class=output_fc, where_clause=query)
# Count the number of selected features
selected_count = arcpy.management.GetCount(output_fc)[0]
# Print message after completing
print(f"Selection completed for feature class: {fc}. Selected {selected_count} features and exported to {output_fc}.")
except Exception as e:
print(f"Error processing feature class {fc}: {e}")
print("Script completed successfully!")
###### Step 14a - Delete Intermediate Data from GDB ######
import arcpy
import os
def delete_selected_feature_classes(input_gdb, input_suffix):
"""
Deletes feature classes in the geodatabase that end with a specific suffix.
Parameters:
input_gdb (str): Path to the geodatabase containing feature classes to delete.
input_suffix (str): Suffix of the feature classes to delete (e.g., "_Diss").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Delete the feature class
arcpy.management.Delete(fc)
print(f"Deleted feature class: {fc}")
except Exception as e:
print(f"Error deleting {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase
input_gdb = r"WettedAreaModels_2022.gdb"
input_suffix = "_SD_Diss" # Suffix of feature classes to delete
print("Starting deletion of intermediate files...")
# Run the deletion process
delete_selected_feature_classes(input_gdb, input_suffix)
print("\n--> Finished script successfully.")
###### Step 15 - Smooth Polygons ######
import arcpy
import os
# Path to the geodatabase
gdb_path = r"WettedAreaModels_2022.gdb"
#
Set the workspace to the geodatabase
arcpy.env.workspace = gdb_path
# Allow overwriting outputs
arcpy.env.overwriteOutput = True
# Define the smoothing parameters
smoothing_algorithm = "PAEK"
smoothing_tolerance = "10 Feet"
topological_errors = "NO_CHECK"
# Get a list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
# Set the output coordinate system (if needed)
output_coordinate_system = arcpy.SpatialReference("NAD 1983 StatePlane California I FIPS 0401 (US Feet)")
# Loop through each feature class
for fc in feature_classes:
# Replace the suffix "SD_Diss_no_MicPonds" with "SmoothPoly" in the output name
if "SD_Diss_no_MicPonds" in fc:
output_fc_name = fc.replace("SD_Diss_no_MicPonds", "SmoothPoly")
else:
output_fc_name = f"{fc}_SmoothPoly" # Fallback if the suffix isn't present
# Define the full path to the output feature class
output_fc = os.path.join(gdb_path, output_fc_name)
# Apply the environment manager for coordinate system
with arcpy.EnvManager(outputCoordinateSystem=output_coordinate_system):
# Run the Smooth Polygon tool
arcpy.cartography.SmoothPolygon(
in_features=fc,
out_feature_class=output_fc,
algorithm=smoothing_algorithm,
tolerance=smoothing_tolerance,
error_option=topological_errors
)
print(f"Smoothed {fc} to {output_fc}")
print("Smoothing operation completed for all feature classes.")
###### Step 15a - Delete Intermediate Data from GDB ######
import arcpy
import os
def delete_selected_feature_classes(input_gdb, input_suffix):
"""
Deletes feature classes in the geodatabase that end with a specific suffix.
Parameters:
input_gdb (str): Path to the geodatabase containing feature classes to delete.
input_suffix (str): Suffix of the feature classes to delete (e.g., "_MicPonds").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Delete the feature class
arcpy.management.Delete(fc)
print(f"Deleted feature class: {fc}")
except Exception as e:
print(f"Error deleting {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase
input_gdb = r"WettedAreaModels_2022.gdb"
input_suffix = "_MicPonds" # Suffix of feature classes to delete
print("Starting deletion of intermediate files...")
# Run the deletion process
delete_selected_feature_classes(input_gdb, input_suffix)
print("\n--> Finished script successfully.")
###### Step 16 - MANUAL - Create Erase Feature for Trimming Smoothed Polygons to D/S Mesh Boundary ######
# The boundary of the SRH2D Mesh Polygon was used to visually locate the extent of the mass points at the downstream end of the wetted area # and based on that, a copy of the largest and widest smoothed polygon (based on11000cfs TA MFF extent) was trimmed using the Split Editing tool in ArcGIS Pro. # The polygon section upstream of the split was deleted, and the remaining pieces merged, creating a single polygon which was used as an Erase feature to "trim"
# the smoothed polygons for all flows to end of the d/s mesh boundary. ###### Step 17 - Trim Wetted Area Polygons ######
import arcpy
# Paths to the geodatabases
input_gdb_path = r"WettedAreaModels_2022.gdb"
erase_feature
= r"cut_poly"
#
Set workspace to the input geodatabase
arcpy.env.workspace = input_gdb_path
# List all feature classes in the input geodatabase
feature_classes = arcpy.ListFeatureClasses()
# Output folder (same as input GDB in this case)
output_path = input_gdb_path
# Check if feature classes were found
if not feature_classes:
print("No feature classes found in the geodatabase.")
else:
# Iterate through each feature class in the input GDB
for fc in feature_classes:
# Skip feature classes that do not end in "_SmoothPoly"
if not fc.endswith("_SmoothPoly"):
print(f"Skipping feature class: {fc}")
continue # Move to the next feature class
# Define the output feature class name
output_fc_name = fc.replace("_SmoothPoly", "") # Remove "_SmoothPoly"
output_fc = f"{output_path}{output_fc_name}"
print(f"Processing: {fc}")
try:
# Perform the Erase operation
arcpy.analysis.Erase(in_features=fc, erase_features=erase_feature, out_feature_class=output_fc)
print(f"Successfully created: {output_fc}")
except Exception as e:
print(f"Failed to process {fc}: {e}")
print("All feature classes processed.")
###### Step 17a - Delete Intermediate Data from GDB ######
import arcpy
import os
def delete_selected_feature_classes(input_gdb, input_suffix):
"""
Deletes feature classes in the geodatabase that end with a specific suffix.
Parameters:
input_gdb (str): Path to the geodatabase containing feature classes to delete.
input_suffix (str): Suffix of the feature classes to delete (e.g., "_SmoothPoly").
"""
arcpy.env.workspace = input_gdb
arcpy.env.overwriteOutput = True
# Get the list of all feature classes in the geodatabase
feature_classes = arcpy.ListFeatureClasses()
if not feature_classes:
print("No feature classes found in the geodatabase.")
return
print(f"Found {len(feature_classes)} feature classes to process.")
# Iterate over feature classes with the specific suffix
for fc in feature_classes:
if fc.endswith(input_suffix):
try:
# Delete the feature class
arcpy.management.Delete(fc)
print(f"Deleted feature class: {fc}")
except Exception as e:
print(f"Error deleting {fc}: {e}")
if __name__ == "__main__":
# Define the input geodatabase
input_gdb = r"WettedAreaModels_2022.gdb"
input_suffix = "_SmoothPoly" # Suffix of feature classes to delete
print("Starting deletion of intermediate files...")
# Run the deletion process
delete_selected_feature_classes(input_gdb, input_suffix)
print("\n--> Finished script successfully.")
###### Step 18 - Remove Model Artifacts ######
# MANUAL
# Despite the edits to the boundary polygon and the addition of points to the Boundary points data, a few modeling artifacts remained. At locations where the model formed a pond within the tributary # channels or their floodplains, plus at one location in Junction City near the firehouse where the water surface continued to associate inappropriately with upstream points due to a sharp bend in the river, but # bedrock would likely block groundwater from reaching that elevation. Polygons with ponding that did not connect to the river at Bucktail (at flows <= 500 cfs) and Lowden Meadows (TA flows) and Junction # City (flows < 1000 cfs) were removed. Above 500 cfs the Bucktail polygons and above1000 cfs, the Junction City polygons became small ponds that appeared more likely tp be correct and including # association with parallel or downstream data.
SRH2D flow modeling output (Bradley 2022). Specifically, the Habitat feature class mesh polygons were used to generate mass points and provide water depth and water surface elevation for All FLows
Trinity River 40 Mile Hydraulic Model: Update with 2022 Topography
2018-02-28
D. Nathan Bradley
U.S. Bureau of Reclamation, Technical Service Center, Sedimentation and River Hydraulics Group
Technical Report No. SRH-2018-11
2022 Topographic Model of the Trinity River from Lewiston Dam to the North Fork Trinity River
Data Package for 2022Trinity River Focal Reach Digital Terrain
Eric Peterson
Trinity River Restoration Program, U.S. Bureau of Reclamation
Data Steward
530-623-1800
Weaverville
CA
96093
US
ebpeterson@usbr.gov
GMA (GMA Hydrology, Inc.). 2017. Data Package for 2016 Trinity River Focal Reach Digital Terrain. Minor correction made September 1. Data Package for the Trinity River Restoration Program (TRRP) under U.S. Bureau of Reclamation contract R14PC00122. GMA Hydrology, Arcata, California.
A polygon with a single elevation was imported from previous terrain creation dataset for inclusion in WSE terrain dataset.
Modifications performed by Eric Peterson.
A polygon to be used as an Erase feature was created for trimming wetted area polygons to D/S extend of the SRH2D Mesh boundary
The boundary of the SRH2D Mesh Polygon was used to visually locate the extent of the mass points at the downstream end of the wetted area and based on that, a copy of the largest and widest smoothed polygon (based on11000cfs TA MFF extent) was trimmed using the Split Editing tool in ArcGIS Pro. The polygon section upstream of the split was deleted, and the remaining pieces merged, creating a single polygon which was used as an Erase feature to "trim"
the smoothed polygons for all flows to end of the d/s mesh boundary.
Step 01. (scripted)
Select wetted polygons (WaterDepth>0) and export to a new feature class
Select Wetted Cells from SRH2D Mesh Cells
Polygons with water depth greater than zero (therefore wetted) were selected and exported to a temporary file within a GDB.
Creation of water surface elevation data
2025-01-14T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 02. (scripted)
Export wetted polygons using field mapping to retain only relevant fields - Reduce fields to only include WSE
Export Feature Classes using field mapping to reduce the number of fields for faster processing, retaining only the water surface elevation (WSE)( field, which will be used in the next step. Step 02a. (scripted)
Delete intermediate data from GDB.
2025-01-15T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 03. (scripted)
Export wetted cell centroids as point feature classes in a Feature Dataset
Data preparation: Create Wetted Points from Mesh Polygon Cell Centroids
Export wetted cell centroids as point feature classes in a Feature Dataset using Feature to Point
Step 03a. (scripted)
Delete intermediate data from GDB.
Creation of the Water Surface Elevation data.
2025-01-16T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 04. (scripted)
Assign WSE value as Z for all point feature classes in FDS using Feature to 3D By Attribute
Data preparation: Create 3D Points using WSE as Z
Assign WSE value as Z for all point feature classes in FDS using Feature to 3D By Attribute, Height = WSE
Step 04a. (scripted)
Delete intermediate data from GDB.
Creation of the Water Surface Elevation data.
2025-01-17T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 05. (scripted)
Remove Extraneous Points from non-TA feature classes
Preparation of mass point files for inclusion in terrains included manual removal of a few modeling artifact polygon centroids where major tributaries approach the river and at the tight river bend in Junction City. Point feature classes for flows that do not include tributary accretions end up with several extraneous points not connected to the mainstem. Script deleted visually identified points with Orig_FID numbers 101245, 101246, 203705, 255637, 255638, 255639, 443203.
Remove artifacts prior to terrain creation
2025-01-18T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 06. (manual)
Add Hard Line and Soft Clip feature classes to Feature Dataset containing Mass Points
Before proceeding this step was performed manually imported feature classes into the FDS:
Hard Line = DtDEM_Beyond_OG Soft Clip = T2022_Terrain_Extent_polygon
These are feature classes that will be used in terrain creation.
AllFlows_FDS is the FDS that contains the set of all mass point files created in the steps above.
Step necessary for terrain creation
2025-01-19T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 07. (scripted)
Create Terrains
The WSE surface (terrain format) was created within the temporary feature dataset. The script iterated through the list of point feature classes prepared in the steps above and created an empty terrain for each flow using ArcGIS Create Terrain tool.
This provides a water surface extrapolated outward from the actual river.
2025-01-20T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 08. (scripted)
Add Terrain Pyramid Levels
Adds pyramid level 1 1000 to each terrain
Necessary step in terrain building process
2025-01-21T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 09. (scripted)
Add Feature Class to Terrain
This script ran the ArcGIS Add Feature Class To Terrain tool to add 3 feature classes to each of the terrains described above, with the Terrain_Extent (boundar)y polygon as an outward soft clip and the buffered boundary polygon with a single elevation (DtDEM_Beyond_OG) as hard line in every instance. Thee SRH2D wetted points (from Step 02) with WSE = Z elevations for each flow added to the corresponding terrain as mass points
Necessary step in terrain building process
2025-01-22T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 10. (scripted)
Build Terrains
Script iterates through all terrains in the FDS and runs ArcGIS Build Terrain tool.*
Step 10a. (alternative) Build Terrains from List If processing is interrupted, prompts for path and builds remaining terrains from dynamically generated list so processing can continue
from where it left off.
*It is recommended that terrain properties be checked after interruption to ensure that all built terrains are valid.
Necessary step in terrain building process.
2025-01-23T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 11. (scripted)
Surface Difference
The WSE surface was intersected with the 2022 terrain surface using the ArcGIS Surface Difference tool. Polygons flagged with a positive value are where the WSE surface is above the terrain surface.
Step 11a. (alternative)
Surface Difference from List
If processing is interrupted, prompts for path and builds remaining terrains from dynamically generated list so processing can continue
from where it left off.
Locates where the water surface is above ground.
2025-01-24T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 12. (scripted)
Select wetted area and micro-islands and export
Filters to Wetted Area. Selected all features from items processed in the last step, including above-ground polygons with shape area >5 ft 2 and exported to a new FC. Select tool used to select rough wetted area and "micro-islands" and exported to a new feature class. Surface difference (rough wetted area) results were filtered to retain only polygons flagged as above ground, plus polygons below ground with an area of <= 5 ft2. ( "Code = 1 Or Code = -1 And Shape_Area <= 5")
First step in filtering to provide a wetted area without excessive complexity from topographic noise. Inclusion of the small below-ground polygons ensures that tiny islands caused by noise in the terrain model data are omitted from the final product.
2025-01-25T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 13 (scripted)
Dissolve lines and absorb Micro-Islands
The dissolve tool was run on the results of the last step to absorb micro-islands (small below-ground polygons with shape area (<=5 ft2) and remove lines to create a single-part polygon
Polygons were dissolved without creating multipart features.
This was a necessary step for data maintenance, as the surface difference tool breaks the dataset into chunks for processing. Each chunk is left in the results as separate polygons. Without this step, there would be an appearance of seamlines at various intervals along the river.
2025-01-26T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 14 (scripted)
Select and Export to Remove MicroPonds
Second step in filtering to provide a wetted area without excessive complexity from topographic noise. Inclusion of the small (shape area = < 5sft) polygons ensures that tiny "micro-ponds"s caused by noise in the terrain model data are omitted from the final product.
Step 14a (scripted)
Delete Intermediate Data from GDB
Removes model artifacts. This ensures that tiny puddles caused by noise in the terrain model data are omitted from the final product.
2025-01-27T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 15. (scripted)
Smooth Polygons
Polygons were smoothed with the ArcGIS Smooth Polygon tool using the PAEK algorithm and a window of 10.
Step 15a (scripted)
Delete Intermediate Data from GDB
This alters harsh angular forms of the wetted area polygon (due to TIN facets of the topographic Terrain data) to a more organic form that is realistic to true water edges. This processing was reviewed compared to high resolution aerial photography (0.5 ft pixels) and found to be as qualitatively accurate, or better, than the angular data from the TIN facets.
2025-01-28T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 16 (manual)
Create Erase Feature for Trimming Smoothed Polygons to D/S Mesh Boundary
The boundary of the SRH2D Mesh Polygon was used to cut off the end of the largest and widest smoothed polygon (based on11000cfs MFF extent). The remaining polygon piece was
was exported to the GDB containing the smoothed polygons and used as an Erase feature to "trim" smoothed polygons for all flows, so that they do not extend past the end of the d/s mesh boundary.
Creates polygons to be used as Erase Feature in next step
2025-01-29T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 17 (scripted)
Trim Wetted Area Polygons Use Erase to trim ends of wetted area polygons to mass point extent using"cut_polygons" feature class created in previous step as Erase Feature
17a (scripted)
Delete Intermediate Data from GDB
Ensures that wetted area polygons extent does not exceed model area
2025-01-30T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
Step 18 (manual)
Remove (Clear) Model Artifacts
Despite the edits to the boundary polygon and the addition of points to the Boundary points data, a few modeling artifacts remained. At locations where the model formed a pond within the tributary channels or their floodplains, plus at one location in Junction City near the firehouse where the water surface continued to associate inappropriately with upstream points due to a sharp bend in the river, but bedrock would likely block groundwater from reaching that elevation. Polygons with ponding that did not connect to the river at Hatchery (at flows <= 300 cfs) Bucktail (at flows <= 500 cfs) and Lowden Meadows (for TA flows) and Junction City (at flows <= 1000 cfs) were removed. Above 500 cfs the Bucktail polygons and above 1000 cfs, the Junction City polygons became small ponds that appeared more likely to be correct and including association with parallel or downstream data.
Cleaning of model artifacts.
2025-01-31T00:00:00
Jeanne McSloy
Trinity River Restoration Program
Natural Resource Specialist
0
EPSG
5.3(9.0.0)
6.18.3(9.3.1.2)
Estimated wetted area.
Simple
FALSE
0
TRUE
FALSE