Preprocessing Data

Preprocessing Data#

This automated script is designed to take surface water outputs from HEC-RAS. This section shows the necessary preprocessing steps for the surface water data before it can be utilized in the modeling packages of ModFlow6 with FloPy.

Hide 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 papermill as pm

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

#---------------------- Directories ----------------------#
# Data Paths
# Can be updated using the batch file "run_build_models.bat" or manually
water_surface_elevation_raster = None
terrain_elevation_raster = None
ground_water_domain_shapefile = None
left_boundary_floodplain = None
right_boundary_floodplain = None

## Projection File Path
projection_file = None

### NOTE: This script can be adapted to run groundwater modeling simulations for other locations by changing the directories above to designated inputs

# Parameters
water_surface_elevation_raster = "C:\\Users\\u4eeevmq\\Documents\\Python\\HyporheicFloPy\\CH00365\\GMS\\WSE (Max).USGS_1m_FebApril2018_DEM.U_USGS-3DEP_dtm_hyrdo_flattened_20240412.2041_1.tif"
terrain_elevation_raster = "C:\\Users\\u4eeevmq\\Documents\\Python\\HyporheicFloPy\\CH00365\\GMS\\BlendedTerrain.LoweredTerrain.tif"
ground_water_domain_shapefile = "C:\\Users\\u4eeevmq\\Documents\\Python\\HyporheicFloPy\\CH00365\\InputShapefiles\\GWDomain.shp"
left_boundary_floodplain = "C:\\Users\\u4eeevmq\\Documents\\Python\\HyporheicFloPy\\CH00365\\InputShapefiles\\L_FPL.shp"
right_boundary_floodplain = "C:\\Users\\u4eeevmq\\Documents\\Python\\HyporheicFloPy\\CH00365\\InputShapefiles\\R_FPL.shp"
projection_file = "C:\\Users\\u4eeevmq\\Documents\\Python\\HyporheicFloPy\\CH00365\\RAS\\GIS_Data\\102739_TX_central.prj"

# Print the parameters to verify
print(f"Water Surface Elevation Raster: {water_surface_elevation_raster}")
print(f"Terrain Elevation Raster: {terrain_elevation_raster}")
print(f"Groundwater Domain Shapefile: {ground_water_domain_shapefile}")
print(f"Left Boundary Floodplain: {left_boundary_floodplain}")
print(f"Right Boundary Floodplain: {right_boundary_floodplain}")
print(f"Projection File: {projection_file}")

#--------------------------------- Load Projection FIle -----------------------------------------#
## HEC-RAS Projection File
hec_ras_crs = CRS.from_string(open(projection_file).read().strip())

## Store the CRS
%store hec_ras_crs
### NOTE: If the projection file is not available, the CRS can be defined manually here.

Water Surface Elevation Raster: C:\Users\u4eeevmq\Documents\Python\HyporheicFloPy\CH00365\GMS\WSE (Max).USGS_1m_FebApril2018_DEM.U_USGS-3DEP_dtm_hyrdo_flattened_20240412.2041_1.tif
Terrain Elevation Raster: C:\Users\u4eeevmq\Documents\Python\HyporheicFloPy\CH00365\GMS\BlendedTerrain.LoweredTerrain.tif
Groundwater Domain Shapefile: C:\Users\u4eeevmq\Documents\Python\HyporheicFloPy\CH00365\InputShapefiles\GWDomain.shp
Left Boundary Floodplain: C:\Users\u4eeevmq\Documents\Python\HyporheicFloPy\CH00365\InputShapefiles\L_FPL.shp
Right Boundary Floodplain: C:\Users\u4eeevmq\Documents\Python\HyporheicFloPy\CH00365\InputShapefiles\R_FPL.shp
Projection File: C:\Users\u4eeevmq\Documents\Python\HyporheicFloPy\CH00365\RAS\GIS_Data\102739_TX_central.prj
Stored 'hec_ras_crs' (CRS)

Terrain Elevation Processing#

This section of the workflow handles the processing and reprojection of the terrain elevation raster. The terrain raster is opened and read using rasterio. If the raster’s coordinate reference system (CRS) does not match the target CRS (e.g., the HEC-RAS CRS), it is reprojected so that all spatial data aligns correctly for subsequent modeling steps.

Steps Performed#

  1. Open the Terrain Raster
    The input raster specified by terrain_elevation_raster is opened and read into memory.

  2. Extract Raster Information
    The script retrieves raster data, transformation matrix, and current CRS.

  3. Calculate New Transform and Dimensions
    The script uses calculate_default_transform to determine the transformation matrix and dimensions for the target CRS.

  4. Update Raster Metadata
    Metadata is updated to reflect the new CRS, transform, width, and height.

  5. Reproject and Save
    The raster is reprojected to the target CRS (e.g., hec_ras_crs) and saved as a new file (reprojected_terrain_raster.tif).

  6. Output
    The reprojected raster is saved and ready for use in subsequent modeling or analysis.

#------------------------------Model Grid--------------------------------#
# Open the terrain raster
with rasterio.open(terrain_elevation_raster) as src:
    terrain_elevation = src.read(1)
    raster_transform = src.transform
    raster_crs = src.crs

    # Calculate the new transform and dimensions
    transform, width, height = calculate_default_transform(
        raster_crs, hec_ras_crs, src.width, src.height, *src.bounds
    )

    # Define metadata for the new raster
    new_meta = src.meta.copy()
    new_meta.update({
        "crs": hec_ras_crs,
        "transform": transform,
        "width": width,
        "height": height
    })

    # Define output file name
    terrain_output_raster = "reprojected_terrain_raster.tif"

    # Reproject and save
    with rasterio.open(terrain_output_raster, "w", **new_meta) as dst:
        reproject(
            source=rasterio.band(src, 1),
            destination=rasterio.band(dst, 1),
            src_transform=raster_transform,
            src_crs=raster_crs,
            dst_transform=transform,
            dst_crs=hec_ras_crs,
            resampling=Resampling.nearest
        )

print(f"Reprojected raster saved as {terrain_output_raster}")

Reprojected raster saved as reprojected_terrain_raster.tif
Hide code cell source 
## Store the terrain elevation raster
%store terrain_elevation
%store raster_transform
%store transform
%store raster_crs
%store terrain_output_raster

Stored 'terrain_elevation' (ndarray)
Stored 'raster_transform' (Affine)
Stored 'transform' (Affine)
Stored 'raster_crs' (CRS)
Stored 'terrain_output_raster' (str)

Water Serface Elevation Processing#

This section of the workflow handles the processing and reprojection of the water surface elevation raster. The raster is opened and read using rasterio. If the raster’s coordinate reference system (CRS) does not match the target CRS (e.g., the HEC-RAS CRS), it is reprojected so that all spatial data aligns correctly for subsequent modeling steps.

Steps Performed#

  1. Open the Water Serface Elevation Raster
    The input raster specified by water_surface_elevation_raster is opened and read into memory.

  2. Extract Raster Information
    The script retrieves raster data, transformation matrix, and current CRS.

  3. Calculate New Transform and Dimensions
    The script uses calculate_default_transform to determine the transformation matrix and dimensions for the target CRS.

  4. Update Raster Metadata
    Metadata is updated to reflect the new CRS, transform, width, and height.

  5. Reproject and Save
    The raster is reprojected to the target CRS (e.g., hec_ras_crs) and saved as a new file (reprojected_water_surface_elevation_raster.tif). The raster is cropped to the extent of the terrain (e.g., terrain_output_raster) and saved as a new file (cropped_water_surface_elevation_raster.tif).

  6. Output
    The cropped and reprojected raster is saved and ready for use in subsequent modeling or analysis.

#--------------------------------- Water Surface Elevation -----------------------------------------#
# Open the input raster
with rasterio.open(water_surface_elevation_raster) as src:
    surface_elevation = src.read(1)
    raster_transform = src.transform
    raster_crs = src.crs

    # Calculate the new transform and dimensions
    transform, width, height = calculate_default_transform(
        raster_crs, hec_ras_crs, src.width, src.height, *src.bounds
    )
    
    # Define metadata for the new raster
    new_meta = src.meta.copy()
    new_meta.update({
        "crs": hec_ras_crs,
        "transform": transform,
        "width": width,
        "height": height
    })

    # Define output file name
    output_raster = "reprojected_water_surface_elevation_raster.tif"

    # Reproject and save
    with rasterio.open(output_raster, "w", **new_meta) as dst:
        reproject(
            source=rasterio.band(src, 1),
            destination=rasterio.band(dst, 1),
            src_transform=raster_transform,
            src_crs=raster_crs,
            dst_transform=transform,
            dst_crs=hec_ras_crs,
            resampling=Resampling.nearest
        )

print(f"Reprojected raster saved as {output_raster}")

# Crop the water surface elevation raster to the terrain extent
with rasterio.open(output_raster) as src:
    terrain_bounds = src.bounds
    terrain_geom = box(*terrain_bounds)

    with rasterio.open(water_surface_elevation_raster) as src:
        out_image, out_transform = mask(src, [terrain_geom], crop=True)
        out_meta = src.meta.copy()
        out_meta.update({
            "driver": "GTiff",
            "height": out_image.shape[1],
            "width": out_image.shape[2],
            "transform": out_transform
        })

        water_surface__output_raster = "cropped_water_surface_elevation_raster.tif"
        with rasterio.open(water_surface__output_raster, "w", **out_meta) as dst:
            dst.write(out_image)

print(f"Cropped water surface elevation raster saved as {water_surface__output_raster}")

Reprojected raster saved as reprojected_water_surface_elevation_raster.tif
Cropped water surface elevation raster saved as cropped_water_surface_elevation_raster.tif

Floodplain and Groundwater Domain Processing#

This section of the workflow handles the processing and reprojection of shapefiles that define the groundwater domain and floodplain boundaries. Accurate reprojection ensures all spatial data layers are aligned in the same coordinate reference system (CRS) for reliable groundwater modeling.

Steps Performed#

  1. Load Shapefiles
    The script reads in the groundwater domain and floodplain boundary shapefiles using geopandas.read_file().

  2. Reproject to Target CRS
    Each shapefile is reprojected to the target CRS used for modeling (e.g., hec_ras_crs) using the .to_crs() method. This ensures all spatial inputs are consistent and compatible.

    • ground_water_domain = gpd.read_file(ground_water_domain_shapefile).to_crs(hec_ras_crs)

    • left_boundary = gpd.read_file(left_boundary_floodplain).to_crs(hec_ras_crs)

    • right_boundary = gpd.read_file(right_boundary_floodplain).to_crs(hec_ras_crs)

  3. Output
    The reprojected shapefiles are ready for use in subsequent modeling, ensuring that domain and boundary definitions match the model grid and other raster data.

#--------------------------------- Shapefiles -----------------------------------------#
## Load and reproject shapefiles to match the raster CRS
ground_water_domain = gpd.read_file(ground_water_domain_shapefile).to_crs(hec_ras_crs)
left_boundary = gpd.read_file(left_boundary_floodplain).to_crs(hec_ras_crs)
right_boundary = gpd.read_file(right_boundary_floodplain).to_crs(hec_ras_crs)
Hide code cell source 
## Save the reprojected files
%store water_surface__output_raster
%store ground_water_domain
%store left_boundary
%store right_boundary

Stored 'water_surface__output_raster' (str)
Stored 'ground_water_domain' (GeoDataFrame)
Stored 'left_boundary' (GeoDataFrame)
Stored 'right_boundary' (GeoDataFrame)
Contents