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!
No products are associated with this question.
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.