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

[Luca Delucchi]

Hi Paulo

On 22 January 2014 23:20, Paulo van Breugel <p.vanbreugel at gmail.com> 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 if the result is what you are looking for but with python
I do something like this

from grass.pygrass.raster import RasterNumpy
import scipy.stats as sstats
map = RasterNumpy('YOURMAP')
map.open()
newmap =  RasterNumpy('YOURNEWMAP', 'w')
newmap.open()
i = 0
for row in map:
    newmap[i] = sstats.rankdata(row)
    i+=1
newmap.close()
map.close()

I tested it, but when I close the map it return the error "Bus error"

> Cheers
>
> Paulo
>

-- 
ciao
Luca

http://gis.cri.fmach.it/delucchi/
www.lucadelu.org