The resulting maps from a r.mapcalc script is not computed properly

Question

This post is a continuation of my other post Why is my controlling shell script not working?. I created another post because this contains a separate issue that I encountered with the same scripts of the previous post.

I'm trying to do the Hydrologic Simulation seen in page 10 of this document.

Here is the script:

#!/bin/sh

g.region rast=elevation.dem -p

r.mapcalc "water = 120" #12 inches of water in each cell

g.copy rast=elevation.dem,elev

d.mon x0
d.rast water
i=1

while [ $i != 100 ]
do
  n=1
  while [ $n != 10 ]
  do
    r.mapcalc < water.mapcalc # run the simulation
    d.rast water
    n=`expr $n + 1`
  done
  g.copy rast=water,water.$i
  i=`expr $i + 1`
done

The following is the water.mapcalc script:

water = water + eval(x = elev + water, \
if (x > (y = elev[-1, 0] + water[-1, 0]), \
   -.15 * if (elev > y, water, x - y), \
   .15 * if (elev[-1, 0] > x, water[-1, 0], y - x))+ \
 if (x > (y = elev[1, 0] + water[1, 0]), \
   -.15 * if (elev > y, water, x - y), \
   .15 * if (elev[1, 0] > x, water[1, 0], y - x))+ \
 if (x > (y = elev[0, -1] + water[0, -1]), \
   -.15 * if (elev > y, water, x - y), \
   .15 * if (elev[0, -1] > x, water[0, -1], y - x))+ \
 if (x > (y = elev[0, 1] + water[0, 1]), \
   -.15 * if (elev > y, water, x - y), \
   .15 * if (elev[0, 1] > x, water[0, 1], y - x))+ \
 if (x > (y = elev[-1, 1] + water[-1, 1]), \
   -.10 * if (elev > y, water, x - y), \
   .10 * if (elev[-1, 1] > x, water[-1, 1], y - x))+ \
 if (x > (y = elev[1, 1] + water[1, 1]), \
   -.10 * if (elev > y, water, x - y), \
   .10 * if (elev[1, 1] > x, water[1, 1], y - x))+ \
 if (x > (y = elev[1, -1] + water[1, -1]), \
   -.10 * if (elev > y, water, x - y), \
   .10 * if (elev[1, -1] > x, water[1, -1], y - x))+ \
 if (x > (y = elev[-1, -1] + water[-1, -1]), \
   -.10 * if (elev > y, water, x - y), \
   .10 * if (elev[-1, -1] > x, water[-1, -1], y - x)))

and I used it with the Spearfish dataset.

The result of the code is that the map (elevation.dem) is getting smaller every succeeding iteration. I cannot figure out why it is happening. Isn't it that the result for every iteration (using the script above) should just be the transfer of excess height and not the cropping of the edges of the map?

Answer 1

How are the novalues borders handled? You are using a FCELL map type, right (in that case it would be the NaN)? The above script uses boundary cells, so if you have an extracted basin, most probably that might mean that at every iteration there is one more cell that has to deal with border novalues (which might be the reason for eating your map at every turn, calculations on novalues should produce novalues). The above example was done in 1992, where I think GRASS used to handle only CELL type maps, so the novalue were just zeros (not NaN). Also I think that the script was meant to be an example of the power of the mapcalc and not a sophisticated model to run.