diff --git a/notebooks/python/graphflood/README.md b/notebooks/python/graphflood/README.md new file mode 100644 index 0000000..bfc954d --- /dev/null +++ b/notebooks/python/graphflood/README.md @@ -0,0 +1,9 @@ +# `graphflood` examples and demos + +You will find in this directory the demos and tutorial for using [graphflood](https://esurf.copernicus.org/articles/12/1295/2024/) in `pytopotoolbox`. + +- `graphflood.ipynb`: the example from [topotoolbox 3 paper](https://egusphere.copernicus.org/preprints/2026/egusphere-2026-2478/). Run graphflood for multiple discharges, compute width, swath, plot... Mostly to reproduce the paper figures and analysis. +- `Workshop_06_2026.ipynb`: a complete tutorial presented during the EGU ECR webinar of June 2026 - probably the best starting point so far (will add link when the recording is available on youtube!) + + + diff --git a/notebooks/python/graphflood/Workshop_06_2026.ipynb b/notebooks/python/graphflood/Workshop_06_2026.ipynb new file mode 100644 index 0000000..f0e2e1b --- /dev/null +++ b/notebooks/python/graphflood/Workshop_06_2026.ipynb @@ -0,0 +1,963 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "556e3d82-b6d3-4238-a984-8d6ecf2d253a", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01a20902-328a-44be-bf7a-dadea6097c63", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [], + "source": [ + "# Ghost cell to set font style for the presentation\n", + "from IPython.display import HTML\n", + "\n", + "HTML(\"\"\"\n", + "\n", + "\"\"\")" + ] + }, + { + "cell_type": "markdown", + "id": "272dd83c-d6f6-4901-abde-5a7dd7854406", + "metadata": {}, + "source": [ + "
\n", + "\n", + "# 🌊 Graphflood in `pytopotoolbox`\n", + "\n", + "`graphflood` (Gailleton et al., 2024) simulates how water moves across terrain surfaces using a fast 2D stationary diffusive wave approximation.\n", + "\n", + "Give it an elevation map (DEM) and some water (rainfall fields or point discharge paths), and it routes flow until it outputs balanced, equilibrated 2D parameter maps:\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Dynamic MetricVariable
πŸ”Ή Flow Depth$h$
πŸ”Ή Velocity$u$
πŸ”Ή Shear Stress$\\tau$
πŸ”Ή Discharge$Q$
\n", + "\n", + "---\n", + "\n", + "### πŸ“¦ Prerequisites\n", + "\n", + "Install dependencies via your terminal:\n", + "```bash\n", + "pip install topotoolbox\n", + "pip install ipympl # Optional: For interactive figures\n", + "```\n", + "\n", + "\n", + "### Offline assistance?\n", + "\n", + "- [Github Discussion](https://github.com/orgs/TopoToolbox/discussions)\n", + "- Last resort: grop me an email (boris.gailleton@univ-rennes.fr)\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e0c8bce1-f57b-47ea-bff9-807ef99bc117", + "metadata": {}, + "outputs": [], + "source": [ + "# General imports\n", + "import topotoolbox as ttb\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import time\n", + "import scipy.stats as st\n", + "from scipy.ndimage import binary_closing, distance_transform_edt\n", + "from skimage.morphology import skeletonize\n", + "\n", + "# Set matplotlib backend\n", + "# %matplotlib inline\n", + "## If you want interactive plots\n", + "%matplotlib ipympl\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f82d1868-9643-4118-931a-6b782d3283bb", + "metadata": { + "jupyter": { + "source_hidden": true + } + }, + "outputs": [], + "source": [ + "# --- Webinar & High-Contrast Readable Settings ---\n", + "plt.rcParams['figure.figsize'] = (10, 6) # Large default size for screen sharing\n", + "plt.rcParams['figure.dpi'] = 110 # Crisp rendering on high-res streams\n", + "plt.rcParams['font.size'] = 14 # Large text for readability\n", + "plt.rcParams['axes.labelsize'] = 16\n", + "plt.rcParams['axes.titlesize'] = 18\n", + "plt.rcParams['xtick.labelsize'] = 12\n", + "plt.rcParams['ytick.labelsize'] = 12\n", + "\n", + "# --- Monospace Font Strategy ---\n", + "plt.rcParams['font.family'] = 'monospace'\n", + "plt.rcParams['font.monospace'] = ['Fira Code', 'JetBrains Mono', 'Consolas', 'Courier New', 'monospace']\n", + "\n", + "# --- Cyberpunk / Geek Dark Mode Theme ---\n", + "plt.rcParams['figure.facecolor'] = '#121214' # Deep charcoal background\n", + "plt.rcParams['axes.facecolor'] = '#1a1a1e' # Slightly lighter plot area\n", + "plt.rcParams['axes.edgecolor'] = '#444444' # Subtle borders\n", + "plt.rcParams['axes.linewidth'] = 1.5\n", + "\n", + "# High-contrast text colors\n", + "plt.rcParams['text.color'] = '#e0e0e0'\n", + "plt.rcParams['axes.labelcolor'] = '#e0e0e0'\n", + "plt.rcParams['xtick.color'] = '#a0a0a0'\n", + "plt.rcParams['ytick.color'] = '#a0a0a0'\n", + "\n", + "# Grid settings\n", + "plt.rcParams['axes.grid'] = True\n", + "plt.rcParams['grid.color'] = '#2d2d34'\n", + "plt.rcParams['grid.linestyle'] = '--'\n", + "plt.rcParams['grid.linewidth'] = 0.8\n", + "plt.rcParams['grid.alpha'] = 0.4\n", + "\n", + "# --- Neon / Synthwave Color Palette ---\n", + "plt.rcParams['axes.prop_cycle'] = plt.cycler(color=[\n", + " '#00ffcc', '#ff007f', '#9d4edd', '#3a86ff', '#ffb703', '#70e000'\n", + "])\n", + "\n", + "# Image defaults\n", + "plt.rcParams['image.cmap'] = 'magma' # Dark-mode default colormap\n" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c7f1587a-9f73-4b47-bbe6-48107d8578ec", + "metadata": {}, + "source": [ + "
\n", + "\n", + "# ⚑ Quickstart\n", + "\n", + "Provide a raw DEM and your water inputs to calculate balanced flow conditions and track fluid behavior maps instantly.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5871f3b8-82a7-44f2-ac37-39043f5b46f3", + "metadata": {}, + "outputs": [], + "source": [ + "# Load DEM and crop to a specific spatial subset via percentage thresholds\n", + "dem = ttb.load_dem('greenriver')\n", + "dem = dem.crop(0.35, 0.68, 0.3, 0.85, 'percent')\n", + "\n", + "# if you want to load your DEM:\n", + "# dem = ttb.read_tif('mydem.tif')\n", + "\n", + "# Plot the cropped topography with an overlaid hillshade\n", + "fig, ax = plt.subplots()\n", + "cb = ax.imshow(dem.z, cmap='terrain', extent=dem.extent)\n", + "ax.imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.6)\n", + "plt.colorbar(cb, label='elevation (m)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "c0887efe-039e-4bb9-a623-30502aff46a2", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## πŸ› οΈ `GFObject` (Graphflood Object)\n", + "\n", + "The primary interface utilized to initialize and run the `graphflood` engine.\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
ArgumentRequirementDescription
❗ gridRequiredCore topography/DEM input grid layout.
βš™οΈ bcsOptionalPixel boundary condition codes. Defaults to unconstrained borders.
βš™οΈ pOptionalPrecipitation supply rates (scalar value or 2D array matrix in $m/s$).
βš™οΈ manningOptionalManning surface roughness index ($n$) used during flow calculations.
\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f197b47-6c5b-487a-a41e-9eee4faba57e", + "metadata": {}, + "outputs": [], + "source": [ + "# Instantiate the Graphflood object with 100 mm/h uniform precipitation converted to m/s\n", + "## 100 mm/h is super high, I chose it to make sure even in our v. small demo site we generate enough flooding\n", + "gfo = ttb.GFObject(dem, p=(100e-3)/3600, manning=0.05)" + ] + }, + { + "cell_type": "markdown", + "id": "400f33ad-fdc7-42f1-824c-d916f8e3801b", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## πŸ”„ Running the Model\n", + "\n", + "The solver routes surface water incrementally across cells until it settles into equilibrium.\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Control ParameterImpact on ConvergenceConstraint / Risk
πŸ”’ IterationsHigher execution runs bring the solver closer to true mathematical flow alignment.Increases computational compute time.
⏱️ Time Step (dt)Larger steps speed up convergence rates towards a steady-state layout.Excessive values trigger numerical instability and cause solution degradation.
\n", + "\n", + "
\n", + " ⚠️ Key Distinction: The parameter dt functions purely as a mathematical tuning knob to speed up calculations. It does not reflect actual physical timeline tracking.\n", + "
\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a6debe4-8763-4f30-8312-1197359935d5", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the steady-state solver loop for 100 fixed steps\n", + "gfo.run_n_iterations(dt=1e-2, n_iterations=100)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "beffa732-3bf6-4907-8792-91c91e270fd5", + "metadata": {}, + "outputs": [], + "source": [ + "# Display calculated water depths layered over the relief hillshade\n", + "fig, ax = plt.subplots()\n", + "cb = ax.imshow(gfo.hw, cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "ax.imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.4)\n", + "plt.colorbar(cb, label='Flow depth (m)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "22cf6330-b203-4194-9efa-0a7c284eec46", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## πŸ“Š Computed Metrics\n", + "\n", + "### πŸ’₯ Shear Stress ($\\tau$)\n", + "\n", + "The total dragging friction force exerted by shifting water directly along the surface channel bed.\n", + "\n", + "$$\\tau = \\rho g h S_w$$\n", + "\n", + "
\n", + " βš™οΈ Engine Execution: Extracted directly using gfo.compute_tau() based on calculated local depth ($h$) and hydraulic water surface slope profiles ($S_w$).\n", + "
\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4cc2eb88-541e-40cb-99d6-1cac9ef22c42", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute and map boundary shear stress across the channel bed\n", + "shear_stress = gfo.compute_tau()\n", + "\n", + "# \n", + "fig, ax = plt.subplots()\n", + "cb = ax.imshow(shear_stress, cmap='magma', extent=dem.extent, vmax=100)\n", + "ax.imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.4)\n", + "plt.colorbar(cb, label=r'Shear Stress ($\\tau$)')\n", + "plt.show()\n" + ] + }, + { + "cell_type": "markdown", + "id": "5b6bac0c-329d-481e-a2da-38ebf31555b3", + "metadata": {}, + "source": [ + "
\n", + "\n", + "### πŸ† Flow Velocity\n", + "\n", + "Calculates flow velocity magnitudes ($u$) across your target grid landscape layout.\n", + "\n", + "
\n", + " βš™οΈ Execution: gfo.get_u()\n", + "
πŸš€ Roadmap: Directional vector component tracking layouts will follow in an upcoming update.\n", + "
\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8315c146-d941-4bbc-ba06-14b6f15d3a4b", + "metadata": {}, + "outputs": [], + "source": [ + "# Retrieve depth-averaged flow velocity magnitudes\n", + "flow_velocity = gfo.get_u()\n", + "\n", + "fig, ax = plt.subplots()\n", + "cb = ax.imshow(flow_velocity, cmap='viridis', extent=dem.extent, vmax = 2.)\n", + "ax.imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.4)\n", + "plt.colorbar(cb, label=r'$u$ flow velocity $\\frac{m}{s}$')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "ee2215e8-0674-47fa-9f23-f6eed2cc1fc7", + "metadata": {}, + "source": [ + "
\n", + "\n", + "### πŸ“‰ Discharge & Hydraulic Slope\n", + "\n", + "A custom hydraulic twist on classical geomorphic **slope-area plotting**, replacing upstream drainage basin surface area ($A$) with explicitly calculated flow discharge tracking values ($Q$).\n", + "\n", + "
\n", + " πŸ” Core Objective: Compares log-scaled surface slope records ($S_w$) against discharge rates ($Q$) to distinguish distinct flow channels and landscape scaling domains.\n", + "
\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8757e904-7fca-4254-91be-8ca6e6cb643d", + "metadata": {}, + "outputs": [], + "source": [ + "# Flatten volumetric discharge and hydraulic water surface slope grids into log space\n", + "Q_in = np.log10(gfo.get_qvol_i().z.ravel())\n", + "Sw = np.log10(gfo.get_sw().z.ravel())\n", + "\n", + "# Filter out NaN and infinite elements\n", + "mask = np.isfinite(Q_in) & np.isfinite(Sw)\n", + "Q_in = Q_in[mask]\n", + "Sw = Sw[mask]\n", + "\n", + "# 1. Compute 2D Density (Binned Counts)\n", + "bins2D = 60\n", + "stat_2d, x_edges, y_edges, _ = st.binned_statistic_2d(Q_in, Sw, None, 'count', bins=bins2D)\n", + "# \n", + "# 2. Compute 1D Profile (Median of Y binned by X)\n", + "bins = 30\n", + "bin_means, bin_edges, _ = st.binned_statistic(Q_in, Sw, statistic='median', bins=bins)\n", + "bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2\n", + "\n", + "# 3. Plot\n", + "fig, ax = plt.subplots()\n", + "\n", + "# 2D Density Map distribution\n", + "ax.imshow(stat_2d.T, origin='lower', cmap='magma', aspect='auto',\n", + " extent=[x_edges[0], x_edges[-1], y_edges[0], y_edges[-1]])\n", + "\n", + "# 1D Overlaid Median Line trend line\n", + "ax.scatter(bin_centers, bin_means, facecolor='cyan', edgecolor='gray', lw=3, s=200)\n", + "\n", + "ax.set_ylim(-2, -0.5)\n", + "ax.set_xlabel('log Q')\n", + "ax.set_ylabel('log Sw')\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bef654ff-afbd-4604-af4b-f146bf122695", + "metadata": {}, + "outputs": [], + "source": [ + "# Initialize empty mask grid matching the topography geometry\n", + "domain = np.zeros_like(dem.z)\n", + "# \n", + "text_arrays = \"\"\"\n", + "Extract raw arrays for discharge and water surface slopes\n", + "\"\"\"\n", + "Qi2d = np.log10(gfo.get_qvol_i().z)\n", + "Sw2d = np.log10(gfo.get_sw().z)\n", + "# \n", + "# Isolate channel cells using joint flow thresholds\n", + "domain[(Qi2d > -1.7) & (Sw2d < -1.4) & (gfo.hw.z > 0.05)] = 1\n", + "# Clean up gaps and edge noise via morphological binary closing operations\n", + "domain = binary_closing(domain, iterations=3)\n", + "# \n", + "# Visualize the calculated binary river network mask\n", + "fig, ax = plt.subplots()\n", + "cb = ax.imshow(domain, cmap='cividis', extent=dem.extent)\n", + "ax.imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "34fcaa2c-a077-461b-a080-1619c5e3b695", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## πŸ“ Automated Channel Width Extraction\n", + "\n", + "A clean, efficient post-processing workflow executed entirely via standard `scipy` and `scikit-image` sequences to map true channel width paths straight from our 2D hydraulic footprint.\n", + "\n", + "
\n", + " πŸ”— Reference Code: Review the complete original script logic inside the official TopoToolbox Paper Gallery Notebook.\n", + "
\n", + "\n", + "### πŸ”„ Processing Pipeline\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
StepOperationCore Target Function
1Clean Channel MaskIsolate active river zones (based on slope-discharge limits) and remove interior holes using scipy.ndimage.binary_closing.
2Extract Centerline SpinesThin the continuous 2D mask tracks down into an individual 1-pixel wide skeleton network via skimage.morphology.skeletonize.
3Calculate Edge DistanceCompute exact Euclidean coordinate distances from centerline nodes to the nearest boundary edge using scipy.ndimage.distance_transform_edt.
\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "916bc996-43dd-4a52-9e55-095e16724ddd", + "metadata": {}, + "outputs": [], + "source": [ + "# =============================================================================\n", + "# 1. MORPHOLOGICAL ANALYSIS & DISTANCE TRANSFORM\n", + "# =============================================================================\n", + "\n", + "# Map the exact Euclidean distance from every channel pixel to the nearest dry bank boundary\n", + "dt_edt = distance_transform_edt(domain)\n", + "\n", + "# Thin the 2D binary channel mask down to a continuous, 1-pixel wide centerline network\n", + "skel = skeletonize(domain)\n", + "\n", + "\n", + "# =============================================================================\n", + "# 2. CHANNEL WIDTH COMPUTATION & FILTERING\n", + "# =============================================================================\n", + "\n", + "# Calculate absolute physical channel width (m) along the centerline skeleton\n", + "# Formula tracks full-width (2 * half-width) and adjusts for discrete pixel resolution\n", + "width = 2 * dt_edt[skel] * dem.cellsize - dem.cellsize\n", + "\n", + "# Generate a boolean filtering mask to isolate features >= 2 meters wide\n", + "mask = width >= 2\n", + "\n", + "\n", + "# =============================================================================\n", + "# 3. SPATIAL COORDINATE TRANSFORMATION\n", + "# =============================================================================\n", + "\n", + "# Extract the 2D array matrix row and column indices for all skeleton nodes\n", + "rows_skel, cols_skel = np.where(skel)\n", + "\n", + "# Project pixel index positions into true geographic map coordinates (X, Y)\n", + "xskel, yskel = ttb.transform_coords(\n", + " dem, \n", + " rows_skel, \n", + " cols_skel,\n", + " input_mode='indices2D', \n", + " output_mode='coordinates'\n", + ")\n", + "\n", + "\n", + "# =============================================================================\n", + "# 4. GEOGRAPHIC VISUALIZATION\n", + "# =============================================================================\n", + "\n", + "fig, ax = plt.subplots()\n", + "\n", + "# Base Layer: Render the topography hillshade in grayscale\n", + "ax.imshow(dem.hillshade(), cmap='gray', extent=dem.extent, vmax=1.2)\n", + "ax.imshow(domain, cmap='cividis', extent=dem.extent, alpha=0.25)\n", + "\n", + "# Data Overlay: Plot the filtered centerline coordinates colored by calculated width\n", + "cb = ax.scatter(\n", + " xskel[mask], \n", + " yskel[mask], \n", + " lw=0, \n", + " c=width[mask], \n", + " cmap='plasma'\n", + ")\n", + "\n", + "# Render colorbar interface and output graphic area\n", + "plt.colorbar(cb, label='width (m)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "fbda712b-b850-42cf-8bd7-e0dbbe78a7e0", + "metadata": {}, + "source": [ + "# 🏞️🏞️🏞️🏞️🏞️ Done with the Quickstart 🏞️🏞️🏞️🏞️🏞️\n", + "\n", + "
\n", + "
\n", + " ⛰️ 🌲 🌲 🌊 🌲 ⛰️ 🌲 🌲 🌲 ⛰️ 🌲 🌲 🌊 🌲 ⛰️ \n", + "
\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "a62df845-7666-4521-8daf-f0e1a12ca93a", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## 🏁 Checking If the Model Is Done (Convergence)\n", + "\n", + "How can we verify when the simulation has finished its job? We track our system metrics to find where local water levels stop changing and lock down into a stable state.\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
What We CheckIdeal GoalWhat Actually Happens
πŸ“Š Water BalanceWater In = Water Out
($Q_{in} = Q_{out}$)		The total volume exiting map parameters will slowly approach absolute parity with inflow records over time.
🌊 Stable Water DepthDepth stops changing
($\\frac{dh}{dt} = 0$)		While water depths rarely freeze completely to absolute zero change, grid values level off enough to stop mattering very quickly.
\n", + "\n", + "---\n", + "\n", + "### βš™οΈ Quick Troubleshooting Strategy\n", + "\n", + "* **Run it longer:** If water surfaces are actively moving, extend processing loops using standard calls to `gfo.run_n_iterations(...)`.\n", + "* **Best way to check:** Monitor the structural changes in depth profile maps between cycles ($dh$). Once variance adjustments shrink past a nominal target baseline, your model is balanced.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af063efa-0aca-4282-a9d5-61b76313a0f2", + "metadata": {}, + "outputs": [], + "source": [ + "# Re-instantiate model instance to benchmark step-wise iteration increments\n", + "gfo = ttb.GFObject(dem, p=100e-3/3600, manning=0.05)\n", + "dhs = [gfo.hw.z.copy()]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd2f656c-4339-44f0-a915-61e9e573e09e", + "metadata": {}, + "outputs": [], + "source": [ + "# \n", + "N = 10\n", + "N_it = 10\n", + "# Stepwise solver tracking loop to monitor water depth grid adjustments over execution\n", + "for i in range(N):\n", + " print(f'running {(i+1) * N_it}/{N * N_it}', end='\\r', flush=True)\n", + " gfo.run_n_iterations(dt=5e-2, n_iterations=N_it)\n", + " # Enforce boundary zero-depth conditions at perimeter edges\n", + " gfo.hw[[0, -1], :] = 0\n", + " gfo.hw[:, [-1, 0]] = 0\n", + " dhs.append(gfo.hw.z.copy())\n", + "dhs_arr = np.array(dhs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dae27cde-2d7c-466f-8f43-b6738bafca80", + "metadata": {}, + "outputs": [], + "source": [ + "# Close any lingering plots and check depth distribution snapshots at steps 1, 2, 5, 10\n", + "# plt.close('all')\n", + "fig, axes = plt.subplots(2, 2)\n", + "cb = axes[0, 0].imshow(dhs_arr[1], cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "cb = axes[0, 1].imshow(dhs_arr[2], cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "cb = axes[1, 0].imshow(dhs_arr[5], cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "cb = axes[1, 1].imshow(dhs_arr[10], cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75f79be6-4e09-46b6-bce0-059979480e0d", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the log-scaled progression of the 90th percentile depth variance delta to assess convergence\n", + "# plt.close('all')\n", + "fig, ax = plt.subplots()\n", + "ax.plot(np.arange(len(dhs)-1)*N_it, np.percentile(np.abs(np.diff(dhs_arr, axis=0)), 90, axis=(1, 2)), lw=3, color='w')\n", + "ax.set_xlabel('N iterations')\n", + "ax.set_ylabel(f'90th perc. of h increment every {N_it} iterations')\n", + "ax.set_yscale('log')\n", + "time.sleep(0.2) # Mitigation step preventing rendering race condition bugs inside JupyterLab interface\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "59c3d1e1-e219-4429-80d6-a10ea7f1d3b8", + "metadata": {}, + "source": [ + "
\n", + "\n", + "## πŸ—ΊοΈ Boundary Conditions: The \"Reach\" Case\n", + "\n", + "Isolate explicit river paths by defining routing directions: driving inbound stream channels through a selected boundary zone (South) and establishing exit points along another (North).\n", + "\n", + "### πŸ”’ Boundary Condition (BC) Grid Codes\n", + "\n", + "Pixel actions are set using an integer array layer that exactly mirrors our topography layout geometry:\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
BC CodeTypeFlow Behavior
0🧱 No Data / WallSolid impermeable boundary. Flow is completely blocked and cannot pass or exit.
1🌊 Normal Active AreaStandard interior routing. Flow moves continuously across cells but cannot scale over borders.
3πŸšͺ Outlet GateUnconstrained drainage point. Flow permanently empties out of the simulation at this node.
\n", + "\n", + "
\n", + " ⚠️ Critical Rule: Active simulation runs require at least one pixel configured as an outlet gate (3). If an exit path is missing, water pools infinitely and the solver cannot balance metrics.\n", + "
\n", + "\n", + "### πŸ“ Spatial Input Mapping (Point Source Discharge)\n", + "\n", + "Instead of applying uniform rainfall across the entire map footprint, we switch to a **spatially variable input array**. This zeroes out fluid additions everywhere except at our designated channel source node, precisely replicating a localized upstream stream inflow.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd94a41c-0221-44ab-bcd4-539973d909a0", + "metadata": {}, + "outputs": [], + "source": [ + "# Set default processing mask states to internal routing (1) \n", + "bcs = np.ones_like(dem.z, dtype=np.uint8)\n", + "# Impose solid impermeable boundaries (0) around the immediate exterior rows and columns\n", + "bcs[[-1, 0], :] = 0\n", + "bcs[:, [-1, 0]] = 0\n", + "\n", + "# Construct empty array map to manage explicit flow source locations\n", + "prec_in = np.zeros_like(dem.z, dtype=np.float32)\n", + "\n", + "# Draw background reference surface map\n", + "fig, ax = plt.subplots()\n", + "ax.imshow(dem.hillshade(), cmap='gray')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "586b4322-c822-44f6-a72b-9698c693728d", + "metadata": {}, + "outputs": [], + "source": [ + "# Create an unconstrained drainage exit point (3) on a segment of the boundary wall\n", + "bcs[20:50, -1] = 3\n", + "\n", + "# Inject a localized discharge point source of 5 m3/s scaled to a precipitation equivalent flux\n", + "prec_in[-2, 30] = 5 / (dem.cellsize**2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36f54b85-f3a0-4c20-8a80-6e32f39ddba3", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# Draw background reference surface map\n", + "fig, ax = plt.subplots()\n", + "ax.imshow(dem.hillshade(), cmap='gray')\n", + "ax.imshow(bcs, cmap='jet', alpha = 0.7)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c599dbeb-067d-4ac9-baa7-af4e79dab8b4", + "metadata": {}, + "outputs": [], + "source": [ + "# Initialize a new Graphflood model configuration using the boundary matrix and localized flow input\n", + "gfo = ttb.GFObject(dem, p=prec_in, manning=0.05, bcs=bcs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e89efaa-268d-4b39-ab5d-5c6c89994268", + "metadata": {}, + "outputs": [], + "source": [ + "# Execute processing loop and mask out non-contributing dry terrain nodes\n", + "gfo.run_n_iterations(dt=1e-2, n_iterations=100)\n", + "gfo.hw.z[gfo.get_qvol_i().z == 0] = 0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0fab88c2-87d0-41a6-9a68-529ce0e26ac7", + "metadata": {}, + "outputs": [], + "source": [ + "# View inundation track profile resulting from the isolated reach boundary settings\n", + "fig, ax = plt.subplots()\n", + "cb = ax.imshow(gfo.hw, cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "ax.imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.4)\n", + "plt.colorbar(cb, label='Flow depth (m)')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a550609a-df4b-4e32-9c6d-6d840f9644c4", + "metadata": {}, + "outputs": [], + "source": [ + "# Archive current water depth output, scale point discharge by 10x, and update flow fields\n", + "thw = gfo.hw.z.copy()\n", + "gfo.precipitations = prec_in * 10\n", + "gfo.run_n_iterations(dt=1e-3, n_iterations=100)\n", + "gfo.hw.z[gfo.get_qvol_i().z == 0] = 0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "05a925d3-bf76-4c82-a9df-d64a7909460b", + "metadata": {}, + "outputs": [], + "source": [ + "# Compare flow depth signatures between the 5 m3/s and 50 m3/s discharge inputs side-by-side\n", + "fig, axes = plt.subplots(1, 2)\n", + "axes[0].imshow(gfo.hw, cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "axes[1].imshow(thw, cmap='Blues', extent=dem.extent, vmax=0.5)\n", + "axes[0].imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.4)\n", + "axes[1].imshow(dem.hillshade(), cmap='gray', extent=dem.extent, alpha=0.4)\n", + "\n", + "axes[0].set_title(r'50 $m^3.s^{-1}$')\n", + "axes[1].set_title(r'5 $m^3.s^{-1}$')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b72d3b81-072b-4268-809c-4c3f8737a837", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "43090742-a4eb-415d-80da-85ab240e1c74", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}