Flood Plain Boundaries

This section defines the spatial boundaries of the flood plain within the model domain and establishes which model grid cells are considered “active” (within the groundwater domain) for the simulation.

Importing Libraries and Retrieving Variables

All necessary libraries are imported to support geospatial analysis and model configuration. Key variables, model parameters, and spatial data—including previously processed boundary lines and model grid information—are retrieved using the %store magic command.

Show code cell source Hide code cell source

# Import Libraries
%run /Users/u4eeevmq/Documents/Python/HyporheicFloPy/VQuintana/notebooks/common_imports.py

# Retrieve stored variables
%store -r md6_exe_path
%store -r md7_exe_path
%store -r sim_name
%store -r workspace
%store -r figs_path
%store -r gwf_name
%store -r mp7_name
%store -r gwf_ws
%store -r mp7_ws
%store -r headfile
%store -r head_filerecord
%store -r budgetfile
%store -r budget_filerecord
%store -r write
%store -r run
%store -r plot
%store -r plot_show
%store -r plot_save

# Retrieve model parameters
%store -r length_units
%store -r time_units
%store -r nper
%store -r cell_size_x
%store -r cell_size_y
%store -r gw_mod_depth
%store -r z
%store -r kh
%store -r kv
%store -r gw_offset
%store -r porosity
%store -r rch_iface
%store -r rch_iflowface
%store -r recharge_rate
%store -r nstp
%store -r perlen
%store -r tsmult

# Retrieve spatial data
%store -r hec_ras_crs
%store -r terrain_elevation
%store -r raster_transform
%store -r transform
%store -r raster_crs
%store -r terrain_output_raster
%store -r water_surface_output_raster
%store -r ground_water_domain
%store -r left_boundary
%store -r right_boundary

# retrieve model domain data
%store -r terrain_elevation
%store -r raster_transform
%store -r raster_crs
%store -r raster_bounds_box
%store -r bed_elevation
%store -r raster_width
%store -r raster_height
%store -r ncol
%store -r nrow
%store -r top
%store -r nlay
%store -r grid_x
%store -r grid_y
%store -r grid_points
%store -r intersecting_points
%store -r xorigin
%store -r yorigin
%store -r xmin
%store -r ymin
%store -r xmax
%store -r ymax
%store -r tops
%store -r botm

Defining the Flood Plain Boundaries

We begin by extracting the start and end coordinates from both the left and right boundary lines of the flood plain, which are typically derived from HEC-RAS or other GIS-based delineations. These coordinates are used to define the upstream and downstream boundaries of the flood plain. The upstream boundary is created by connecting the starting points of the left and right boundaries, while the downstream boundary connects their respective endpoints.

These boundaries are then represented as LineString objects, which are converted into GeoDataFrames for further spatial analysis and visualization. This formalizes the spatial limits of the flood plain for subsequent modeling steps.

#--------------------------- Define Flood Plain Boundaries -----------------------------#
## Extract the first (start) and last (end) coordinates from each boundary geometry
left_start = left_boundary.geometry.iloc[0].coords[0]  # First point
left_end = left_boundary.geometry.iloc[-1].coords[-1]  # Last point

right_start = right_boundary.geometry.iloc[0].coords[0]  # First point
right_end = right_boundary.geometry.iloc[-1].coords[-1]  # Last point

## Print start and end coordinates
print(f"Left Boundary Start: {left_start}, Left Boundary End: {left_end}")
print(f"Right Boundary Start: {right_start}, Right Boundary End: {right_end}")

## Upstream boundary coordinates (use left and right start points)
upstream_start_x, upstream_start_y = left_start  # Start of upstream from left boundary
upstream_end_x, upstream_end_y = right_start  # End of upstream from right boundary

## Downstream boundary coordinates (use left and right end points)
downstream_start_x, downstream_start_y = left_end  # Start of downstream from left boundary
downstream_end_x, downstream_end_y = right_end  # End of downstream from right boundary

## Print extracted coordinates
print(f"Upstream Start: ({upstream_start_x}, {upstream_start_y}), Upstream End: ({upstream_end_x}, {upstream_end_y})")
print(f"Downstream Start: ({downstream_start_x}, {downstream_start_y}), Downstream End: ({downstream_end_x}, {downstream_end_y})")

## Use the extracted start and end points
upstream_line = LineString([(upstream_start_x, upstream_start_y), (upstream_end_x, upstream_end_y)])
downstream_line = LineString([(downstream_start_x, downstream_start_y), (downstream_end_x, downstream_end_y)])

## Convert to GeoDataFrame
upstream_boundary = gpd.GeoDataFrame(geometry=[upstream_line], crs=left_boundary.crs)
downstream_boundary = gpd.GeoDataFrame(geometry=[downstream_line], crs=left_boundary.crs)

## Print boundary lines
print(f"Upstream Line: {upstream_line}")
print(f"Downstream Line: {downstream_line}")

## Check if boundaries were created successfully
if not upstream_boundary.empty and not downstream_boundary.empty:
    print("✅ Upstream and Downstream boundaries created successfully!")
else:
    print("❌ Error: One or both boundaries are empty. Check input data.")

Left Boundary Start: (2407358.236295765, 10515701.162934666), Left Boundary End: (2407570.6507961783, 10515177.034472264)
Right Boundary Start: (2407065.5187525107, 10515580.276633618), Right Boundary End: (2407188.132000717, 10515108.82005953)
Upstream Start: (2407358.236295765, 10515701.162934666), Upstream End: (2407065.5187525107, 10515580.276633618)
Downstream Start: (2407570.6507961783, 10515177.034472264), Downstream End: (2407188.132000717, 10515108.82005953)
Upstream Line: LINESTRING (2407358.236295765 10515701.162934666, 2407065.5187525107 10515580.276633618)
Downstream Line: LINESTRING (2407570.6507961783 10515177.034472264, 2407188.132000717 10515108.82005953)
✅ Upstream and Downstream boundaries created successfully!

Defining Active Cells Using the Groundwater Domain

To accurately represent the model domain, we generate polygons for each model grid cell based on the established grid resolution. Each cell is checked for intersection with the groundwater domain polygon (representing the modeled flood plain). Through a spatial join, we determine which cells fall inside the groundwater domain.

An idomain array is initialized, marking all cells as inactive by default. Cells that intersect the groundwater domain are marked as active across all model layers. This step is essential for MODFLOW 6, as it determines where groundwater flow is simulated within the model grid.

#--------------------------- Define Active Cells using Groundwater Domain -----------------------------#
## Define grid cell polygons based on grid resolution
grid_cells = []
for row in range(nrow):
    for col in range(ncol):
        x_min = grid_x[row, col] - (cell_size_x / 2)
        x_max = grid_x[row, col] + (cell_size_x / 2)
        y_min = grid_y[row, col] - (cell_size_y / 2)
        y_max = grid_y[row, col] + (cell_size_y / 2)
        grid_cells.append(Polygon([(x_min, y_min), (x_min, y_max), (x_max, y_max), (x_max, y_min)]))

## Convert to GeoDataFrame
grid_gdf = gpd.GeoDataFrame(geometry=grid_cells, crs=ground_water_domain.crs)

## Perform Spatial Join
grid_gdf["inside_domain"] = grid_gdf.geometry.intersects(ground_water_domain.unary_union)

## Initialize IDOMAIN array
idomain = np.zeros((nlay, nrow, ncol), dtype=int)

## Assign active cells where grid intersects groundwater domain
for idx, inside in enumerate(grid_gdf["inside_domain"]):
    row, col = divmod(idx, ncol)  # Convert flat index to row, col
    if inside:
        idomain[:, row, col] = 1  # Mark as active

## Debugging: Print active/inactive cell count
print(f"✅ Total Active Cells: {np.sum(idomain == 1)}")
print(f"✅ Total Inactive Cells: {np.sum(idomain == 0)}")

✅ Total Active Cells: 317960
✅ Total Inactive Cells: 6319440

C:\Users\u4eeevmq\AppData\Local\Temp\ipykernel_28112\709847143.py:16: DeprecationWarning: The 'unary_union' attribute is deprecated, use the 'union_all()' method instead.
  grid_gdf["inside_domain"] = grid_gdf.geometry.intersects(ground_water_domain.unary_union)

Show code cell source Hide code cell source

# Store new variables
%store left_start
%store left_end
%store right_start
%store right_end
%store upstream_start_x
%store upstream_start_y
%store upstream_end_x
%store upstream_end_y
%store downstream_start_x
%store downstream_start_y
%store downstream_end_x
%store downstream_end_y
%store upstream_line
%store downstream_line
%store upstream_boundary
%store downstream_boundary
%store grid_cells
%store grid_gdf
%store idomain

Stored 'left_start' (tuple)
Stored 'left_end' (tuple)
Stored 'right_start' (tuple)
Stored 'right_end' (tuple)
Stored 'upstream_start_x' (float)
Stored 'upstream_start_y' (float)
Stored 'upstream_end_x' (float)
Stored 'upstream_end_y' (float)
Stored 'downstream_start_x' (float)
Stored 'downstream_start_y' (float)
Stored 'downstream_end_x' (float)
Stored 'downstream_end_y' (float)
Stored 'upstream_line' (LineString)
Stored 'downstream_line' (LineString)
Stored 'upstream_boundary' (GeoDataFrame)
Stored 'downstream_boundary' (GeoDataFrame)

Stored 'grid_cells' (list)

Stored 'grid_gdf' (GeoDataFrame)
Stored 'idomain' (ndarray)