[GRASS-dev] [GRASS GIS] #2359: r.stream.distance with a large map

[GRASS GIS]

#2359: r.stream.distance with a large map
-------------------------------+--------------------------------------------
 Reporter:  hcho               |       Owner:  grass-dev@…              
     Type:  defect             |      Status:  new                      
 Priority:  normal             |   Milestone:  7.1.0                    
Component:  Raster             |     Version:  svn-trunk                
 Keywords:  r.stream.distance  |    Platform:  All                      
      Cpu:  All                |  
-------------------------------+--------------------------------------------
 r.stream.distance creates FCELL output maps, but sometimes with a large
 map, the precision of FCELL is not enough. For example, to create a
 longest flow path map:
 {{{
 r.stream.distance -o stream_rast=outlet direction=drain method=upstream
 distance=flus
 r.stream.distance -o stream_rast=outlet direction=drain method=downstream
 distance=flds
 r.mapcalc expression="flusds=flus+flds"
 eval `r.info -r flusds`
 r.mapcalc expression="lfp=if(flusds>=$max-0.1, flusds, null())"
 }}}

 lfp should be a single connected stream representing the longest flow
 path, but because of rather big floating-point errors caused by FCELL, it
 can have more than one segment of streams.

 In line 434 in distance_calc.c:
 {{{
 cur_dist = tmp_dist + G_distance(...);
 }}}
 As tmp_dist becomes large, I get the following result:
 {{{
 cur_dist = tmp_dist + G_distance(...);
 double double_dist = (double)tmp_dist + G_distance(...);
 fprintf(stderr, "float: %f = %f + %f\n", cur_dist, tmp_dist,
 G_distance(...));
 fprintf(stderr, "double: %f = %f + %f\n", double_dist, tmp_dist,
 G_distance(...));

 Output:
 float: 117685.031250 = 117642.601562 + 42.426407
 double: 117685.027969 = 117642.601562 + 42.426407
 }}}

 These errors get accumulated and can cause unexpected problems. I suggest
 to change the default type of output maps to DCELL and add a flag (-f) for
 FCELL outputs, if needed.

-- 
Ticket URL: <http://trac.osgeo.org/grass/ticket/2359>
GRASS GIS <http://grass.osgeo.org>