[GRASS-user] average calculation

[Glynn Clements]

Seb wrote:

> If one were to calculate an average of several rasters, one could simply
> do:
> 
> r.mapcalc "ave = (A + B + C) / 3"
> 
> But how can we get around the problem of null values in any of the
> rasters, which would propagate it to the result?  What is an efficient
> way to calculate both the numerator and denominator for each pixel so
> that it corresponds only to rasters with non-null values.

As Hamish says, "r.series method=average ...".

For r.mapcalc:

	r.mapcalc <<-EOF
		ave = eval(\
		n_A = isnull(A), n_B = isnull(B), n_C = isnull(C),\
		((n_A ? 0 : A) + (n_B ? 0 : B) + (n_C ? 0 : C))\
		 / (!n_A + !n_B + !n_C))
	EOF

>  A second
> problem is how to script this (shell) so that a large number of rasters
> can be included in this calculation.

	maps="map1 map2 map3 map4"
	(
	  echo "ave = eval(\\"
	  for map in $maps ; do
	    echo "n_$map = isnull($map),\\"
	  done
	  echo "(0\\"
	  for map in $maps ; do
	    echo " + (n_$map ? 0 : $map)\\"
	  done
	  echo ") / (0\\"
	  for map in $maps ; do
	    echo " + !n_$map\\"
	  done
	  echo "))"
	) | r.mapcalc

Using a real language would produce more legible code; e.g. in Python:

	import subprocess
	maps = ['map1', 'map2', 'map3', 'map4']
	nulls = ["n_%s = isnull(%s)" % (m, m) for m in maps]
	numer = ["(n_%s ? 0 : %s)" % (m, m) for m in maps]
	denom = ["!%s" % m for m in maps]
	expr = "ave = eval(%s,(%s) / (%s))" % (",".join(nulls), " + ".join(numer), " + ".join(denom))
	subprocess.call(['r.mapcalc', expr])

-- 
Glynn Clements <glynn at gclements.plus.com>