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

[Hugo Costa]

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
when a tile includes just some of the strata (raster r2 in the toy data has
2 strata, but my real data has 7). When tiles include all my strata,
Vectorize seems to work, although slowly (it is processing while writing
this message). Anyway, your message was important to spot this behaviour, I
might need to include some workaround to process the tiles with less
strata, while trying to speed up the processing.

Thank you very much
Hugo

Rafael Wüest <rafael.wueest using gmail.com> escreveu no dia quinta, 9/08/2018
à(s) 14:33:

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

	[[alternative HTML version deleted]]