While looking at a roadside cutting in order to evaluate options for treating it for erosion, I thought it would be interesting to
gather the data and develop a 3D model in R. The process can be confusing the first time, as the data need to be massaged to meet the requirements of the interpolative algorithm that converts scattered heights to a regular grid, which can then be visualized.

Below is the process I used, and it was relatively quick.

Step 1. First I walked a transect every 25m along the road collecting distances and slopes at every breakpoint in the slope. I can collect distances and slopes quickly this way,
and then use the slopes (s1, s2, …) and distances (d1, d2, …) to calculate heights at each point. E.g. the data looks like this.

x  s1  d1 s2  d2  s3  d3
0 -0.1 2.2  0.33  5.2  0.25 11.5
25 -0.1  2.2  0.32  5.4  0.25  12.3
...

Step 2. Then (using Excel though I could have used R) I converted the slopes and distances to an x, y, and height at each point, e.g.

Dist	Batter	Height
0	2.2	-0.0308
25	2.2	-0.0308
50	2.7	-0.0378
75	3	-0.042
100	2.9	-0.0406
125	3.2	-0.0448
150	3.5	-0.049
175	5.5	-0.077
200	5.2	-0.0728
225	5.2	-0.0728
250	5.2	-0.0728
0	2.2	-0.0308
25	7.3	1.2952
50	8.9	2.0888
75	10.4	2.4592
100	10.8	1.9344
...

The data above can be read into R and plotted up using the short R script below. Read.table reads in the tabular data, interp creates a regular grid from the triples (the xi$z variable) and persp creates the image of the cutting below. The parameters to persp need to be adjusted to give a point of view that works for you.

xs< -read.table("Embankment-Xsection.txt",header=T)
xi<-interp(xs$Dist,xs$Batter,xs$Height)
persp(xi$x,xi$y,xi$z,scale=F,expand=2,phi=5,theta=80)

Next, how to calculate and map flow accumulation lines.