r.topidx

r.topidx is a GRASS GIS module that generates the topographic index map used for modeling TOPMODEL.

1   Use the correct resolution

Re: topidx calculation

Since r.topidx tries to calculate topographic index values based on the resolution of the DEM, the DEM has to have the right resolution. If the resolution is greater than the actual grid size of data, r.topidx cannot know whether cells within a grid are sinks or not.

When your elevation data is (assuming its resolution is 3 m/grid):

123
456
789

In this map, each cell composes one “data grid.”

If the above data has a resolution greater than 3 m/grid (let’s say 1 m/grid here), an imported raster will look like the following:

111222333
111222333
111222333
444555666
444555666
444555666
777888999
777888999
777888999

In this map, one-cell data grids have been imported as nine-cell data grids, and nine cells compose one data grid.

Now, GRASS cannot say whether those duplicated (because of the greater resolution) nine cells per each data grid really compose one data grid or not. Similarly, r.topidx treats those cells as sinks.

To find the correct resolution of your DEM data, run the following script (find_res.sh dem):

#!/bin/sh
r.stats -1n "$1" | awk '{
	v = $0
	if(NR > 1 && v != v0){
		if(min_n == 0 || n < min_n)
			min_n = n
		n = 1
		v0 = v
		next
	}
	v0 = v
	n++
}
END{
	print min_n
}'

The above script counts the minimum number of duplicate cell values. Now you want to multiply the current resolution by the printed number. For example,

$ g.region -p | grep nsres | awk '{print $2}'
  1
$ find_res.sh dem
  3
$ g.region res=3	# 1 * 3
$ d.erase
$ r.topidx input=dem output=topidx

2   Be careful with flat areas

Re: TOPIDX flow direction algorithm

r.topidx has a known issue that flat areas have topographic index values of NULL since contributing areas cannot be calculated with the D8 algorithm when area is flat. If you know any way to work around this shortcoming, please let me know.