3

I am calculating Topographic Wetness Index for a study area using SRTM DEM (90 and 30m).

I've seen several versions of TWI (CTI). I have been using: TWI = ln(As +1 / tanB) where As = catchment area, B = slope (degrees/radians??)

When I run this function, the output contains negative values (using slope in degrees and radians). Are negative values acceptable? Or does this indicate an error in the analysis??

I understand the output should be scaled to 1-10, but negative values seems counter intuitive. Could someone explain what negative values could mean (due to flat areas, no data values?)

On a side not, the calculations are being run on on a very large area, which overaps several watersheds. Is it appropriate to run TWI over such a large extent?

3
  • 1
    Where does this formula for TWI come from? As stated, there's nothing to prevent it being negative: all you need is to measure catchment area in (say) millions of square kilometers, so that As is tiny, and have a slope appreciably greater than 45 degrees (so that 1/tan(B) is less than 1). At a minimum, then, the formula has to assume some specific units of measurement for the catchment area and that the the area is not too small and the slope is not too steep. (BTW, if you use degrees where radians are intended, you are likely to get many undefined values for the TWI.)
    – whuber
    Commented Dec 5, 2012 at 21:19
  • The original reference for the CTI is: Moore, ID., P.E. Gessler, G.A. Nielsen, and G.A. Petersen (1993) Terrain attributes: estimation methods and scale effects. In Modeling Change in Environmental Systems, edited by A.J. Jakeman M.B. Beck and M. McAleer Wiley , London, pp. 189 - 214. Another common one that is in a journal: Gessler, P.E., I.D. Moore, N.J. McKenzie, and P.J. Ryan. (1995). Soil-landscape modeling and spatial prediction of soil attributes. International Journal of GIS. 9(4):421-432. Commented Dec 6, 2012 at 18:59
  • 2
    I just noticed that this question was asked again recently here gis.stackexchange.com/questions/104543/… . I provided a detailed explanation for why this occurs at that site as well. Perhaps there is a way of linking the two questions? Commented Sep 11, 2014 at 16:40

1 Answer 1

10

You should not be seeing negative values in the CTI. Since you did not provide a reproducible example I cannot speculate as to why you are getting incorrect results. The expected range is not limited 1-10. The range will be defined by flow accumulation which is influenced by the size of the basins that are accumulating flow. The index does not rely on washed boundaries so the extent of the DEM should not matter.

The common implementation of CTI can be shown as: CTI = ln(a / (tan (beta)) where; a=Upland contributing area [(flow accumulation + 1 ) * (cellsize)]and beta=slope in radians.

Here are the raster algebra steps to calculate CTI (you will have to translate syntax to your version of ArcGIS):

cellsize=30
fd = flowdirection(dem)
sca = flowaccumulation(fd)
slope = ( slope(dem) * 1.570796 ) / 90
tan_slp = con( slope > 0, tan(slope), 0.001 )
sca_scaled = ( sca + 1 ) * cellsize
cti = ln ( sca_scaled / tan_slp ) 
9
  • This looks good overall, but I am puzzled by the formula you quote: if I measure cellsize in kilometers, say (or worse, in decimal degrees), then a / tan(beta) can easily be less than 1, leading to a negative CTI. Perhaps your formula assumes cellsize is given in some standard, fixed units such as meters?
    – whuber
    Commented Dec 5, 2012 at 21:22
  • Sorry, I should have specified that the assumed projection units are meters or feet. Commented Dec 5, 2012 at 23:30
  • Thanks for the reply, everyone. Jeffrey, I found the Geomorphometry Arc toolbox online and some of your other blogs. I will investigate further today. Thanks again.
    – derelict
    Commented Dec 6, 2012 at 18:18
  • One more question regarding the ArcToolbox! In the python script, flow accumulation is as follows: "corr_sca = ( sca +1 ) * 30 " but in the help file it states cell size should be in sq. meters: As = (flow accumulation + 1) * (pixel area in m2) So, should cell size be area, such as square meters (30x30 m), or spatial resolution (30m)? Also, will the function automatically recognize different cell sizes, or should the cell size be changed in the script before running?
    – derelict
    Commented Dec 6, 2012 at 19:10
  • 2
    The units of upslope contributing area are m2/m or m (in fact A is As or specific contributing area, contributing area per unit contour width) So that's why you multiply by cell resolution = 30 m and not by cell area (30*30) A more correct way to calculate this would be to take into account the flow direction (if it is diagonal or perpendicular), but for most applications, the approximated calculus given here is enough
    – user23115
    Commented Oct 21, 2013 at 10:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.