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

[Rafael Wüest]

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 using gmail.com> 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
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo



--
Rafael Wüest Karpati
http://www.rowueest.net