Note
This page was generated from a jupyter notebook.
Introduction to priority flood component¶
The priority flood flow director is designed to calculate flow properties over large scale grids. In the following notebook we illustrate how flow accumulation can be calculated for a real DEM downloaded with the BMI_topography data component. Moreover, we demonstrate how shaded relief can be plotted using the imshowhs_grid function.
First we will import all the modules we need.
[ ]:
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
from bmi_topography import Topography
from landlab import imshowhs_grid
from landlab.components import FlowAccumulator, PriorityFloodFlowRouter
from landlab.io import esri_ascii
Create a function to download and save SRTM images using BMI_topography.
[ ]:
def get_topo(buffer, north=40.16, south=40.14, east=-105.49, west=-105.51):
params = Topography.DEFAULT.copy()
params["south"] = south - buffer
params["north"] = north + buffer
params["west"] = -105.51 - buffer
params["east"] = -105.49 + buffer
params["output_format"] = "AAIGrid"
params["cache_dir"] = Path.cwd()
dem = Topography(**params)
name = dem.fetch()
# props = dem.load()
# dim_x = props.sizes["x"]
# dim_y = props.sizes["y"]
# cells = props.sizes["x"] * props.sizes["y"]
with open(name) as fp:
grid = esri_ascii.load(fp, name="topographic__elevation", at="node")
return grid
Make function to plot DEMs and drainage accumulation with shaded relief.
[ ]:
def plotting(
grid, topo=True, DA=True, hill_DA=False, flow_metric="D8", hill_flow_metric="Quinn"
):
if topo:
azdeg = 200
altdeg = 20
ve = 1
plt.figure()
plot_type = "DEM"
ax = imshowhs_grid(
grid,
"topographic__elevation",
grid_units=("deg", "deg"),
var_name="Topo, m",
cmap="terrain",
plot_type=plot_type,
vertical_exa=ve,
azdeg=azdeg,
altdeg=altdeg,
default_fontsize=12,
cbar_tick_size=10,
cbar_width="100%",
cbar_or="vertical",
bbox_to_anchor=[1.03, 0.3, 0.075, 14],
colorbar_label_y=-15,
colorbar_label_x=0.5,
ticks_km=False,
)
if DA:
# %% Plot first instance of drainage_area
grid.at_node["drainage_area"][grid.at_node["drainage_area"] == 0] = (
grid.dx * grid.dx
)
plot_DA = np.log10(grid.at_node["drainage_area"] * 111e3 * 111e3)
plt.figure()
plot_type = "Drape1"
drape1 = plot_DA
thres_drape1 = None
alpha = 0.5
cmap1 = "terrain"
ax = imshowhs_grid(
grid,
"topographic__elevation",
grid_units=("deg", "deg"),
cmap=cmap1,
plot_type=plot_type,
drape1=drape1,
vertical_exa=ve,
azdeg=azdeg,
altdeg=altdeg,
thres_drape1=thres_drape1,
alpha=alpha,
default_fontsize=12,
cbar_tick_size=10,
var_name="$log^{10}DA, m^2$",
cbar_width="100%",
cbar_or="vertical",
bbox_to_anchor=[1.03, 0.3, 0.075, 14],
colorbar_label_y=-15,
colorbar_label_x=0.5,
ticks_km=False,
)
props = dict(boxstyle="round", facecolor="white", alpha=0.6)
textstr = flow_metric
ax.text(
0.05,
0.95,
textstr,
transform=ax.transAxes,
fontsize=10,
verticalalignment="top",
bbox=props,
)
if hill_DA:
# Plot second instance of drainage_area (hill_drainage_area)
grid.at_node["hill_drainage_area"][grid.at_node["hill_drainage_area"] == 0] = (
grid.dx * grid.dx
)
np.log10(grid.at_node["hill_drainage_area"] * 111e3 * 111e3)
plt.figure()
plot_type = "Drape1"
# plot_type='Drape2'
drape1 = np.log10(grid.at_node["hill_drainage_area"])
thres_drape1 = None
alpha = 0.5
cmap1 = "terrain"
ax = imshowhs_grid(
grid,
"topographic__elevation",
grid_units=("deg", "deg"),
cmap=cmap1,
plot_type=plot_type,
drape1=drape1,
vertical_exa=ve,
azdeg=azdeg,
altdeg=altdeg,
thres_drape1=thres_drape1,
alpha=alpha,
default_fontsize=10,
cbar_tick_size=10,
var_name="$log^{10}DA, m^2$",
cbar_width="100%",
cbar_or="vertical",
bbox_to_anchor=[1.03, 0.3, 0.075, 14],
colorbar_label_y=-15,
colorbar_label_x=0.5,
ticks_km=False,
)
props = dict(boxstyle="round", facecolor="white", alpha=0.6)
textstr = hill_flow_metric
ax.text(
0.05,
0.95,
textstr,
transform=ax.transAxes,
fontsize=10,
verticalalignment="top",
bbox=props,
)
Compare default Landlab flow accumulator with priority flood flow accumulator¶
For small DEMs (small buffer size, in degrees), the default flow accumulator is slightly faster than the priority flood flow accumulator. For large DEMs, the priority flood flow accumulator outperforms the default flow accumulator by several orders of magnitude. To test the performance for larger DEM’s increase the buffer size (e.g. with 1 degree = 111 km).
Default flow director/accumulator¶
[ ]:
# Download or reload topo data with given buffer
# dim_x, dim_y, cells, grid_LL, z_LL, dem = get_topo(0.05)
grid_LL = get_topo(0.05)
fa_LL = FlowAccumulator(
grid_LL, flow_director="D8", depression_finder="DepressionFinderAndRouter"
)
fa_LL.run_one_step()
# Plot output products
plotting(grid_LL)
[ ]:
north = 40.16
south = 40.14
east = -105.49
west = -105.51
buffer = 0.05
params = Topography.DEFAULT.copy()
params["south"] = south - buffer
params["north"] = north + buffer
params["west"] = -105.51 - buffer
params["east"] = -105.49 + buffer
params["output_format"] = "AAIGrid"
params["cache_dir"] = Path.cwd()
dem = Topography(**params)
name = dem.fetch()
# props = dem.load()
name
Priority flood flow director/accumulator¶
Calculate flow directions/flow accumulation using the first instance of the flow accumulator
[ ]:
# Download or reload topo data with given buffer
# dim_x, dim_y, cells, grid_PF, z_PF, dem = get_topo(0.05)
grid_PF = get_topo(0.05)
# Here, we only calculate flow directions using the first instance of the flow accumulator
flow_metric = "D8"
fa_PF = PriorityFloodFlowRouter(
grid_PF,
surface="topographic__elevation",
flow_metric=flow_metric,
suppress_out=False,
depression_handler="fill",
accumulate_flow=True,
)
fa_PF.run_one_step()
# Plot output products
plotting(grid_PF)
Priority flood flow director/accumulator¶
Calculate flow directions/flow accumulation using the second instance of the flow accumulator
[ ]:
# 3. Priority flow director/accumualtor
# Download or reload topo data with given buffer
# dim_x, dim_y, cells, grid_PF, z_PF, dem = get_topo(0.05)
grid_PF = get_topo(0.05)
# For timing compare only single flow
flow_metric = "D8"
hill_flow_metric = "Quinn"
fa_PF = PriorityFloodFlowRouter(
grid_PF,
surface="topographic__elevation",
flow_metric=flow_metric,
suppress_out=False,
depression_handler="fill",
accumulate_flow=True,
separate_hill_flow=True,
accumulate_flow_hill=True,
update_hill_flow_instantaneous=False,
hill_flow_metric=hill_flow_metric,
)
fa_PF.run_one_step()
fa_PF.update_hill_fdfa()
# 4. Plot output products
plotting(grid_PF, hill_DA=True, flow_metric="D8", hill_flow_metric="Quinn")
# Remove downloaded DEM. Uncomment to remove DEM.
# os.remove(dem.fetch())