[GRASS-stats] r.mapcalc, median and NULL

[Hamish]

Agustin Lobo wrote:
> The grass doc is very clear on the behavior of r.mpacalc regarding NULL:
> "NULL-values in function arguments should result in NULL"
> 
> I have to calculate the median of a set of 12 monthly layers, i.e.
> mean(l1,l2,l3,...l12)
> 
> the problem is that some layers have some NULL values. I can
> replace NULL by, for example, 0, but this will be having a large effect
> in the value of the mean. Is not there any way of setting something
> like the na.rm=T of R?

in matlab this is called nanmean(). they (& grass) don't let the default
mean() do nice tricks like that automatically because it lets you silently
fool yourself*. If "n" is not consistent across the mean-map the result
is somewhat suspect. at minimum you'd also want a std dev map.
* (but excel does let you do that ;)

An idea from a smart friend of mine to fill in those gaps: if data is
multi-year you can make a null-filled average for each month without
introducing the temporal bias, then take a propagating-null average
for the full year from those filled-monthly maps.

see also the sea surface temperature hole-filling and 28-day rolling
composite ideas here:
  http://grass.osgeo.org/wiki/MODIS#Removing_holes


> Alternatively, could I calculate (in GRASS) a raster of the number of
> months with NULL value present for each i, j position?

sure, you could do that with r.mapcalc's if( isnull() ).

> In that way I could replace NULL by 0 but then add up the layers and
> divide by a raster that would be 12 in most cells, 11 in those having
> 1 NULL, 10 in those having 2 NULL etc

yeah, but consider why those cells are NULL. Maybe this is from satellite
surface temperature data and an area is often cloudy some time of year
(eg monsoon) so satellite can't see the ground and no data. e.g. if winter
months are often cloudy you get lots of NULLs in your winter month maps,
and by skipping those you are biasing the end result to be warmer than it
should be.


Hamish