[R-sig-Geo] raster::overlay with a function that calls a list of vectors

[Hugo Costa]

Dear list

I wonder if it’s possible to use overlay function of raster package with a
function that calls a list of vectors to perform some calculations based on
two rasters. So far I just saw examples of functions performing some raster
algebra without calling external data.

Hereafter I provide some toy code to illustrate what I’m trying to do, but
can also provide some context on my real problem. Specifically, I need to
classify each pixel as either zero (absence) or one (presence) of housings.
The likelihood of housing presence is related to the percentage of built-up
area covering the pixel (raster ‘r1’ below), and the land cover type
(raster ‘r2’ below). This likelihood is known based on reference data,
which is stored in a list like ‘probs’ below.

library(raster)
# continuous and categorical maps
r1<-r2<-raster()
r1[]<-round(runif(ncell(r1))*100)
r2[]<-1
r2[1:30000]<-2
# probability of housing presence in each stratum
prob1<-1:100/100
prob2<-log(1:100)/max(log(1:100))
# list of probabilities to be used in overlay
probs<-list(prob1,prob2)
# overlay - not working
o<-overlay(r1,r2,fun=function(x,y,...){return(rbinom(n=1, size=1,
prob=probs[[y]][x]))})

the error is
cannot use this formula, probably because it is not vectorized

Alternatively to the toy code above, I thought to process each categorical
class separately and use function calc rather than function overlay (see
below). However, this is extremely slow (if not impossible) for large
rasters, so I though overlay would be better.

# alternative: loop across categorical classes (extremely slow for
large rasters)
r<-list()for(i in 1:2){
  stratum<-r2
  stratum[Which(stratum !=i)]<-NA
  r[[i]]<-calc(r1, fun=function(x,...){return(rbinom(n=1, size=1,
prob=probs[[i]][x]))})
  r[[i]]<-mask(r[[i]],stratum)}

r<-stack(r)
r<-sum(r,na.rm=T)

par(mfrow=c(1,3))
plot(r1)
plot(r2)
plot(r)

	[[alternative HTML version deleted]]