[GRASS-user] average calculation

Glynn Clements glynn at gclements.plus.com
Thu May 14 20:05:11 EDT 2009

```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>
```