MATLAB Answers


How can I (quickly) buffer a river centerline with constant width to compute banks?

Asked by Jon
on 5 Feb 2013
Latest activity Answered by Chad Greene
on 26 Feb 2015

I am developing a meandering river model. At each time step, centerline nodes are adjusted by a dx and dy value and the river moves about its floodplain. At some point, the bend of a river runs into another bend of the same river, causing a cutoff. In reality, the cutoff will occur when the banks intersect, not the centerline.

My problem is that my bank calculation method fails for segments with very high curvature (e.g. those that have just undergone a cutoff), but works sufficiently well elsewhere. I cannot figure out how else to generate my right and left banks to avoid this problem. for some images of the problem.

Here's what I'm doing:

Given a vector of X,Y centerline coordinates, I first calculate the downstream angle between the river centerline and the (arbitrary) valley centerline:

atan((Y2-Y1))/(X2-X1)). Actual code is

Xang = Xs(1:end-1);
Xang_plus1 = Xs(2:end);
Yang = Ys(1:end-1);
Yang_plus1 = Ys(2:end);
angles = atan2(Yang_plus1-Yang,Xang_plus1-Xang);

Once I compute these angles, I calculate the left and right banks using the following code:

Xl = Xs(1:end-1)-sin(pi-angles)*B;
Yl = Ys(1:end-1)-cos(pi-angles)*B; 
Xr = Xs(1:end-1)+sin(pi-angles)*B;
Yr = Ys(1:end-1)+cos(pi-angles)*B; 

Like I said before, this code works very well except where the curvature is high, ie where the river changes directions drastically.

I've tried different thresholding, smoothing, and regridding techniques but can't find a robust method. I've looked at bufferm, but not only did it not work correctly (user error, I'm sure), it took ~10 seconds. This is a routine I need to run at least twice per time step over 1000s of time steps, so a 10s subroutine is unacceptable.

Does anyone have any suggestions for how I can generate my riverbanks accurately and quickly, or do I have to resort to curvature-based thresholding?

Thanks for any help!


Image Analyst
on 5 Feb 2013

Some screenshots would help illustrate your situation.

on 5 Feb 2013

Thanks for the advice; I uploaded a few screenshots here:


No products are associated with this question.

2 Answers

Answer by Cedric Wannaz
on 5 Feb 2013
Edited by Cedric Wannaz
on 5 Feb 2013
 Accepted answer

Your approach generates narrower streams than what you want, and will fail each time the angle is greater than pi/2, as at least one point delineating your buffer will essentially switch side(s).

bufferm() was known to have flaws but most of them were corrected. You could have a look at bufferm2 though. I've never used them to be honest (I am using the ArcGIS geoprocessor and more recently arcpy to do these things), but if I had to place a wild guess at what might not work, I would bet on the fact that your segments are not polygons, but lines. Some packages for building buffers are able to deal with points/lines, but some aren't (they need, among other things, to be able to define polygons interior(s) and an exterior(s)).

Otherwise, what you want is (I guess) the intersect between parallels to each segment of your center line. You could get them easily by solving linear systems (intersecting lines) with only one particular case to manage: co-linearity between consecutive segments of the center line (note that there are two sub-cases [aligned and superimposed] that cannot be managed the same way). You might also want to truncate boundaries (the exterior ones) for angles greater than pi/2.


Cedric Wannaz
on 6 Feb 2013

The issue that even parallels won't solve is that your B is much greater than the length of segments from the center line. I have no easy solution for that (but I am still thinking about it). I tested bufferm; it works well but it is true that it is quite slow.

on 6 Feb 2013

Yes, I got the bufferm2 to work and it did work very nicely, but it takes much too long. I'm working on discriminating between erroneous bank intersections and actual cutoff bank intersections now.

on 25 Feb 2015

If anyone else has a similar problem, I ultimately ended up using rangesearch (Stats toolbox) to find cutoffs. Check out the supplementary code with this paper: 10.1002/2014JF003252

Answer by Chad Greene
on 26 Feb 2015

The xy2sn function may be useful. Convert your river xy coordinates to along-flow and cross-flow components, make two lines as the original line +/- some cross-flow offset, then convert those two lines back to xy coordinates.


Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply today