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.
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.
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.
To do this with a grid given in decimal degrees (such as an SRTM DEM), perform the following steps:
Create a latitude grid. (This is just the y-coordinate grid.)
Compute its cosine.
Project both the DEM and the cosine of the latitude using a Mercator projection in which scale is true at the Equator.
If necessary, convert the elevation units to agree with the units of the projected coordinates (usually meters).
Compute the slope of the projected DEM either as a pure slope or a percent (not as an angle).
Divide this slope by the projected cosine(latitude) grid.
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.
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).
John P. Snyder, Map Projections--A Working Manual. USGS Professional Paper 1395, 1987.
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).
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.
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.
-s 111120 -compute_edges and SHAZAM!
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. :)
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.
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.
