I wanted to create a DEM of constant slope say 1/3000.
How can I proceed in ArcGIS Desktop 9.3?
The area may be 1000m x 1000m or any other as per convenience.
Add a comment
|
1 Answer
A solution to this is discussed on the ESRI forum here:
Every planar surface has a simple formula, Jesse, in the form
(something)*x + (something else)*y + (constant).
In the Raster Calculator, x is $$XMap and y is $$YMap. The trick is figuring out the coefficients. Often you want to specify a plane in terms of its slope (s), aspect (a), and height at some given point ("base elevation"). Its equation will be in the form
-tan(a)sin(s)[x] - tan(a)cos(s)[y] + (constant).
To work out that constant, plug in the coordinates where you know the base elevation, say (x0, y0), and then adjust the constant to make it work out. This produces the formula
Base elevation + tan(a)sin(s)[x0] + tan(a)cos(s)[y0] = constant.
For example, to produce a plane sloping s=30 degrees with an northeastern aspect a=45 degrees, compute
tan(a) = tan(45 degrees) = 1 sin(s) = sin(30 degrees) = 1/2 cos(s) = cos(30 degrees) = Sqrt(3)/2 = about 0.866.
Therefore the plane's height will be in the form
$$XMap * (-0.5) - $$YMap * (0.866) + (constant)
and you can work out the constant as described. (The Raster Calculator will be REALLY picky about whether you use decimals and where you put spaces, so watch out! Use an expression as close as possible to this one.)
As a double check, you can compute the slope and aspect of your surface after you have generated it. Each should be approximately constant to about five decimal places and equal to what you intended.
The most common error made in ArcGIS is to compute the trig functions using radians instead of degrees, because ArcGIS uses sometimes one and sometimes the other. Because this is so quick and easy to do, first reproduce the example here so you know you have your angular units correct, then generate the plane you really want.
--Bill Huber
You can take the logic and work it up into a Python script. Some useful discussion and links are found from this ArcGIS Idea, according to which the $$XMAP were special variables not available after 9.3.1.
