Global DEM to Slope calculation

Question

I have downloaded SRTM GDEM (~90 km resolution). I am using ArcGIS 10. I have tried to use spatial analyst to compute for slope. However, I can not compute for the slope. The output values has only two ranges 0 and 0.1-90. I am not really sure what is the problem? Thanks.

Answer 1

I'm guessing the horizontal units for your raster are in either degrees or arcseconds. You need to reproject this raster to a spatial projection where your horizontal and vertical units are the same (i.e., if the vertical units are in metres, then I suggest using UTM, which has horizontal units of metres).

To reproject a raster with ArcCatalog/ArcGIS, look in:

ArcToolbox > Data Management Tools > Projections and Transformations > Raster > Project Raster

Choose a projected spatial reference that covers your region of interest, e.g., try a UTM zone. There are many other options which are best documented in the manual. Note, you cannot create a slope dataset for the whole Earth (if that is what you are trying to do).

Answer 2

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. 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.

Answer 3

Simply put, there isn't one. By definition a coordinate system based on degrees is un-projected. In common parlance we say WGS84 is a "geographic" projection, but that's untrue, just for convenience.

I think I remember reading about a software or process for accurately working with elevation models in un-projected geographic space but I can't locate it right now. In any case it would have been an experimental or build it yourself from code kind of process.

---

Ahhh, found it: Development of a Global Slope Dataset for Estimation of Landslide Occurrence Resulting from Earthquakes (USGS). Page 4 describes the problem well

...the length of one degree varies depending on its latitudinal location. At the equator, a one-degree by one-degree block is reasonably square when converted to units of meters (111,321 meters in the x-direction by 110,567 meters in the y-direction ... but closer to the poles the distances in the x-direction grow smaller as a function of the cosine of latitude, owing to convergence of the meridians. Most GIS packages, ArcGIS included, operate only on square pixels, and so using a factor to adjust the x, y, or z dimensions to a common unit is not possible.

The paper goes on to describe the specific calculations and software tools (gdal, python, numpy) they used to workaround this fundamental issue. The paper doesn't include the code, but if asked nicely they might share. In any case though I'd probably just ask where the results are, being the USGS it's probably already online somewhere. :)