Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

share|improve this question
add comment

1 Answer

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.

share|improve this answer
    
Is there a way for me to create an r.mapcalc script that ignores null values for the comparison? I think it will resolve the problem of null cells being taken into consideration. –  teamAngat Mar 1 '11 at 9:34
    
You could for example do a check on the resulting value after the neighbourhood calculation. If that is NaN, then put it to the original value. That would at least not eat up your raster, even if the border would be of inconsistent values. If you do the script on a wider enough area than the one of interest, that should not affect you. –  moovida Mar 1 '11 at 11:15
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.