R (Raster): how to record ties using which.max

1
votes

I am new to using the R program. I have a stack of 10 rasters where each raster represents the probability (between 0 and 1) of a ecosystem type (e.g conifer forest, grassland, wetland) occurring at a given point (pixel) on the landscape. All of the raster layers have the same numbers of rows and columns.

I have been using the "which.max" command in the raster package to find out which raster layer (ecosystem type) has the highest probability of occurring at a point on the landscape. I save the resulting raster as an ascii raster file and look at it in ArcMap or SAGA to see how the 10 ecosystem types are distributed across the study area. The resulting raster identifies the ecosystem type at each pixel by identifying which layer in the raster stack has the highest probability so if the value is "1", then the first raster in the stack has the highest probability of occurring at that point on the landscape.

However, from what I understand, if two or more rasters have the same value, the "which.max" command identifies only one raster layer (the last raster with the tied value?). This means that the resulting raster could be strongly influenced by the order in which I assign the 10 ecosystem types in the raster stack.

What I would like to do is have a resulting raster record which layers tied if more than one raster layer has the same maximum probability. For example, I was thinking that if layers 4 and 8 in the raster stack have the same maximum value, then the resulting raster layer could identify this by recording the value "48". Or if three raster layers all had the same maximum value, the resulting value could be "259"

Does anyone have any suggestions on how to do this?

Thanks.

2
rr-raster

2 Answers

2
votes

try this:

> # create a small sample set 
> raster <- matrix(sample(3, 25, TRUE), 5)
> raster
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    1    2    3    2
[2,]    2    1    2    1    3
[3,]    2    1    1    2    1
[4,]    1    3    2    3    2
[5,]    1    3    3    1    3
> # create a new value representing the column number (assuming 1-9)
> cols <- apply(raster, 1, function(.row){
+     mCols <- which(.row == max(.row))
+     sum(mCols * 10 ^ (rev(seq_along(mCols)) - 1))
+ })
> cbind(raster, cols)
               cols
[1,] 3 1 2 3 2   14
[2,] 2 1 2 1 3    5
[3,] 2 1 1 2 1   14
[4,] 1 3 2 3 2   24
[5,] 1 3 3 1 3  235
0
votes

Here an example with a RasterBrick based on Data Munger's nice example

whiches <- function(i, stat=max) {
    m <- which(i == stat(i, na.rm=TRUE))
    sum(m * 10^(rev(seq_along(m)) - 1))
}

library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
x <- calc(b, whiches)

freq(x)
plot(x)

Note that you can only use this if you have less than 10 layers!