6

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.

4

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.

  • Nice, straightforward answer and, as always, quite informative. – Jeffrey Evans Nov 10 '15 at 17:49
  • Thanks for replying. Apologies in advance if i am wrong, this is quite new to me, but this looks like a way to calculate the maximum slope and then apply an angular correction (using Aspect)? – XMV Nov 11 '15 at 10:04
  • Is this the same than the slope of the original DTM in one direction: as the gradient of the cells adjacent in one direction to the cell we want to calculate slope? – XMV Nov 11 '15 at 10:06
  • @XMV Not quite. It's usually not a good idea to use only the cells adjacent in a given direction: the estimate is too variable and arbitrary (it depends heavily on how you interpolate). The slope and aspect grids are computed in many ways, but the most common (and the one used in ArcGIS) is to fit a least-squares plane to the nine-point neighborhood and find its gradient. However the slope and aspect might be computed, they are equivalent to fitting some plane (linear function) at each point, and the procedure given in my answer returns the slope of that plane in the given direction. – whuber Nov 11 '15 at 14:48
  • Ok. That must explain the differences observed between the method i was using and the one you describe. Basically, i generate a duplicate of the DTM, namely DTM2, then shift the cells a predefined distance (e.g. 30m to the west), then use raster calculator to calculate the linear slope by using the formula: arctan(abs(DTM-DTM2)/ 30). Not very elegant but robust i think! – XMV Nov 12 '15 at 17:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.