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 |
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 _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
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 <[hidden email]> 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 <[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 > > [[alternative HTML version deleted]] _______________________________________________ R-sig-Geo mailing list [hidden email] https://stat.ethz.ch/mailman/listinfo/r-sig-geo |
Free forum by Nabble | Edit this page |