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

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.

http://imgur.com/a/uHvRM 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]

Some screenshots would help illustrate your situation.

[Jon]

Thanks for the advice; I uploaded a few screenshots here: http://imgur.com/a/uHvRM

Answer 1

[Cedric]

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.

[Jon]

Unfortunately the bufferm command is just too time consuming. According to the help file, it accepts lines defined by lat,lon (or X,Y, I assume) which is what I'm inputting. The output is a field of about 10 randomly sized and spaced polygons that don't cover the domain of the inputs...??

Anyway, I'm interested in your intersecting parallels suggestion. That seems exactly like what I'd want, but I'm not sure on how to go about coding it. I already am implementing a code that finds intersections (to identify the cutoffs) so perhaps I can recycle that bit. And I guess generating parallels should be no problem...I'll give it a shot! Thanks for the idea.

[Cedric]

Coding polygons vertices coordinates is a little subtle.. there is a convention in the way you number external and internal boundaries (clockwise/couterclockwise), that allows the presence of holes for example. I don't know bufferm, but I guess that it will work that way. You should send coordinates of a "square" in lon/lat, and see what you get in return. If it needs a polygon and is not too sensitive to null area polygons, you could just send, e.g P1, P2, P3, P4, P3, P2, P1 and see what you get.

If you really need to implement something by yourself, you already compute slopes of all center line segments. Parallels will have same slope +/- but some offset that you can compute (managing "vertical" cases separately). Then, for two consecutive segments, you find the intersect just by solving a system of 2 eqs with two unknowns. This can certainly be vectorized.

[Jon]

Nice. Thanks for the help!

[Cedric]

Looking at your figures, I had another idea if you don't want to re-implement code.. it seems that you already implemented code that detects intersecting/overlapping banks. You could reuse that to cut loops in banks after your cutoff operations, typically, when segment [p to p+1] of the bank intersects segment [p+n to p+n+1], you eliminates points p+1 to p+n and replace them by a new point that is the intersect..

[Jon]

Yes, I thought about that. I'm not modeling bank migration so I have to recompute the banks at each time step. The problem is that for the next time step after a cutoff, the bank is recomputed and the error is still there even though there wasn't a cutoff that time step. So the code would think the error is an actual cutoff and remove nodes from the centerline. Hope that makes sense.

edit for clarity:

1. migrate centerline

2. compute banks

3. if banks intersect themselves, cutoff the centerline

4. recompute banks, check 3 again (multiple cutoffs in the same time step)

5. if no cutoffs, go back to 1.

The problem is in step 3--the erroneous banks will be interpreted as cutoffs and the centerline will be modified when it shouldn't.

[Cedric]

What about the following (pseudo-code)..

 for t in timesteps
    move center line
    recompute banks
    detect/flag intersects/contacts between banks segments
    cut chunks of center line between flagged segments and merge
    % Until this point, I guess that this is what you are doing.
    % Add this part..
    detect/flag intersects/contacts between banks segments (2nd time)
    replace all banks points between flagged segments/points with the  
      coordinates of the intersect.
 end

This way you would keep a number of points in banks that would match the number of points in the center line, and all points that define the loop of the bank in your 3rd figure would be replaced by the intersect (superimposed).

[Jon]

The problem is that it's possible and not entirely improbable to have two cutoffs at different sections of the river occur on the same timestep. I can work around this, though, but the problem I mentioned still remains: at the time step following a cutoff, when the banks are recomputed near the cutoff, they still erroneously intersect themselves. So the code thinks there's a cutoff when there's actually not a cutoff, and I have no robust way of informing the code that it's not a real cutoff...although typing this has given me an idea.

I would rather just compute the banks properly, because not only would it solve my problem of erroneous cutoffs, it would also make the results look more realistic for presentation purposes.

Also, I just tried the intersecting parallel lines method and the problem still remains at the regions of high curvature (cutoff nodes). You can hand draw a sharp bend on some paper and see why.

edit: added a screenshot of the failed method here: http://imgur.com/a/3VvYg

[Cedric]

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.

[Jon]

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.

[Jon]

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 2

[Chad Greene]

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.