On 01/23/2014 09:37 AM, Luca Delucchi wrote:
Hi Paulo
On 22 January 2014 23:20, Paulo van Breugel <[email protected]> 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@data, 1, function(x){
rank(x, na.last="keep", ties.method="average")
})
lyrs@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
_______________________________________________
grass-dev mailing list
[email protected]
http://lists.osgeo.org/mailman/listinfo/grass-dev
