Downwinding Workflow Example 2:
Irregular domain boundary and no river mask
Calculating slopes with a downwinding approach to be consistent with ParFlow’s OverlandFlow boundary conditon.
This requires two inputs: (1) a DEM and (2) a watershed mask
This example usese the test domain from Condon and Maxwell (2019) (https://doi.org/10.1016/j.cageo.2019.01.020) the datasets for this domain are provided with the PriorityFlow R package to use your own datasets you should have a DEM and mask files formatted as a matrices with [i,j] corresponding to the x and y axes of the domain (i.e. DEM[0,0] is the lower left corner of the domain and DEM[nx - 1, ny - 1] is the upper right)
import numpy as np
import matplotlib.pyplot as plt
from priority_flow import (
init_queue,
d4_traverse_b,
load_dem,
load_watershed_mask,
drainage_area,
calc_subbasins,
slope_calc_upwind,
)
from parflow.tools.io import read_pfb, write_pfb
# -------------------------------------------------------------------------
# Settings
# -------------------------------------------------------------------------
# DEM processing
ep = 0.01 # epsilon for D4TraverseB (PriorityFlow processing)
# Slope scaling
maxslope = 0.5 # maximum slope; set to -1 to disable
minslope = 1e-5 # minimum slope; set to -1 to disable
scale = 0.1 # max ratio of |secondary| / |primary| (secondaryTH)
# River and subbasin size for slope calculations
sub_th = 100 # area threshold (cells) for subbasin delineation
riv_th = sub_th # optional: threshold for river mask for slope processing
riv_method = 3 # 0: none, 1: scale river secondary, 2: basin mean, 3: reach mean
mrg_th = 10 # merge threshold for small subbasins
# Grid dimensions for slopes
dx = 1000.0
dy = 1000.0
# Run name for outputs
runname = "Test"
# -------------------------------------------------------------------------
# Load DEM and watershed mask
# -------------------------------------------------------------------------
DEM = load_dem()
watershed_mask = load_watershed_mask()
nx, ny = DEM.shape
print(f"Domain size: nx={nx}, ny={ny}")
Domain size: nx=172, ny=215
# Plot the inputs:
fig, axes = plt.subplots(1, 2, figsize=(10, 3))
im0 = axes[0].imshow(DEM, cmap='viridis', origin='lower')
axes[0].set_title("DEM")
plt.colorbar(im0, ax=axes[0])
im1 = axes[1].imshow(watershed_mask, cmap='viridis', origin='lower')
axes[1].set_title("Watershed mask")
plt.colorbar(im1, ax=axes[1])
<matplotlib.colorbar.Colorbar at 0x15332dfe4a40>
Process the DEM:
initialize the queue with all the border cells for the irrgular boundary
Process the DEM so that all cells drain to the borders
init = init_queue(DEM, domainmask=watershed_mask)
No init mask provided all border cells will be added to queue
No border provided, setting border using domain mask
trav_hs = d4_traverse_b(
DEM,
init["queue"].copy(),
init["marked"].copy(),
basins=init["basins"].copy(),
epsilon=ep,
mask=watershed_mask,
)
def _plot_trav_hs(trav_hs):
"""Plot DEM processing results."""
fig, axes = plt.subplots(2, 2, figsize=(10, 8))
im0 = axes[0, 0].imshow(trav_hs["dem"], cmap='viridis', origin='lower')
axes[0, 0].set_title("Processed DEM")
plt.colorbar(im0, ax=axes[0, 0])
im1 = axes[0, 1].imshow(trav_hs["marked"], cmap='viridis', origin='lower')
axes[0, 1].set_title("Marked")
plt.colorbar(im1, ax=axes[0, 1])
im2 = axes[1, 0].imshow(trav_hs["step"], cmap='viridis', origin='lower')
axes[1, 0].set_title("Step")
plt.colorbar(im2, ax=axes[1, 0])
im3 = axes[1, 1].imshow(trav_hs["basins"], cmap='viridis', origin='lower')
axes[1, 1].set_title("Basins")
plt.colorbar(im3, ax=axes[1, 1])
plt.tight_layout()
plt.show()
_plot_trav_hs(trav_hs)
Calculate the slopes
Note this step also fixes the directions of the borders because directions are not provided when the queue is initialized
Option 1: just calculate the slopes for the entire domain with no distinction between river and hillslope cells
In this example secondary slope scaling is turned on and the secondary Slopes in the secondary direction are set to a maximum of 0.1*primary flow direction To calculate only slopes in the primary flow direction set the secondaryTH to 0 Additionally primary slopes are limited by min slope and max slope thresholds
slopes_uw = slope_calc_upwind(
dem=trav_hs["dem"].copy(),
direction=trav_hs["direction"].copy(),
dx=dx,
dy=dy,
mask=watershed_mask.copy(),
minslope=minslope,
maxslope=maxslope,
secondary_th=scale,
)
upwinding slopes
Limiting slopes to maximum absolute value of 0.5
Limiting the ratio of secondary to primary slopes 0.1
WARNING: 3 Flat cells found
After processing: 0 Flat cells left
Limiting slopes to minimum 1e-05
Option 2: If you would like to handle river cells differently from the rest of the domain
# Calculate the drainage area
area = drainage_area(trav_hs["direction"], printflag=False)
Define subbasins for calcualting river reach slopes the riv_th here is the drainage area threshold for splitting the river network branches when you do this you can still end up with subbasins with drainage areas less than the riv_th when multiple branches come together in a small area. To fix this you can set a merge threshold (merge_th) so that subbains with areas < merge_th autmoatically get merged with their downstream neighbor.
subbasin = calc_subbasins(
trav_hs["direction"].copy(),
area=area.copy(),
mask=watershed_mask.copy(),
riv_th=sub_th,
merge_th=mrg_th,
)
WARNING: non-zero merge thresholds are not compatible with the RiverSmooth function
temp = subbasin['RiverMask'].copy()
temp[temp == 0] = np.nan
fig, ax = plt.subplots(figsize=(6, 4))
# Base layer
im0 = ax.imshow(subbasin['subbasins'],
cmap='viridis',
origin='lower')
# Overlay (rivers)
im1 = ax.imshow(temp * 2,
cmap='Reds', # different colormap helps visibility
origin='lower',
alpha=0.8) # transparency like add=TRUE overlay
ax.set_title("Subbasins + Rivers")
plt.colorbar(im0, ax=ax, fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()
# Calculate the slopes
# The "river_method' flag here determines how the river cells will be handeled
# (e.g. using subbasin averages along reaches). Refer to the top of this script or the function for details.
slopes_uw = slope_calc_upwind(
dem=trav_hs["dem"].copy(),
direction=trav_hs["direction"].copy(),
dx=dx,
dy=dy,
mask=watershed_mask.copy(),
minslope=minslope,
maxslope=maxslope,
secondary_th=scale,
river_method=riv_method,
rivermask=subbasin["RiverMask"].copy(),
subbasins=subbasin["subbasins"].copy(),
)
upwinding slopes
Limiting slopes to maximum absolute value of 0.5
Limiting the ratio of secondary to primary slopes 0.1
River Method 3: assigning average river slope to river cells by watershed
Scaling secondary slopes along river mask to 0 * primary slope
WARNING: 3 Flat cells found
After processing: 0 Flat cells left
Limiting slopes to minimum 1e-05
Option 2b: Alternate more advanced approach:
Define a river mask separate from the subbasin river mask and use this for the slope calculations. If you do this the average slopes will still be calculated by subbasin using the sub_th, but you can apply those average sloeps to more river cells by setting a lower threshold here. This is the ‘riv_th’ set at the top if you set riv_th=sub_th at the top this will have the same effect as just running the slope calc with the subbasin['RiverMask'].
rivers = np.zeros_like(area)
rivers[area < riv_th] = 0
rivers[area >= riv_th] = 1
# Plot the subbasins with the new river mask to check that the threshold is good
temp = rivers.copy()
temp[temp == 0] = np.nan
fig, ax = plt.subplots(figsize=(6, 4))
# Base layer
im0 = ax.imshow(subbasin['subbasins'],
cmap='viridis',
origin='lower')
# Overlay (rivers)
im1 = ax.imshow(temp * 2,
cmap='Reds', # different colormap helps visibility
origin='lower',
alpha=0.8) # transparency like add=TRUE overlay
ax.set_title("Subbasins + Rivers")
plt.colorbar(im0, ax=ax, fraction=0.046, pad=0.04)
plt.tight_layout()
plt.show()
slopes_uw = slope_calc_upwind(
dem=trav_hs["dem"].copy(),
direction=trav_hs["direction"].copy(),
dx=dx,
dy=dy,
mask=watershed_mask.copy(),
minslope=minslope,
maxslope=maxslope,
secondary_th=scale,
river_method=riv_method,
rivermask=rivers,
subbasins=subbasin["subbasins"].copy(),
)
upwinding slopes
Limiting slopes to maximum absolute value of 0.5
Limiting the ratio of secondary to primary slopes 0.1
River Method 3: assigning average river slope to river cells by watershed
Scaling secondary slopes along river mask to 0 * primary slope
WARNING: 3 Flat cells found
After processing: 0 Flat cells left
Limiting slopes to minimum 1e-05
#Look at the slopes and directions
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
im0 = axes[0].imshow(slopes_uw['slopex'], cmap='viridis', origin='lower')
axes[0].set_title("slopex")
plt.colorbar(im0, ax=axes[0])
im1 = axes[1].imshow(slopes_uw['slopey'], cmap='viridis', origin='lower')
axes[1].set_title("slopey")
plt.colorbar(im1, ax=axes[1])
im2 = axes[2].imshow(slopes_uw['direction'], cmap='viridis', origin='lower')
axes[2].set_title("direction")
plt.colorbar(im2, ax=axes[2])
plt.show()
# Calculate the drainage area - if you went with option 1 for slopes and you didn't do this already
area = drainage_area(trav_hs["direction"], mask=watershed_mask, printflag=False)
# Save slope data as PFB files:
write_pfb("slopex.pfb", slopes_uw["slopex"])
write_pfb("slopey.pfb", slopes_uw["slopey"])