raster::overlay with a function that calls a list of vectors

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

raster::overlay with a function that calls a list of vectors

Hugo Costa
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
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[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
Reply | Threaded
Open this post in threaded view
|

Re: 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 <[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