[GRASS-dev] Re: comparing r.cost and r.terracost [was: Re: [GRASS-user] affecting cat of closest feature as pixel value]

[Moritz Lennert]

Hello Laura,

On 15/07/09 00:01, Laura Toma wrote:
> Hi Markus, Moritz,
> 
> I could reproduce the problem. It was happening because the cost of the 
> source points were incorrectly assigned to be zero, instead of the 
> actual cost of the points.  Because of  this, the immediate neighbors of 
> the source points had "half" costs.  I fixed one line in initialize.cc 
> and checked in the change in grass-addons/raster/r.terracost   
> (hopefully it is  the right place).

yes, perfect. Thanks for fixing this so quickly. Now results are almost 
identicial to r.cost without the knight's move. Differences are very 
small, and (my guess) probably due to differences in rounding:

g.region n=10 s=0 w=0 e=10 res=1 -a
r.mapcalc one=1.0
echo "5|5" | v.in.ascii out=start --o
v.to.rast start use=val val=1 out=start --o
time r.cost one start_rast=start out=r_cost --o
time r.terracost one start_rast=start out=terra_cost --o
r.mapcalc diff3=r_cost-terra_cost

r.univar diff

n: 100
minimum: -1.90012e-07
maximum: 1.91819e-07
range: 3.81831e-07
mean: -4.1721e-09

And image of diff:
http://geog-pc40.ulb.ac.be/grass/misc/comparison_rcost_rterracost_4.png

This seems to indicate that r.cost "overestimates" diagonals and 
"underestimates" the straighter lines when compared to r.terracost.

The same, but with a 100x100 region:

r.univar diff:
n: 10000
minimum: -3.04995e-05
maximum: 9.16214e-06
range: 3.96616e-05
mean: 6.50277e-07

And image of diff:
http://geog-pc40.ulb.ac.be/grass/misc/comparison_rcost_rterracost_3.png

This seems to confirm the above, but not for farther distances...


Second test:

g.region rast=elev_state_500m
v.to.rast nc_state use=val val=500.0 out=cost --o
v.to.rast urbanarea use=val val=1 out=urbanarea --o
r.cost in=cost start_rast=urbanarea out=dist_urban --o
r.terracost in=cost start_rast=urbanarea out=dist_urban_terracost --o
r.mapcalc diff=dist_urban-dist_urban_terracost

r.univar diff

n: 513065
minimum: -0.0125342
maximum: 0.0684331
range: 0.0809673
mean: 0.00173129

Third test:

g.region rast=landclass96
r.reclass landclass96 out=cost_landcover --o << EOF
1 = 200
2 = 20
3 = 5
4 = 50
5 = 30
6 = 1000
7 = 5
EOF
echo "629980|228527.25" | v.in.ascii out=point --o
v.to.rast point use=val val=1 out=point --o
r.cost in=cost_landcover out=dist_rcost start_rast=point --o
r.terracost in=cost_landcover out=dist_terracost start_rast=point --o
r.mapcalc diff2=dist_rcost-dist_terracost

r.univar diff2

n: 249323
minimum: -0.554688
maximum: 0.523438
range: 1.07812
mean: -0.0104599

In all cases r.terracost takes about 1/4-1/3 the time of r.cost.

In the idea of replacing r.cost by r.terracost, I would think that the 
latter would first have to implement the options in r.cost, e.g.

- equivalents of the start_points and coordinates options of r.cost
- stop_points/coordinates
- max_cost
- the knight's move
- ...

Moritz