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@gclements.plus.com>