Model Domain#
This section focuses on defining the spatial extent, grid resolution, and vertical structure of the groundwater model domain. The model domain setup leverages geospatial data to ensure that the model grid accurately represents the physical terrain.
Importing Libraries and Retrieving Variables#
Before constructing the model domain, we import required libraries for geospatial data processing, numerical analysis, and visualization. We also retrieve all previously stored variables and model parameters, such as simulation names, model settings, and spatial data.
Show code cell source
#-----------------------Importing Libraries-----------------------#
# import os
# import flopy
# import pathlib
# import matplotlib.pyplot as plt
# import numpy as np
# import geopandas as gpd
# import pandas as pd
# import rasterio
# import pyproj
# import shutil
# import random
# import jupyter_book
# import pathlib as pl
# from pathlib import Path
# from rasterio.crs import CRS
# from rasterio.plot import show
# from rasterio.warp import calculate_default_transform, reproject, Resampling
# from rasterio.transform import from_bounds
# from rasterio.transform import rowcol
# from rasterio.mask import mask
# from shapely.geometry import box, Point, Polygon, LineString
# from flopy.utils.binaryfile import HeadFile
# from scipy.interpolate import griddata
# from pprint import pformat
# from flopy.plot.styles import styles
# from matplotlib.lines import Line2D
# from flopy.mf6 import MFSimulation
# from matplotlib import cbook, cm
# from matplotlib.colors import LightSource
# from modflow_devtools.misc import get_env, timed
# 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
Reading Terrain Raster and Calculating Grid Parameters#
The terrain elevation raster, previously processed and reprojected, forms the basis for defining the model grid. We read the raster data to extract elevation values, calculate the bed elevation (minimum surface value), and determine the spatial extent (width and height) of the raster. The number of model columns (ncol) and rows (nrow) are then calculated based on the cell size parameters.
#------------------------------Model Grid--------------------------------#
## Plot Terrain Elevation
reprojected_terrain_elevation_raster = terrain_output_raster
# Read raster data and extract elevation values
with rasterio.open(reprojected_terrain_elevation_raster) as src:
raster_array = src.read(1) # Read the first band
raster_transform = src.transform
raster_crs = src.crs
raster_bounds_box = box(*src.bounds) # Create a shapely box for raster bounds
terrain_elevation = np.ma.masked_equal(raster_array, src.nodata) # Mask no-data values
# Find Bed Elevations
bed_elevation = np.min(terrain_elevation) # Bed elevation is the minimum value of the cropped surface elevation
# Print the bed elevation
print(f"Bed elevation: {bed_elevation}")
# Calculate the extent of the raster
raster_width = src.width * transform[0] # Width in feet
raster_height = src.height * abs(transform[4]) # Height in feet
# Determine the number of rows and columns based on the cell size
ncol = int(raster_width / cell_size_x)
nrow = int(raster_height / cell_size_y)
# Print the calculated grid dimensions
print(f'Number of columns: {ncol}')
print(f'Number of rows: {nrow}')
top = np.full((nrow, ncol), gw_mod_depth) # feet below bed surface (depth of groundwater model)
nlay = int(top.max() / z) # Number of groundwater layers based on default depth)
# Simulation object
cfg = SimpleNamespace(
sim_name= sim_name,
gwf_ws= gwf_ws
)
Bed elevation: 1599.8125
Number of columns: 431
Number of rows: 385
Initializing Grid and Model Layers#
With the grid dimensions established, we generate arrays representing the center coordinates of each model cell. These coordinates are stored in a GeoDataFrame for spatial analyses. The vertical structure of the model is defined by calculating the top elevation of the uppermost layer from the surface raster, while subsequent layer elevations are computed based on the specified layer thickness.
# Calculate grid cell centers
grid_x, grid_y = np.meshgrid(
np.arange(0, ncol) * cell_size_x + (cell_size_x / 2),
np.arange(0, nrow) * cell_size_y + (cell_size_y / 2),
)
# Convert cell centers to Points for intersection checks
grid_points = gpd.GeoDataFrame(
{"geometry": [Point(x, y) for x, y in zip(grid_x.ravel(), grid_y.ravel())]},
crs= hec_ras_crs, # Replace with the actual CRS of your grid
)
# Read raster data and extract elevation values
with rasterio.open(reprojected_terrain_elevation_raster) as src:
raster_array = src.read(1) # Read the first band
raster_transform = src.transform
raster_crs = src.crs
raster_bounds_box = box(*src.bounds) # Create a shapely box for raster bounds
terrain_elevation = np.ma.masked_equal(raster_array, src.nodata) # Mask no-data values
# Create a GeoDataFrame for raster bounds
raster_bounds_gdf = gpd.GeoDataFrame(
{"geometry": [raster_bounds_box]}, crs=hec_ras_crs
)
# Reproject grid points to match raster CRS
# Use the raster bounds to define the grid extent
minx, miny, maxx, maxy = raster_bounds_box.bounds
grid_x, grid_y = np.meshgrid(
np.linspace(minx, maxx, ncol),
np.linspace(miny, maxy, nrow),
)
# Recreate the grid points
grid_points = gpd.GeoDataFrame(
{"geometry": [Point(x, y) for x, y in zip(grid_x.ravel(), grid_y.ravel())]},
crs=raster_crs,
)
# Check intersection between grid points and raster bounds
intersecting_points = grid_points[grid_points.geometry.intersects(raster_bounds_box)]
# Debugging: Print details about the GeoDataFrames
print(f"Raster CRS: {raster_bounds_gdf.crs}")
print(f"Grid CRS: {grid_points.crs}")
print(f"Number of grid points: {len(grid_points)}")
print(f"Number of intersecting points: {len(intersecting_points)}")
# Set x and y origin
xmin, ymin, xmax, ymax = raster_bounds_box.bounds # Extract bounding box extent
xorigin = xmin # Set xorigin to the left-most boundary
yorigin = ymin # Set yorigin to the bottom-most boundary
# Store xmin, ymin, xmax, ymax
%store xmin
%store ymin
%store xmax
%store ymax
# Extract raster extent before looping
transform = raster_transform
xmin = transform.c
ymax = transform.f
xmax = xmin + (terrain_elevation.shape[1] * transform.a)
ymin = ymax + (terrain_elevation.shape[0] * transform.e)
print(f"✅ Raster Extent: X = ({xmin}, {xmax}), Y = ({ymin}, {ymax})")
Raster CRS: PROJCS["NAD_1983_StatePlane_Texas_Central_FIPS_4203_Feet",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",2296583.33333333],PARAMETER["False_Northing",9842500.0],PARAMETER["Central_Meridian",-100.333333333333],PARAMETER["Standard_Parallel_1",30.1166666666667],PARAMETER["Standard_Parallel_2",31.8833333333333],PARAMETER["Latitude_Of_Origin",29.6666666666667],UNIT["US survey foot",0.304800609601219]]
Grid CRS: PROJCS["NAD83 / Texas Central (ftUS)",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.257222101004,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4269"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["latitude_of_origin",29.6666666666667],PARAMETER["central_meridian",-100.333333333333],PARAMETER["standard_parallel_1",30.1166666666667],PARAMETER["standard_parallel_2",31.8833333333333],PARAMETER["false_easting",2296583.33333333],PARAMETER["false_northing",9842500],UNIT["US survey foot",0.304800609601219,AUTHORITY["EPSG","9003"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
Number of grid points: 165935
Number of intersecting points: 165935
Stored 'xmin' (float)
Stored 'ymin' (float)
Stored 'xmax' (float)
Stored 'ymax' (float)
✅ Raster Extent: X = (2406093.5149272773, 2408251.1279889615), Y = (10514475.360625861, 10516401.62506268)
Handling Missing Data and Layer Generation#
Any missing or invalid elevation data in the top layer are interpolated using spatial interpolation techniques to ensure a continuous surface. We then iteratively generate the tops and bottoms of all groundwater model layers, maintaining a consistent thickness for each layer.
# Function to interpolate NA values
def interpolate_na(terrain):
# Get the coordinates of the non-masked values
valid_mask = ~terrain.mask
valid_coords = np.array(np.nonzero(valid_mask)).T
valid_values = terrain[valid_mask]
# Get the coordinates of the masked values
invalid_mask = terrain.mask
invalid_coords = np.array(np.nonzero(invalid_mask)).T
# Interpolate the values at the masked coordinates
interpolated_values = griddata(valid_coords, valid_values, invalid_coords, method='nearest')
# Fill the masked values with the interpolated values
terrain[invalid_mask] = interpolated_values
return terrain
# Initialize the top array
top = np.full((nrow, ncol), np.nan)
# Update "top" values for each cell in the first layer based on surface elevation
for i in range(nrow):
for grid_col in range(ncol):
# Calculate the x, y coordinates of the cell center
point_x = grid_x[i, grid_col]
point_y = grid_y[i, grid_col]
# Convert the grid cell center coordinates to raster indices
col, row = ~raster_transform * (point_x, point_y)
col, row = int(col), int(row)
# Check if the indices are within raster bounds
if 0 <= row < terrain_elevation.shape[0] and 0 <= col < terrain_elevation.shape[1]:
elevation_value = terrain_elevation[row, col]
# Update "top" based on the raster value
top[i, grid_col] = elevation_value
# Interpolate any remaining NA values in the top array
top = interpolate_na(np.ma.masked_invalid(top))
## Initialize `tops` and `botm` lists
tops = [top] # Add the top layer (surface elevation or default)
botm = []
# First layer bottom is calculated from the updated "top" values
first_layer_botm = np.full_like(top, bed_elevation) # Subtract 0.5 ft for the first layer
botm.append(first_layer_botm)
# Create remaining layers with a constant thickness of 0.5 ft
for layer in range(1, 40): # Layers 2 to 40
next_layer_top = botm[-1] # The top of the current layer is the bottom of the previous layer
next_layer_botm = next_layer_top - z # Subtract thickness from the top
tops.append(next_layer_top) # Add the top of the current layer
botm.append(next_layer_botm) # Add the bottom of the current layer
# Debugging: Check the top and bottom elevations for all layers
print("Top layer elevation (max, min):", tops[0].max(), botm[0].min())
for layer_idx in range(len(tops)):
print(f"Layer {layer_idx + 1} top (max, min):", tops[layer_idx].max(), tops[layer_idx].min())
print(f"Layer {layer_idx + 1} botm (max, min):", botm[layer_idx].max(), botm[layer_idx].min())
## Visualization: Plot the top and last bottom layers
fig, axs = plt.subplots(1, 2, figsize=(15, 8))
# Plot the top layer elevation
im1 = axs[0].imshow(top, cmap="terrain", interpolation="nearest", origin="lower") # 🔹 Ensures (0,0) is bottom-left
axs[0].set_title("Top Layer Elevation")
axs[0].set_xlabel("Column")
axs[0].set_ylabel("Row")
fig.colorbar(im1, ax=axs[0], label="Elevation (ft)")
# Plot the bottom layer elevation (Layer 40)
im2 = axs[1].imshow(botm[-1], cmap="terrain", interpolation="nearest", origin="lower") # 🔹 Ensures bottom-left
axs[1].set_title("Bottom Layer Elevation (Layer 40)")
axs[1].set_xlabel("Column")
axs[1].set_ylabel("Row")
fig.colorbar(im2, ax=axs[1], label="Elevation (ft)")
plt.tight_layout()
plt.show()
C:\Users\u4eeevmq\AppData\Local\Temp\ipykernel_14668\2273222409.py:39: UserWarning: Warning: converting a masked element to nan.
top[i, grid_col] = elevation_value
Top layer elevation (max, min): 1642.875 1599.8125
Layer 1 top (max, min): 1642.875 1599.875
Layer 1 botm (max, min): 1599.8125 1599.8125
Layer 2 top (max, min): 1599.8125 1599.8125
Layer 2 botm (max, min): 1599.3125 1599.3125
Layer 3 top (max, min): 1599.3125 1599.3125
Layer 3 botm (max, min): 1598.8125 1598.8125
Layer 4 top (max, min): 1598.8125 1598.8125
Layer 4 botm (max, min): 1598.3125 1598.3125
Layer 5 top (max, min): 1598.3125 1598.3125
Layer 5 botm (max, min): 1597.8125 1597.8125
Layer 6 top (max, min): 1597.8125 1597.8125
Layer 6 botm (max, min): 1597.3125 1597.3125
Layer 7 top (max, min): 1597.3125 1597.3125
Layer 7 botm (max, min): 1596.8125 1596.8125
Layer 8 top (max, min): 1596.8125 1596.8125
Layer 8 botm (max, min): 1596.3125 1596.3125
Layer 9 top (max, min): 1596.3125 1596.3125
Layer 9 botm (max, min): 1595.8125 1595.8125
Layer 10 top (max, min): 1595.8125 1595.8125
Layer 10 botm (max, min): 1595.3125 1595.3125
Layer 11 top (max, min): 1595.3125 1595.3125
Layer 11 botm (max, min): 1594.8125 1594.8125
Layer 12 top (max, min): 1594.8125 1594.8125
Layer 12 botm (max, min): 1594.3125 1594.3125
Layer 13 top (max, min): 1594.3125 1594.3125
Layer 13 botm (max, min): 1593.8125 1593.8125
Layer 14 top (max, min): 1593.8125 1593.8125
Layer 14 botm (max, min): 1593.3125 1593.3125
Layer 15 top (max, min): 1593.3125 1593.3125
Layer 15 botm (max, min): 1592.8125 1592.8125
Layer 16 top (max, min): 1592.8125 1592.8125
Layer 16 botm (max, min): 1592.3125 1592.3125
Layer 17 top (max, min): 1592.3125 1592.3125
Layer 17 botm (max, min): 1591.8125 1591.8125
Layer 18 top (max, min): 1591.8125 1591.8125
Layer 18 botm (max, min): 1591.3125 1591.3125
Layer 19 top (max, min): 1591.3125 1591.3125
Layer 19 botm (max, min): 1590.8125 1590.8125
Layer 20 top (max, min): 1590.8125 1590.8125
Layer 20 botm (max, min): 1590.3125 1590.3125
Layer 21 top (max, min): 1590.3125 1590.3125
Layer 21 botm (max, min): 1589.8125 1589.8125
Layer 22 top (max, min): 1589.8125 1589.8125
Layer 22 botm (max, min): 1589.3125 1589.3125
Layer 23 top (max, min): 1589.3125 1589.3125
Layer 23 botm (max, min): 1588.8125 1588.8125
Layer 24 top (max, min): 1588.8125 1588.8125
Layer 24 botm (max, min): 1588.3125 1588.3125
Layer 25 top (max, min): 1588.3125 1588.3125
Layer 25 botm (max, min): 1587.8125 1587.8125
Layer 26 top (max, min): 1587.8125 1587.8125
Layer 26 botm (max, min): 1587.3125 1587.3125
Layer 27 top (max, min): 1587.3125 1587.3125
Layer 27 botm (max, min): 1586.8125 1586.8125
Layer 28 top (max, min): 1586.8125 1586.8125
Layer 28 botm (max, min): 1586.3125 1586.3125
Layer 29 top (max, min): 1586.3125 1586.3125
Layer 29 botm (max, min): 1585.8125 1585.8125
Layer 30 top (max, min): 1585.8125 1585.8125
Layer 30 botm (max, min): 1585.3125 1585.3125
Layer 31 top (max, min): 1585.3125 1585.3125
Layer 31 botm (max, min): 1584.8125 1584.8125
Layer 32 top (max, min): 1584.8125 1584.8125
Layer 32 botm (max, min): 1584.3125 1584.3125
Layer 33 top (max, min): 1584.3125 1584.3125
Layer 33 botm (max, min): 1583.8125 1583.8125
Layer 34 top (max, min): 1583.8125 1583.8125
Layer 34 botm (max, min): 1583.3125 1583.3125
Layer 35 top (max, min): 1583.3125 1583.3125
Layer 35 botm (max, min): 1582.8125 1582.8125
Layer 36 top (max, min): 1582.8125 1582.8125
Layer 36 botm (max, min): 1582.3125 1582.3125
Layer 37 top (max, min): 1582.3125 1582.3125
Layer 37 botm (max, min): 1581.8125 1581.8125
Layer 38 top (max, min): 1581.8125 1581.8125
Layer 38 botm (max, min): 1581.3125 1581.3125
Layer 39 top (max, min): 1581.3125 1581.3125
Layer 39 botm (max, min): 1580.8125 1580.8125
Layer 40 top (max, min): 1580.8125 1580.8125
Layer 40 botm (max, min): 1580.3125 1580.3125
Saving Domain Variables#
All relevant variables related to the model domain—including elevation arrays, grid geometry, and layer boundaries—are stored for use in subsequent modeling steps.
Show code cell source
# Save Variables
%store terrain_elevation
%store raster_transform
%store raster_crs
%store raster_bounds_box
%store bed_elevation
%store raster_width
%store raster_height
%store ncol
%store nrow
%store top
%store nlay
%store grid_x
%store grid_y
%store grid_points
%store intersecting_points
%store xorigin
%store yorigin
%store tops
%store botm
Stored 'terrain_elevation' (MaskedArray)
Stored 'raster_transform' (Affine)
Stored 'raster_crs' (CRS)
Stored 'raster_bounds_box' (Polygon)
Stored 'bed_elevation' (float32)
Stored 'raster_width' (float)
Stored 'raster_height' (float)
Stored 'ncol' (int)
Stored 'nrow' (int)
Stored 'top' (MaskedArray)
Stored 'nlay' (int)
Stored 'grid_x' (ndarray)
Stored 'grid_y' (ndarray)
Stored 'grid_points' (GeoDataFrame)
Stored 'intersecting_points' (GeoDataFrame)
Stored 'xorigin' (float)
Stored 'yorigin' (float)
Stored 'tops' (list)
Stored 'botm' (list)