Querying rasters to create result raster

Question

I have two environmental rasters for the same extent etc, and dataframe of model results (which contains all possible combinations of env1 and env2). I am trying to create a third raster filled with model results for each cell.

library(raster)
## Mock env rasters
env.param1 <- seq(1:100)
env.param2 <- seq(101:200)

ext <- extent(1,10,1,10)
env1 <- raster(ext,nrow=10,ncol=10)
values(env1) <- env.param1
env2 <- raster(ext,nrow=10,ncol=10)
values(env2) <- env.param2

## Mock model results
param.set <- expand.grid(env.param1,env.param2)
res.v <- seq(10001:20000)
res <- cbind(param.set, res.v)
res.rast <- raster(ext, nrow=10, ncol=10)
## 
for(i in 1:ncell(env1)){
  res.rast[i] <- res[which(res$Var1==env1[i] & res$Var2==env2[i]),"res.v"]  
}

This appears to work, but it fails when a value in env* does not occur in res which will happen in my actual data set. I've come up with a solution, but it seems slow.

for(i in 1:ncell(env1)){

  res.rast[i] <-res[which(ifelse(abs(res$Var1-env1[i])==min(abs(res$Var1-env1[i])),TRUE,FALSE)
            & ifelse(abs(res$Var2-env2[i])==min(abs(res$Var2-env2[i])),TRUE,FALSE)),"res.v"]  
}

Is there a way to optimise this to A) increase speed as I use higher resolution env* rasters (I have a new dataset that is 225x higher resolution than my old one)? B) How can I scale when I begin to include more parameters in my analysis.

score 0 · Answer 1 · answered Dec 04 '16 at 20:24

This should work:

for(env1ind in unique(param.set$Var1)) {
  for(env2ind in unique(param.set$Var2)) {

    cells           <- which(env.param1 == env1ind &  env.param2 == env2ind)
    res.rast[cells] <- as.numeric(res[which(res$Var1 == env1ind & res$Var2 == env2ind),][3])

  }
}

and should be faster than cycling on pixels, but I cannot test because your "working" solution gives me an error.

This is however not at all efficient nor optimized. And if you have big images, having to load all data in memory beforehand via values can be problematic.

If your "model" is not complex (i.e., computationally expensive), I'd consider putting all the input rasters together via raster::stack , and then using calc to apply the model. That would allow you to work by "image blocks", and also to use parallel processing if needed.

HTH

Querying rasters to create result raster

1 Answers1