A flood simulation from terrain and rain

This page takes the most common task from start to finish: you have a terrain model, a roughness map and a storm, and you want the peak flood depth. It uses the terrain bundled with the examples, a real 8 x 8 km patch of Cook County, Illinois at 10 m, so every block below runs as written from the repository root.

The whole script

import numpy as np
from geoswe import Mesh2D, Config, Solver2D, RainfallForcing

# 1. terrain and roughness: arrays indexed [x, y], elevations in metres
case = np.load("examples/data/cookcounty_mini.npz")
win = np.s_[336:464, 336:464]                # a 1.3 km window; use np.s_[:, :] on a GPU
bed, manning = case["bed"][win], case["manning"][win]
nx, ny = bed.shape
mesh = Mesh2D(nx=nx, ny=ny, dx=float(case["dx"]), dy=float(case["dy"]))

# 2. a storm: 75 mm/h for one hour, then dry
rain = RainfallForcing(time_s=[0, 3600], rate_mm_h=[75, 0])

# 3. the solver: dry start, water leaves freely at the edges
cfg = Config(friction="manning", bc_x="fall", bc_y="fall", rainfall_forcing=rain)
solver = Solver2D(mesh, cfg, np.zeros((3, nx, ny)), bed)
solver.set_manning(manning)

# 4. run two hours and read the result
solver.run(t_end=2 * 3600)
peak = solver.max_depth()                    # NumPy array (nx, ny), metres
print(f"deepest water {peak.max():.2f} m; {100 * (peak > 0.05).mean():.0f}% of the area reached 5 cm")

It prints deepest water 1.05 m; 27% of the area reached 5 cm, in about half a minute on a CPU and two seconds on a GPU. With the window removed, all 640,000 cells run in about eight seconds on one GPU.

The same script runs on both backends. GeoSWE uses the GPU when CuPy and a CUDA device are present and NumPy otherwise; set GEOSWE_BACKEND=numpy or GEOSWE_BACKEND=cupy to choose.

What each part does

Arrays. Every field you pass has shape (nx, ny) and is indexed [i, j] with i along x and j along y. The state is q = [h, hu, hv] with shape (3, nx, ny): depth and the two discharges per unit width, in SI units. A dry start is an array of zeros. The solver adds its own ghost cells around the arrays; you do not pad anything.

Configuration. Config() with no arguments is the scheme of the GeoSWE paper: first-order HLLC fluxes, the well-balanced surface-reconstruction bed treatment, forward Euler at CFL 0.5. A flood run adds three things to it: friction, rain, and boundaries that let water out. Precision follows the backend (single on the GPU, double on the CPU) unless you pass dtype.

Roughness. friction="manning" switches friction on, and solver.set_manning(n) takes Manning’s n as one number or as an (nx, ny) map.

Rain. RainfallForcing(time_s, rate_mm_h) holds each rate from its time until the next one, and the last rate from then on, so [75, 0] at [0, 3600] is one hour of rain. For rain that varies in space, give one (nx, ny) frame per time: rate_mm_h of shape (nt, nx, ny).

Edges. "fall" lets water leave and nothing enter, which is what a rainfall-runoff domain usually wants. "wall" closes an edge and "extrapolate" (the default) is a zero-gradient open edge. See boundary conditions.

Running. run(t_end) chooses each time step from the CFL condition. On a dry bed that condition says nothing, so while it rains run also keeps the step below the CFL step of the film the rain lays down; you do not have to cap the first steps by hand. run(t_end, dt_max=...) adds your own cap, and run(t_end, callback=f) calls f(solver, step) after every step.

Results. solver.depth() is the depth now and solver.max_depth() the largest depth each cell has held, both as NumPy arrays of shape (nx, ny). solver.q_interior is the full state on the device.

Your own terrain

With the io extra (rasterio), a GeoTIFF becomes an array in the layout above, and a result goes back out with its georeferencing:

from geoswe.io_geotiff import read_geotiff, write_geotiff, GeoArray

dem = read_geotiff("dem.tif")                # dem.data is (nx, ny); dem.dx, dem.dy in metres
bed = np.nan_to_num(dem.data, nan=float(np.nanmin(dem.data)))   # fill no-data cells
mesh = Mesh2D(nx=bed.shape[0], ny=bed.shape[1], dx=dem.dx, dy=dem.dy)
# ... build and run the solver as above ...
write_geotiff("peak_depth.tif", GeoArray(data=solver.max_depth(), dx=dem.dx, dy=dem.dy,
                                         x0=dem.x0, y0=dem.y0, crs_wkt=dem.crs_wkt))

The grid must be in a projected coordinate system with metre units, and the no-data cells must be filled before the run.

Adding a water level

A tide, a surge or a river stage is a water-surface elevation over time on a set of cells. Mark the cells with a mask and pass the time series:

from geoswe import StageBoundary

coast = np.zeros((nx, ny), dtype=bool)
coast[0, :] = True                           # here: the whole x = 0 edge
low = float(bed[0].min())                    # lowest ground on that edge
tide = StageBoundary.from_mask(coast, mesh, bed,
                               time_s=[0, 3600, 7200], stage_m=[low, low + 1.0, low])
cfg = Config(friction="manning", bc_x="fall", bc_y="fall",
             rainfall_forcing=rain, stage_boundary=tide)
solver = Solver2D(mesh, cfg, np.zeros((3, nx, ny)), bed)
solver.set_manning(manning)
solver.run(t_end=2 * 3600)

After every step the marked cells are set to depth max(0, stage - bed) with zero momentum, and the stage is interpolated linearly in time. The elevations use the datum of the terrain. StageBoundary.from_noaa_csv reads a NOAA CO-OPS gauge file instead of a list.

When the rectangle is mostly not your problem

If the area you care about is a county, a watershed or a coastline inside a much larger rectangle, the compressed active-cell mesh stores and updates only those cells (GPU only). You build the same dense problem, mark the cells to keep, and convert it before running:

from geoswe import CompressedSolver

ii, jj = np.meshgrid(np.arange(nx), np.arange(ny), indexing="ij")
study_area = np.hypot(ii - nx / 2, jj - ny / 2) < 50      # boolean (nx, ny): a disc here

cfg = Config(friction="manning", bc_x="fall", bc_y="fall")   # the rain goes to `cs`, below
solver = Solver2D(mesh, cfg, np.zeros((3, nx, ny)), bed)
solver.set_manning(manning)
solver.set_inside_mask(study_area)           # the cells to keep
cs = CompressedSolver.from_dense(solver)     # packs them; `solver` is used up
cs.set_rain(rain)
cs.run(t_end=2 * 3600)
depth = cs.depth()                           # zero outside the study area

examples/ex08_compressed_mesh.py runs both solvers on the whole patch and compares them.

Things that go wrong

  • Units. Rain in RainfallForcing is in mm/h; Config.rainfall (a constant rate) is in m/s. Everything else is SI.

  • Array order. Arrays are [x, y]. An image-style [row, column] array read with another library needs np.flipud(a).T to become [x, y], which read_geotiff does for you.

  • No friction. A roughness value has no effect until friction="manning" is set. Config warns when you give one without it, and so does set_manning_table; set_manning switches friction on instead.

  • Speed on the GPU. The fastest path needs single precision (the default there), friction with a roughness map of at most 256 distinct values, and rain given as a field or not at all. GEOSWE_VERBOSE=1 prints which path a run takes.