Opened 5 years ago

Closed 5 years ago

Last modified 3 years ago

#2359 closed defect (fixed)

r.stream.distance with a large map

Reported by: hcho Owned by: grass-dev@…
Priority: normal Milestone: 7.2.0
Component: Raster Version: svn-trunk
Keywords: r.stream.distance Cc:
CPU: All Platform: All

Description

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.

Change History (3)

comment:1 Changed 5 years ago by hcho

I took the liberty to change the type of output maps from FCELL to DCELL for a higher precision (r61206). The above code for calculating the longest flow path should work now.

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())"

comment:2 Changed 5 years ago by neteler

Resolution: fixed
Status: newclosed

(r.stream.distance was moved back to Addons)

Closing as fixed.

comment:3 Changed 3 years ago by neteler

Milestone: 7.1.07.2.0

Milestone renamed

Note: See TracTickets for help on using tickets.