[GIS] Directional slope calculation

arcgis-10.0arcgis-desktoprasterslope

I'm trying to generate rasters of slope and curvature of an elevation DTM. The particularity is that i need to be able to specify the direction (e.g 180, 270 degrees)in this calculation.

I am using ArcGIS Desktop 10, however if the solution comes from another software package or code it's fine.

Best Answer

The direction is best described as a vector v = (v1, v2) = (cos(t), sin(t)) where t is the direction angle. Precompute these two components.

Next, after computing the slope and aspect grids for the original DEM, convert the aspect grid into a grid for its cosine (call it [C]) and another grid for its sine (call it [S]). Compute the linear combination

-v1 * [C] - v2 * [S]

This grid gives the cosines of the angles between v and the aspects. Multiplying them by the slope grid (with values expressed as a tangent or a "rise:run" ratio, not in degrees) gives the desired answer. (If you need degrees, take the inverse tangent of the result.)

(To see why this formula works, notice that wherever v points in the same direction as the aspect, this linear combination equals -cos(t)^2 + -sin(t)^2 = -1, thereby negating the original slope and correctly indicating that v points straight downhill. When v points directly opposite the aspect, the linear combination is cos(t)^2 + sin(t)^2 = 1, resulting in positive slopes to designate uphill directions. When v is perpendicular to the aspect, the linear combination is zero, giving a zero slope to reflect the fact this direction is along a contour line. An easy mathematical proof of the formula begins with the definition of a directional derivative.)

Finally, slopes of zero usually correspond to NoData or special codes in the aspect grid (such as -1). Use a conditional calculation to replace all cells in your answer with zeros at any location where the aspect has such a special code.

Make sure you and the software agree about how angles are measured. Many GISes are schizophrenic: they may report aspects in degrees but require radians as arguments to trig functions. ArcGIS is one of them. You may need to make some conversions. If not, you will get bizarre results.

Related Question