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 <gl...@gclements.plus.com>
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user
