Buffering only the upslope areas of a DEM

OK, here's the solution to the problem proposed in my previous post.  It's not a perfect solution, as it doesn't account for all possible types of terrain (more on that below), but it's pretty close and in the data I have, works pretty well.

The first thing to do is take care of the slope and the 10' elevation change in one swoop. By adjusting the size of our pixels, we can make sure that any change in the slope that meets the 40% rule also meets the elevation change. In other words, what is the minimum Δx (i.e. the resolution) that yields a change of Δy>=10' at an angle θ?

• Because: y=x·sin(θ)
• We substitute: 10= x·sin(tan^-1(.4))
• Solving for x yields x=26.9. We round that to 27 for convenience.

• Now calculate your slope raster with an output resolution of 27ft (or equivalent, depending on your projection).
• Reclass that raster so that everything < 40% is null.
• Using your >=40% layer assign a value to each grouping (in GRASS, this is r.clump; in Arc*, this is the Region Group tool from Spatial Analyst).

This next part is hard because looping pretty much requires a programming language - I used GRASS (wrapped in a shell script), you could use AML or Python if you're an ESRI user.


For each group:

• Find minimum elevation in that group
• Buffer out 300' from the group
• Delete all pixels from the buffered region where slope is null and elevation is less than the minimum for that group.

Ta da!  You have a raster buffered out 300' upslope.

But wait....


What if a pixel is both downslope from one face and upslope from another? (Imagine a ledge or an old lake shoreline like the bench in northern Utah created by Lake Bonneville.)  Well, in this case, we will have eliminated those pixels, which we don't want to do.

We can solve this by taking an intermediary step - creating an empty raster, then instead of deleting the pixels, do a logical OR against the output of the previous iteration of our loop. This allows us to keep our pixels if they're within 300' of any upslope region.

Unfortunately, that loop can add some significant processing time, anyone have another way to do it?