[GRASS-user] r.terraflow

Tue Jul 29 08:21:49 EDT 2008

Thanks Thomas. I think we should put these on the wiki... When I find time.

I understand the coding scheme, but reclassifying from directions
r.terraflow scheme to the r.watershed one yealds significative
differences, and r.water.outlut works bad with the r.terraflow
results... The main problem is that r.terraflow produces lots of
zeroes, where r.watershed resolve the directions.

I think I will go on with r.watershed, but having to tile the whole
raster, I have to solve the problem of patching the tiles (because my
basins cover areas bigger then what my machine can manage in one
shot). The patches produce artifacts along the boundaries.

I will let you know if I find solutions.

2008/7/29 Thomas Adams <Thomas.Adams at noaa.gov>:
> Ivan et al,
>
> Below are some snippets from an email thread in 2005:
>
> Has anyone had any luck using the results from r.terraflow in
> r.water.outlet? More specifically I want to use the output direction grid
> from r.terraflow for the drainage direction map input of r.water.outlet. The
> r.water.outlet documentation only refers to using a drainage direction map
> generated in r.watershed.
>
> I really don't know about r.water.outlet but while r.whatershed generate
> flow direction with a D4 or D8 algorithm while r.terraflow use a Dinf
> algorithm (6 month ago the -s flag was not working, maybe now is fixed,
> maybe not.... try to test....).
>
> r.terraflow uses by default a MFD algorithm (multiple flow path). In short,
> the outflow of one cell is distributed among ALL its downslope neighbors,
> according to the relative differences in elevation. The nb stored in the
> flow direction output of terraflow corresponds to the neighbors id, stored
> as an octet (so, between 0 and 255), clockwise from East (1 for East, 2 for
> SE, 4 for S, 8 for SW...) : for example, if water flows to the North and the
> south simultaneously, the flow direction value is 68 (4 for S + 64 for N).
> NNE ? 64+128, and so on...
>
> r.watershed uses a SFD (single flow direction) algorithm, where water is
> routed to ONE downslope neighbor. The flow directions are coded as an
> integer between 1 and 8.
>
> There should be a trick to get main flow directions (pseudoSF) from
> terraflow MFD: check the sources, direction.cc: there's a series of test to
> get the main direction from the 0-255 value (function findDominant). A
> simple script should do the work. I tried to write such a script smoe time
> ago, but don't ask me where it is, my box is a mess these days...
>
> OK, that makes sense. If I use the -s flag to force r.terraflow to use SFD
> will my flow directions be 1-8 rather than 0-255?
>
> Nope, the values will still be between 0 and 128, but as
> 1,2,4,8,16,32,64,128 only: you'll still have to postprocess your data to get
> it to a number between 1 and 8. And make sure that the way the directions
> are reported in r.watershed are consistent (viz, 1 is E and you go
> clockwise)
>
> I did try r.terraflow initially but in this case r.fill.dir is much faster
> because all I need is flow direction, not flow accumulation, TCI or any of
> the other things r.terraflow provides. Also, I was never able to figure out
> how to get the flow direction r.water.outlet needs from the r.terraflow flow
> direction (I was the guy that started that thread you mentioned). Regardless
> of what I use to get flow direction it would be helpful to know what format
> r.water.outlet expects, ie what direction is 1 North, East or something
> else.
>
> for the record, here's some aspect map format information collected from
> the help pages and source code:
>
> (use d.rast.arrow and d.rast.num to investigate AGNPS, ANSWERS, GRASS
> formats)
>
>
> AGNPS (r.fill.dir)
> category values from 1-8, with 1 facing north and increasing values in
> the clockwise direction.
>
>
> ANSWERS (r.fill.dir)
> category values from 0-360 degrees, with 0 (360) facing east and values
> increasing in the counter clockwise direction at 45 degree increments.
>
>
> GRASS (r.fill.dir)
> broken in 6.0.0. (was old GRASS 4? format) In the old days (pre-floating
> point GRASS rasters) there was a 24-points based GRASS aspect map format
> which still rears its ugly head sometimes:
> 24,1,2   east
> 3,4,5    north east
> 6,7,8    north
> 9,10,11  north west
> 12,13,14 west
> 15,16,17 south west
> 18,19,20 south
> 21,22,23 south east
>
>
> GRASS (r.slope.aspect, others)
> categories represent the number degrees of east and they increase
> counterclockwise: 90deg is North, 180 is West, 270 is South and 360 is
> East. Beware: 0 may be null! Similar to the ANSWERS format AFAICT.
>
>
> r.param.scale
> the direction of maximum gradient (considered downslope) is stored as
> (West is 0 degree, East is ± 180 degree):
> * 0..+180 degree from West to North to East
> * 0..-180 degree from West to South to East
> [Anomaly: I vote for changing this to standard GRASS5 aspect format]
>
>
> r.watershed drainage map (what r.water.outlet wants)
> drainage direction. Provides the "aspect" for each cell. Multiplying
> positive values by 45 will give the direction in degrees that the
> surface runoff will travel from that cell. The value -1 indicates that
> the cell is a depression area (defined by the depression input map).
> Other negative values indicate that surface runoff is leaving the
> boundaries of the current geographic region. The absolute value of these
> negative cells indicates the direction of flow.
> [should we change this to standard AGNPS(+neg flag) aspect format ??]
>
>
> It would be nice to get these down to AGNPS, ANSWERS, and GRASS format
> anyway.
>
>
>
> so AGNPS is not the same as r.watershed format, but 45deg off:
>
> dir AGNPS  r.watershed
> N    1      0 or 8 ?
> NE   2      1
> E    3      2
> SE   4      3
> S    5      4
> SW   6      5
> W    7      6
> NW   8      7
>
>
> fix:
>
> r.mapcalc 'aspect.watershed = if(aspect.AGNPS == 1, 8, \
>  if(aspect.AGNPS > 1 && aspect.AGNPS <= 8, aspect.AGNPS-1, null() ))'
>
>
> if that works we can add it to the help page. (looks ok to me)
>
> note the negative values from the watershed program will be lost, which
> might adversely affect the output of r.water.outlet! (no idea)
>
> FYI, I've whipped up a Matlab script to extract the dominant flow direction
> from the r.terraflow MFD output map and reclass it into
> the GRASS aspect format for use with d.rast.arrow, etc. It should
> be easy enough to rewrite it in something more Free and convert the
> export into any of the above aspect formats by changing asp_mtx[].
> Usually you can just use r.slope.aspect directly as a proxy.
>
> find it here:
> http://grass.gdf-hannover.de/twiki/bin/view/GRASS/GrassAddOns
>
> AFAICT: 1 is NE, moving around the compass clockwise until 8 is North,
> and mind the negative values and information lost by taking the absolute
> value of flow direction.
>
> Drainage direction needed by r.water.outlet:
>
>
>
> N: this means that water is going to North ...
>
> N: 2
> N/E: 1
> E: 8
> S/E: 7
> S: 6
> S/W: 5
> W: 4
> N/W: 3
>
>
> I hope this helps…
>
> Regards,
> Tom
>
>
> ivan marchesini wrote:
>>
>> Il giorno mar, 29/07/2008 alle 02.22 +0200, G. Allegri ha scritto:
>>
>>>
>>> Thanks Luca.
>>> The coding mapping seems correct. So, zeroes are for unsolved pixels...
>>> I'll do some tests to compare the results with r.watershed. At now
>>> they seem quite different. I'm not sure if I should trust onflow
>>> directions from r.terraflow to use them in r.water.outlet...
>>>
>>>
>>
>> hmmmm.... in my opinion the r.terraflow output is not usefull for
>> r.water.outlet but if you find a solution please let we know   ;-)
>>
>>  bye
>>
>> Ivan
>>
>>
>>
>>
>>
>>>
>>> http://lists.osgeo.org/mailman/listinfo/grass-user
>>>
>>>
>
>
> --
> Thomas E Adams
> National Weather Service
> Ohio River Forecast Center
> 1901 South State Route 134
> Wilmington, OH 45177
>
> EMAIL:  thomas.adams at noaa.gov
>
> VOICE:  937-383-0528
> FAX:    937-383-0033
>
>