Which tolerance should I choose when using knickpointfinder?

Posted on

Kai Deng from the GFZ Potsdam commented on my previous post on knickpointfinder and asked the question about the parameter ‘tol’:

[…] one question on the parameter “tol”. You said that “the value of tol should reflect uncertainties that are inherent in longitudinal river profile data”. However, I’m still confused about how to translate this sentence into code, or in another word, how can I validate my choice of tol by some criteria based on codes?

Now that’s an excellent question. The ‘tol’-parameter determines how many knickpoints the finder will detect and we might not want to miss a knickpoint. At the same time, we know that river profile data often has some erroneous values and we wouldn’t want to pick knickpoints that are merely artefacts of the data. So, what value of ‘tol’ should we choose?

I wrote that “the value of tol should reflect uncertainties that are inherent in longitudinal river profile data”. More specifically, I would argue that one should choose tolerance values that are higher than the maximum expected error between the measured and the true river profile. With lower tolerance values we would likely increase the false positive rate, i.e., the probability that we choose a knickpoint that is due to an artefact.

Now, the problem that we immediately face is that we don’t know this error. Of course, there are quality measures for most common DEM products but usually these do not apply to values measured along valley bottoms, at least if you take common global DEM products such as SRTM, ASTER GDEM, etc (Schwanghart and Scherler, 2017). Usually, there are also no repeated measurements. Hence, we cannot estimate the error as we would in the case of measurements in the lab.

But there is another approach: What we know for sure is that river profiles should monotonously decrease downstream. At least, this is true for the water surface, yet not necessarily bed elevations (e.g. due to the presence of riffle-pool sequences). The elevation of a pixel on the stream network should thus be either lower or equal to the elevation of its upstream neighbor. If it is not, than the value is probably wrong (for whatever reason). I say ‘probably’, because the upstream pixel might have an erroneously low value, too.

The following code demonstrates how we could determine these errors: Through downstream minima imposition (imposemin) and upstream maxima imposition (cummaxupstream). The shaded areas in the thus generated plot (here I show a zoomed version) show the differences between these values.

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM);
S  = STREAMobj(FD,'minarea',1000);
S = klargestconncomps(S);
zmax = cummaxupstream(S,DEM);
zmin = imposemin(S,getnal(S,DEM));
plotdzshaded(S,[zmin zmax])
hold on
Shaded areas show the difference upstream maxima and downstream minima imposition.

The maximum difference between the elevations is the maximum relative error that we estimate from the profile.


ans =


Now we might want to take this elevation difference as tolerance value for the knickpointfinder. However, this value might be very high if our profile is characterized by a few positive or negative outliers. This would result in only few if any knickpoints. Thus, we could come up with a more robust estimate of the error which is based on quantiles. Robust in this sense means: robust to outliers. And the proposed technique is  constrained regularized smoothing (CRS) (Schwanghart and Scherler, 2017).

What is constrained regularized smoothing? It is an approach to hydrologically correct DEMs based on quantile regression techniques. While minima or maxima imposition present two extremes of hydrological conditioning – carving and filling, respectively – CRS (and quantile carving, which is similar but doesn’t smooth) is a hybrid technique and returns river longitudinal profiles that run along the tau’s quantile of the data. Why would this be preferred? Simply because we would not want that the shape of the resulting profile is predominantly determined by the errors with the largest magnitude.

Here is how to apply quantile carving using the quantiles 0.1 and 0.9.

z90 = crs(S,DEM,'K',1,'tau',0.9);
z10 = crs(S,DEM,'K',1,'tau',0.1);

ans =


Of course, the maximum difference is lower than the previous value obtained by minima and maxima imposition. But not very much. Let’s see how well the algorithm performs if we add some additional errors.

IX = randlocs(S,10);
IX = find(ismember(S.IXgrid,IX));
z = getnal(S,DEM);
z(IX) = z(IX) + sign(randn(numel(IX),1))*100;


River profile with some additional errors.

Let’s compute the difference between carving and filling.

zmin = imposemin(S,z);
zmax = cummaxupstream(S,z);

ans =


And whatabout our uncertainty estimate based on the CRS algorithm?

z90 = crs(S,z,'K',1,'tau',0.9);
z10 = crs(S,z,'K',1,'tau',0.1);
d = max(z90-z10)

d =


The CRS result is not affected at all by few outliers. That means, it enables us to derive a rather robust estimate about the uncertainty of the data.

Moreover, we can use the CRS-corrected profile as input to the knickpointfinder, because we have filtered the spikes from the data. As elevations are often overestimated, you might use directly the node-attribute list z10.

I have shown that you can use an uncertainty estimate derived from your river profile as direct input to knickpointfinder. There are, however, a few things to note:

  1. Don’t carve or fill your data before using CRS. Otherwise, your uncertainty estimate is flawed.
  2. The uncertainty estimate may not reliably detect the true uncertainty in your data. Alsways cross-check whether the uncertainty estimate is either too low or too high.
  3. I haven’t talked about the smoothness parameter ‘K’ which also controls the uncertainty estimates. Larger values of K will likely increase the uncertainty estimate.
  4. CRS can take a while for very large river networks. Consider estimating the tolerance value from a subset of the data.
  5. You can also use a node-attribute list of tolerances in knickpointfinder and thus account for spatially varying tolerances. In this case, however, add some offset because the CRS might indicate areas with zero uncertainty and thus overestimate the true precision of the data.


Schwanghart, W., Scherler, D., 2017. Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques. Earth Surface Dynamics, 5, 821-839. [DOI: 10.5194/esurf-5-821-2017]


Smoothing the planform geometry of river networks

Posted on

In previous posts (here, here, here), I have covered some of the tools and applications of smoothing in the analysis of river networks. Thereby, I have focused on the analysis of river profiles (see also our paper here). However, these tools can also be used to smooth and analyse the planform patterns of river networks.

Here is how:

DEM = GRIDobj('srtm_bigtujunga30m_utm11.tif');
FD = FLOWobj(DEM,'preprocess','carve');
S  = STREAMobj(FD,'minarea',1000);
S  = klargestconncomps(S);
y = smooth(S,S.y,'k',500,'nstribs',true);
x = smooth(S,S.x,'k',500,'nstribs',true);

[~,~,xs,ys] = STREAMobj2XY(S,x,y);
Upper panel: Original river network data. Lower panel: Smoothed river network.

Now what is this good for? First, simplification by smoothing might help in visualization. Second, the smoothed river coordinates can be used to derive metrics such as planform curvature of river networks. These metrics are strongly affected by the circumstance that river networks derived from digital elevation models are forced to follow the diagonal and cardinal directions of the grid. This results in local flow directions in steps of multiples of 45°. Changes in flow direction thus occur in jumps, and accordingly planform curvature of river networks are highly erratic. These erratic values may hide some of the larger-scale curvature patterns that we might be interested in.

The STREAMobj/curvature function features an example that shows how to accentuate larger-scale curvature patterns using smoothing. Following up on the code above, deriving curvature from smoothed planform coordinates can be done as follows:

y = smooth(S,S.y,'k',100,'nstribs',true);
x = smooth(S,S.x,'k',100,'nstribs',true);
c = curvature(S,x,y);
box on
% center 0-value in colormap
caxis([-1 1]*max(abs(caxis)))
h = colorbar;
h.Label.String = 'Curvature';
Planform curvature of the river networks.

Clearly, the amount of smoothing determines the spatial scale of curvature. Here, the chosen smoothing parameter (K=100) accentuates the bends at a scale of a few hundred meters.

There are a number of applications for this kind of analysis. Lateral erosion of rivers may undercut slopes and it may be interesting to investigate if planform curvature metrics help to detect and predict locations of landsliding. Just an idea…

Paper out now: Bumps in river profiles

Posted on

Different preprocessing techniques and their effect on the longitudinal river profile (Figure 1 from Schwanghart and Scherler (2017)).

Our paper on uncertainty quantification and smoothing of longitudinal river profiles has now been published in Earth Surface Dynamics. The paper describes a new approach to hydrologically correct and smooth profiles using quantile regression. Smoothing is based on curvature regularization with the constraint that elevation must decrease with increasing flow distance. We thus refer to this technique as constrained regularized smoothing (CRS). We compare CRS-derived profiles to profiles obtained from the common methods of filling and carving, and show that CRS outperforms these methods in terms of accuracy and precision.

Check out the new TopoToolbox functions that accompany the paper:

  • STREAMobj/crs
  • STREAMobj/crslin
  • STREAMobj/crsapp
  • STREAMobj/smooth
  • STREAMobj/quantcarve
  • FLOWobj/quantcarve

Applications of these functions were already described in the posts on An alternative to along-river swath profiles, Steepness derived from smoothed river profiles, Smooth operator….


Schwanghart, W., Scherler, D., 2017. Bumps in river profiles: uncertainty assessment and smoothing using quantile regression techniques. Earth Surface Dynamics, 5, 821-839. [DOI: 10.5194/esurf-5-821-2017]