Hi Rafael.

Your suggestion works perfectly for the toy data.

For my real data, it's not that simple. I'm processing data for the entire

Europe, and thus I split data in tiles. Apparently, Vectorize does not work

2 strata, but my real data has 7). When tiles include all my strata,

this message). Anyway, your message was important to spot this behaviour, I

strata, while trying to speed up the processing.

> 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 <

[hidden email]> wrote:

> >

> > 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]]

> >

> > _______________________________________________

> > R-sig-Geo mailing list

> >

[hidden email]
> >

https://stat.ethz.ch/mailman/listinfo/r-sig-geo>

>

>

> --

> Rafael Wüest Karpati

>

http://www.rowueest.net>

>