Computing global DEM topographic indices and compromising projections?

Question

I would like to compute global topographic indices based on SRTM elevation data. These indices are computed along river/channel networks extracted from a DEM and are a functions of both local slope and upstream area obtained from a flow accumulation routine. For those that are familiar them, the indices I am interested in are a stream power and channel steepness index.

To compute these on a global grid (worldwide from ca. 60°N to 60°S) and to be able to compare slopes and areas in different regions of the world, I was thinking about using compromise projections such as a Dymaxion/Fuller map, which would limit both area and angle distortions. These would allow me to treat the earth as a whole and ease the processing. However I have a hard time getting a handle on the distortions these kind of projection induce and the associated uncertainty or bias in slope/area (I would be happy with less than a few percent).

Alternatives could be to split the map into small regions where UTM projections would provide good alternatives, but this would mean a more complicated workflow. Another alternative is to go for a conformal projection to get the slopes right and then run a flow accumulation routine that is weighted by a correction grid which would account for changes in area with changes in latitude.

Does anyone have comments or suggestions about these workflows?

As these global estimates are computationally intense I'd better think twice about the right projection to use.

Answer 1

A conformal projection on a world scale is a good idea, but it won't get the slopes right. Conformality assures that the scale distortion at any point does not vary with direction: that makes the slopes in all directions comparable, resulting in relatively accurate computation of aspect and of flow drainage directions. However, all conformal maps covering large areas must introduce substantial distortion from one point to the other, which means they cannot be used by themselves to compute a consistent slope across the entire study area.

One solution is to correct the computed slopes. This method is described in my answer at Global DEM to slope calculation.

Having done this, perform the followup calculation of flow directions, aspects, etc. using the corrected slopes. Then, for carrying out a flow accumulation calculation (or other calculations that depend on slopes and aspects), consider reprojecting those DEMs to an equal-area projection. There are many global equal-area projections to choose from; for instance, any cylindrical equal area projection will be fine for latitudes up to 60 degrees (or even almost to the poles).

---

One potential problem to anticipate concerns the resampling of values that occurs during reprojection. Most flow accumulation algorithms discretize the flow directions. For instance, ArcGIS uses a "D8" algorithm in which only eight directions are recorded, employing a binary code comprising powers of two from 1 = 2^0 through 128 = 2^7. Terrible things can happen to these values if they are averaged during reprojection. For instance, direction 127 is next to direction 1: during resampling just about any direction between 1 and 127 can result near those locations, which can create a huge mess. One also has to wonder what artifacts a nearest-neighbor resampling procedure (which doesn't average) would produce.

I have not experimented with the following solution and so can only propose it as an idea for further consideration and testing. As with most "circular" data, an effective approach to averaging is to replace directions by their cosines and sines (essentially converting a polar coordinate into Cartesian coordinates). This could be done for flow direction codes merely with a lookup table. For example, ESRI software codes would be transformed as follows:

Code Direction Bearing Cosine Sine
---- --------- ------- ------ ----
   1         E      90      0    1
   2        SE     135     -s    s
   4         S     180     -1    0
   8        SW     225     -s   -s
  16         W     270      0   -1
  32        NW     315      s   -s
  64         N       0      1    0
 128        NE      45      s    s

where s = Sqrt(1/2) = 0.7071... . (I would be tempted to multiply these values by the tangent of the slope, turning them into components of a vector field representing the gradient of the surface. Reprojection would then be averaging the gradient itself.)

After reprojection, the resulting (cosine, sine) pairs can be converted back into flow direction codes according to the same scheme. The inverse tangent of the interpolated (cosine, sine) determines an angle which can be written in degrees between 0 and 360. Dividing by 45 = (360/8) gives a number between 0 and 8. Rounding this to the nearest whole number determines the flow direction. For example, suppose the interpolated cosine equals -0.1 and the interpolated sine equals 0.9. (Notice their squares do not sum to unity! The separate interpolations will usually destroy this relationship.) The inverse tangent of this pair is 96 degrees which, when divided by 45, equals 2.14. That rounds to 2, corresponding to 2 * 45 = 90 degrees, for a direction code of 1.

The same approach will work with any flow direction coding scheme: convert the codes into (cosine, sine) pairs, reproject, then convert back after determining the angle via the inverse tangent.

---

For an assessment of the distortions in flow direction that arise when these issues are not addressed, please read the reply at Calculating flow direction and delineating basins....