SRTM Global DEM Slope – Using SRTM Global DEM for Slope Calculation

arcgis-10.0arcgis-desktopdemsrtm

I have downloaded SRTM GDEM (~90m resolution).

I am using ArcGIS 10.

I have tried to use the Spatial Analyst tool to compute for slope.

However, I cannot compute for the slope.

The output values has only two ranges 0 and 0.1-90.

I am not really sure what the problem is.

Best Answer

This seems like a good place to describe a simple, fast, and more than reasonably accurate way to compute slopes for a globally extensive DEM.

Principles

Recall that the slope of a surface at a point is essentially the largest ratio of "rise" to "run" encountered at all possible bearings from that point. The issue is that when a projection has scale distortion, the values of "run" will be incorrectly computed. Even worse, when the scale distortion varies with bearing--which is the case with all projections that are not conformal--how the slope varies with bearing will be incorrectly estimated, preventing accurate identification of the maximum rise:run ratio (and skewing the calculation of the aspect).

We can solve this by using a conformal projection to ensure that the scale distortion does not vary with bearing, and then correcting the slope estimates to account for the scale distortion (which varies from point to point throughout the map). The trick is to use a global conformal projection that allows a simple expression for its scale distortion.

The Mercator projection fits the bill: assuming scale is correct at the Equator, its distortion equals the secant of the latitude. That is, distances on the map appear to be multiplied by the secant. This causes any slope calculation to compute rise:(sec(f)*run) (which is a ratio), where f is the latitude. To correct this, we need to multiply the computed slopes by sec(f); or, equivalently, divide them by cos(f). This gives us the simple recipe:

Compute the slope (as rise:run or a percent) using a Mercator projection, then divide the result by the cosine of the latitude.

Workflow

To do this with a grid given in decimal degrees (such as an SRTM DEM), perform the following steps:

  1. Create a latitude grid. (This is just the y-coordinate grid.)

  2. Compute its cosine.

  3. Project both the DEM and the cosine of the latitude using a Mercator projection in which scale is true at the Equator.

  4. If necessary, convert the elevation units to agree with the units of the projected coordinates (usually meters).

  5. Compute the slope of the projected DEM either as a pure slope or a percent (not as an angle).

  6. Divide this slope by the projected cosine(latitude) grid.

  7. If desired, reproject the slope grid to any other coordinate system for further analysis or mapping.

The errors in the slope calculations will be up to 0.3% (because this procedure uses a spherical earth model rather than an ellipsoidal one, which is flattened by 0.3%). That error is substantially smaller than other errors that go into slope calculations and so can be neglected.

Fully global calculations

The Mercator projection cannot handle either pole. For work in polar regions, consider using a polar Stereographic projection with true scale at the pole. The scale distortion equals 2 / (1 + sin(f)). Use this expression in place of sec(f) in the workflow. Specifically, instead of computing a cosine(latitude) grid, compute a grid whose values are (1 + sin(latitude))/2 (edit: use -latitude for the South Pole, as discussed in the comments). Then proceed exactly as before.

For a complete global solution, consider breaking the terrestrial grid into three parts--one around each pole and one around the equator--, performing a slope calculation separately in each part using a suitable projection, and mosaicing the results. A reasonable place to split the globe is along circles of latitude at latitudes of 2*ArcTan(1/3), which is about 37 degrees, because at these latitudes the Mercator and Stereographic correction factors are equal to each other (having a common value of 5/4) and it would be nice to minimize the sizes of the corrections made. As a check of the computations, the grids should be in very close agreement where they overlap (tiny amounts of floating point imprecision and differences due to resampling of the projected grids ought to be the only sources of discrepancies).

References

John P. Snyder, Map Projections--A Working Manual. USGS Professional Paper 1395, 1987.

Related Question