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