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

[Paulo van Breugel]

On 01/23/2014 09:37 AM, Luca Delucchi wrote:
> 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
>>
Hi Luca,

Thanks..  I don't know much about Python unfortunately, so it is 
difficult for me to see exactly what it does. I'll try to see if I get 
this to work, and see if I understand the output.

Cheers

Paulo



-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20140123/ea8fd56b/attachment.html>