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

ecodiv · January 22, 2014, 10:20pm

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?

Cheers

Paulo

lucadelu · January 23, 2014, 8:37am

Hi Paulo

On 22 January 2014 23:20, Paulo van Breugel <p.vanbreugel@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@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

--
ciao
Luca

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

ecodiv · January 23, 2014, 9:14am

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

···

On 01/23/2014 09:37 AM, Luca Delucchi wrote:

Hi Paulo

On 22 January 2014 23:20, Paulo van Breugel [<p.vanbreugel@gmail.com>](mailto:p.vanbreugel@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@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

mlennert · January 23, 2014, 10:15am

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@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 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

ecodiv · January 23, 2014, 10:28am

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@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 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.

Cheers,

Paulo

mlennert · January 23, 2014, 10:35am

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@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 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