[GRASS-dev] Determine raster rank order per raster cell

[Moritz Lennert]

On 23/01/14 11:28, Paulo van Breugel wrote:
> On Thu 23 Jan 2014 11:15:17 AM CET, Moritz Lennert wrote:
>> On 22/01/14 23:20, Paulo van Breugel wrote:
>>> Does anybody know a smart way to calculate for X rasters per raster cell
>>> the rank order of those rasters? For example, if I have three rasters
>>> X1, X2 and X3:
>>>
>>> X1   X2   X3
>>> 1      3     2
>>> 2      5     8
>>> 5     1      3
>>> NA  8      2
>>>
>>> would give three new raster layers with the values:
>>>
>>> Y1   Y2   Y3
>>> 1      3     2
>>> 1      2     3
>>> 3      1     2
>>> NA   2    1
>>>
>>> I am now reading the rasters in R and using a combination of apply() and
>>> rank() function to create new raster layers:
>>>
>>> in <- c("X1",  "X2", "X3")
>>> lyrs <- readRAST6(in)
>>> tmp <- apply(lyrs at data, 1, function(x){
>>>      rank(x, na.last="keep", ties.method="average")
>>> })
>>> lyrs at data <- t(tmp)
>>>
>>> I am sure there are better ways in R, but especially when dealing with
>>> very large raster layers, I was hoping there is a way to do this
>>> directly in GRASS GIS. Or otherwise, would it be difficult to create a
>>> function for this?
>>
>> I don't know how efficient this would be, but how about using
>> r.mapcalc with something like this for each map:
>>
>> Y1 = if(X1>X2, if(X1>X3, 1, 2), if(X1>X3, 2, 3))
>>
>> Moritz
>
>
> Thanks Moritz,
>
> I had thought about that, works good for few maps, but it gets a bit
> unwieldy with e.g., 40 maps. I'll try though if I can split the if
> function somehow.

Maybe you can find a way to script the construction of the r.mapcalc calls.

You could then launch each r.mapcalc call as a separate process to allow 
"parallel processing" (see i.pansharpen and others as examples).

Moritz