Hello, I'd really appreciate some help, 

I have four raster layers, each one represents transition 
potential for a particular land use, e.g. 1 urban, 2 irrigated, 3 
non-irrigated, 4 forest. In theoretical terms, the transition potential 
layers represent a competition for space between the 4 land uses, with 
the eventual land use being determined by "the highest bidder", that's 
to say, the land use that has the highest transition potential for that 
cell.   

So, to allocate the final land use, I need to create a new layer
 from the raster stack, in which each cell contains the number of the 
layer in the stack that had the highest value for that cell; e.g. if 
layer r1 had the highest value for a particular cell out of all the 
layers in the stack, that cell should have value 1 in the new layer. So 
far, simple enough. But this is not quite what I need, because I need to
 allocate a fixed quantity of cells (the "demand") from each layer, not 
all the cells that won the competition. This is where it gets slightly 
more complicated. If the cell demand (e.g. 5) for a particular layer, 
e.g. r1, is less than the number of cells that won the competition (e.g.
 21), then the new map needs to contain 1s only where the 5 highest 
value cells from map 1; once these 5 cells have been allocated, no more 
1's can be allocated, so the cells that would have taken value 1 should 
take the value from the layer that came second instead. This is as far as I've got:

#start R script 
require(raster) 
r1 <- raster(nrows=10, ncols=10) 
r2=r3=r4=r1 
r1[]= runif(ncell(r1)) 
r2[]= runif(ncell(r1))+0.2 
r3[]= runif(ncell(r1))-0.2 
r4[]= runif(ncell(r1)) 
rs=stack(r1,r2,r3,r4) 

r1demand <- 5  #this is the number of cells I want to allocate from r1 
r2demand <- 8	#this is the number of cells I want to allocate from r2 
r3demand <- 14	#this is the number of cells I want to allocate from r3 
r4demand <- 7	#this is the number of cells I want to allocate from r4 

miwhichmax <- function(x,na.rm=TRUE) which.max(x) 
rsmaxvals <- stackApply(rs,indices=c(1,1,1,1),max)  #this returns
 a single raster with the highest values from all of the other layers 
allocated <- stackApply(rs,indices=c(1,1,1,1),miwhichmax) #this 
allocates a value to each cell corresponding to the layer which has the 
highest value for that cell. 


Any ideas?

Many thanks, 
Richard 

	[[alternative HTML version deleted]]

