The problem with doing complete filling on flats is that the resulting DEM contains mathematically flat areas with no local gradient. This makes it impossible to determine flow directions for these areas.
Imposing an epsilon gradient during depression-filling is one solution as discussed in Epsilon Filling (Epsilon Filling); however, the gradient produced may appear unaesthetic because drainage takes a least-distance route to the flat’s edges.
We stress the word aesthetic here since, after depression-filling, no information about local gradients remains from the original DEM, so, in a sense, all reconstructed drainage patterns are equally silly. Sometimes, though, this is the best you can do.
This pages discusses alternatives.
Barnes (2014) Flat Resolution¶
Barnes, R., Lehman, C., Mulla, D., 2014a. An efficient assignment of drainage direction over flat surfaces in raster digital elevation models. Computers & Geosciences 62, 128–135. doi:10.1016/j.cageo.2013.01.009
The Barnes (2014) flat resolution algorithm routes flow both towards the edges of flats and away from the high areas surrounding them. The result is an aesthetically pleasing drainge pattern.
It can do so either by adjust the elevations of the DEM’s cells or by adjusting flow metrics derived from the DEM.
In this method, the elevation of a DEM can be adjusted so that every cell in the region is raised some small amount, ε, above cells which are closer to a depression’s spill point and farther from its surrounding high areas.
This must be done carefully. In floating-point DEMs, the value ε is non-constant
and must be chosen using the !`std::nextafter` function. If a depression is too
large, the imposed gradient may result in the interior of the depression being
raised above the surrounding landscape. Using
double instead of
reduces the potential for problems at a cost of twice the space used. If a
problem does arise, RichDEM provides a warning.
Recall from Epsilon Filling (Epsilon Filling) that an epsilon gradient imposed during depression-filling results in an elevation adjustment that looks like this:
import richdem as rd import numpy as np #Load dataset beau = rd.rdarray(np.load('imgs/beauford.npz')['beauford'], no_data=-9999) #Fill the depression entirely beau_filled = rd.FillDepressions(beau, epsilon=False, in_place=False) #Construct the epsilon drainage surface via filling beau_eps = rd.FillDepressions(beau, epsilon=True, in_place=False) diff = beau_eps - beau_filled rd.rdShow(diff, ignore_colours=, axes=False, cmap='jet', figsize=(8,5.5))
In contrast, the Barnes (2014) convergent elevation adjustment looks like this:
#Resolve flats by imposing a convergent epsilon gradient beau_flat_eps = rd.ResolveFlats(beau_filled, in_place=False) diff = beau_flat_eps - beau_filled rd.rdShow(diff, ignore_colours=, axes=False, cmap='jet', figsize=(8,5.5))
The difference versus the depression-filling epsilon adjustment appears as follows. Note the deep V-shaped notches in the flats indicating the increased convergence of the Barnes (2014) method.
diff = beau_flat_eps - beau_eps rd.rdShow(diff, ignore_colours=, axes=False, cmap='jet', figsize=(8,5.5))
Flow Metric Adjustment¶