Using SRTM Global DEM for Slope calculation?

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?

• This depends on where you are analyzing in the world. There are different projections for each location. Where are you examining? – djq Nov 6 '12 at 19:36
• Resolution is actually ~90m, not ~90km. – Akheloes Aug 22 '13 at 10:09
• Just a comment, if your on maintenance for Desktop, you can login to ArcGIS Online and use their elevation services (with no need for NA extension). The slope layer is free to use as a reference layer. In Australia, we have the 1 second SRTM data (~30m res) blogs.esri.com/esri/arcgis/2014/07/11/… – Simon May 8 '15 at 12:39

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.

• I find myself in the position, as I so often do, of thanking whuber once again for describing a solution as well as putting forth the reasoning that builds it. My hat is off to you sir. – matt wilkie Nov 7 '12 at 8:35
• Thanks @matt. I did not mean to suggest earlier that your (now deleted) answer should be rescinded: In fact, I had upvoted it because you shared a link to an interesting USGS reference that could be of use to many readers. (My comment was critical only of a secondary passage in that paper, not the paper itself.) – whuber Nov 7 '12 at 16:14
• ahh. thanks for the clarification. I've restored the answer, trusting people have enough information in front of them now to make an informed choice :) – matt wilkie Nov 9 '12 at 5:28
• Coming from a french background, it took me a while to translate the necessary terminology to better understand this great answer, so I thought dropping this link is a good help for newbies such as myself : webhelp.esri.com/arcgisdesktop/9.2/… – Akheloes Aug 22 '13 at 16:20
• It's a great approach and I already used your solution to generate a global slope raster. One hint from practical experience: Since latitude values south of equator are negative, you have to use the absolute latitude value in the following equation: (1 + sin(latitude))/2 – Saleika Jun 23 '17 at 11:46

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

Better answer, using GDAL with a scale

Now that SRTM data are globally available I can actually see and work with the files. The gdaldem utility from GDAL can compute slope and hillshade using a scale option for a ratio of vertical units to horizontal. The manual recommends 111120 m/° for something like SRTM tiles. So for example, from an OSGeo4W shell:

\$ gdaldem slope -s 111120 -compute_edges N44E007.hgt N44E007_slope.tif

The -compute_edges option makes the edges more seamless, if you want to stitch a few tiles together. Or compute tiles for a large region. The disadvantage with the "scale" technique is that distances in the E-W and N-S directions are not equal, except at the equator, so for tiles closer to the poles, there could be some odd misrepresentations of slope.

• It's worth emphasizing your final comment: this is a poor solution for points not near the Equator. This not a small matter of "odd misrepresentations": the results will be grossly wrong, especially in places closer to the Poles than the Equator. The documentation for gdaldem states "For locations not near the equator, it would be best to reproject your grid using gdalwarp before using gdaldem." Unfortunately that won't work for datasets covering the globe, unless you break them into small pieces (74 UTM zones, perhaps?), project those, compute the slopes, and mosaic the results. – whuber Jan 17 '18 at 21:24

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 (, , ) 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. :)

• That paper's suggestion that an azimuthal equidistant projection could be used to compute slopes is misguided and wrong. It will indeed give correct slopes near the origin of the projection, but they, too, will get progressively less accurate as distance to the origin increases. – whuber Nov 6 '12 at 21:00
• thanks for pointing that. Readers, please be sure to read gis.stackexchange.com/a/40464/108 as well, for balance – matt wilkie Nov 9 '12 at 5:29

Global DEM parameters (where most of formulas are based on assumption of Euclidean space) can be derived efficiently using the EQUI7 GRID system (Bauer-Marschallinger et al. 2014). EQUI7 GRID splits world into 7 land areas, all projected in equidistant projection system with a minimum loss of precision. See en example of global DEM at 250 m resolution in the EQUI7 GRID. Here you can find some sample code that shows how to derive global DEM parameters using SAGA GIS. Once you complete deriving DEM parameters in the EQUI7 GRID system, you can back-transform all maps into WGS84 longlat coordinates and then create a global mosaick using GDAL.

• Could you explain how this answers the question? If you're proposing using equidistant projections for slope calculations, please be aware that's a poor solution because of the large relative metric distortions that occur as one moves out from the center of the projection. Although focusing seven such projections on land masses helps alleviate that problem, it's still not the best choice. – whuber Feb 20 '17 at 17:33
• The paper by Bauer-Marschallinger et al. (2014) explains why were these projections chosen to represent global land masses (they are assumed to have minimum loss in precision). I agree that any 2D projection will eventually lead to deformations, but as far as I know, EQUI7 is a good compromise between loss in precision and comodity (2D algebra). Having said that, hexagons are again being used to represent global land surfaces (although DEM analysis with 3D hexagons is still a cumbersome). – Tom Hengl Feb 21 '17 at 19:54
• Thank you for the reference. Its abstract suggests it solves a quite different problem, that of "minimizing local data oversampling appearing during projection of generic satellite images to a regular raster grid." That does not imply it will perform well for other purposes, such as estimating slopes. – whuber Feb 21 '17 at 20:07
• Of course EQUI7 does not absolutely solve the problem of estimating the local slopes accurately, but it is probably a more elegant solution than using the above suggested Mercator projection. Ultimately, if one wishes to estimate slopes with perfect precision, then the only options are probably to (1) use local (equidistant) projections per smaller size tiles (e.g. 100 by 100 km) with 10-20% overlap, as mentioned also in Verdin et al. (2007), or (2) to use the hexagon grid (dggridR package). – Tom Hengl Feb 24 '17 at 9:57
• The problem isn't precision--it lies in producing systematically biased slopes and aspects. Because the equidistant projections differentially distort directions orthogonal to the geodesics originating at their centers, the aspects will always be wrong (although reasonably accurate close to the centers where all distortion is low) and the errors in the slopes will grow rapidly with distance. The use of many local projections of course will work, but is the very opposite of the elegance you value. – whuber Feb 24 '17 at 14:12

Slope is rise / run. Compute rise and compute run and you have your answer. It is simple to compute the distance between geographic coordinates. This will introduce less resampling error as compared to conversion to UTM, etc.