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)
raster::overlay with a function that calls a list of vectors
3 messages · Rafael Wüest, Hugo Costa
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 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]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Rafael W?est Karpati http://www.rowueest.net
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 at 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 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]]
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-geo
-- Rafael W?est Karpati http://www.rowueest.net