Hi Hugo
Vectorize() sometimes solves this issue.
foo <- function(x,y,...) {return(rbinom(n=1, size=1, prob=probs[[y]][x]))}
fooV <- Vectorize(foo)
Using your toy example, check carefully if
o <- overlay(r1,r2,fun = fooV)
gives you what you are looking for.
HTH
Rafael
On 9 Aug 2018, at 11:51, Hugo Costa <hugoagcosta at gmail.com> wrote:
Dear list
I wonder if it?s possible to use overlay function of raster package with
function that calls a list of vectors to perform some calculations based
two rasters. So far I just saw examples of functions performing some
algebra without calling external data.
Hereafter I provide some toy code to illustrate what I?m trying to do,
can also provide some context on my real problem. Specifically, I need to
classify each pixel as either zero (absence) or one (presence) of
The likelihood of housing presence is related to the percentage of
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
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]]