[GRASS-user] r.mapcalc & neighborhood modifiers

[Ken Mankoff]

Hi List,

I'd like to calculate X (for example, change in elevation) between each raster cell and its downstream neighbor, where downstream is defined by r.watershed flow direction. I think I have achieved this with the following:

r.watershed elevation=z_s direction=fdir
r.mapcalc "d_z_s = z_s * 0"  # initialize delta z_s
for dir in `jot 8`; do 
  if [[ ${dir} == '1' ]]; then r=-1; c=1;  fi # 45 CCW from E; /°
  if [[ ${dir} == '2' ]]; then r=-1; c=0;  fi # 90 |^
  if [[ ${dir} == '3' ]]; then r=-1; c=-1; fi # 135 °\
  if [[ ${dir} == '4' ]]; then r=0;  c=-1; fi # 180 <-
  if [[ ${dir} == '5' ]]; then r=1;  c=-1; fi # 225 ./
  if [[ ${dir} == '6' ]]; then r=1;  c=0;  fi # 270 .|
  if [[ ${dir} == '7' ]]; then r=1;  c=1;  fi # 315 \.
  if [[ ${dir} == '8' ]]; then r=0;  c=1;  fi # 360 ->

  r.mapcalc "d_z_s = if(fdir == ${dir}, z_s[${r},${c}] - z_s, d_z_s)" --o
done

Now, I'd like to do something similar, which is route all stream water 1 grid cell downstream. I can't figure out how to do this. Perhaps this is beyond the ability of r.mapcalc, or reasonable bash scripting? Hopefully someone here can suggest a way to do this in grass.

I think the limitation is that I can't use neighborhood modifiers on the LHS of an r.mapcalc equation like this:

r.mapcalc "down[${r},${c}] = if(fdir == ${dir}, up, down)"

Which means I need to know the upstream cell, which is information I do not have. Looping through and finding the upstream cell(s) is non-trivial because there are likely to be multiple, whereas I know there is only 1 downstream cell because I ran r.watershed with "-s".

I think an alternate method would be to "accumulate" everything 1 cell downstream with r.watershed, and then subtract that raster from the original. But I have not been able to get the max_slope_length option of r.watershed to work.

Any advice will be greatly appreciated.

Thanks,

  -k.