[R-sig-Geo] apply function to lists of rasters

[Robert J. Hijmans]

This mapply problem has now been fixed for the next release of raster (>= 2.1-6)
https://stat.ethz.ch/pipermail/r-devel/2013-February/065847.html
Also thanks to Roman Luštrik for off list suggestions. Robert

On Mon, Feb 11, 2013 at 9:14 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
>
> mapply does indeed not seem to work here, further illustrated below. I do
> not know why. However, in this case I think it is simpler not to use mapply
> (see example at the very end).
>
>> library(raster)
> Loading required package: sp
> raster 2.1-3 (1-February-2013)
>>
>> aa<-matrix(c(2,1,3,4),nrow=2,ncol=2)
>> bb<-matrix(c(5,7,4,6),nrow=2,ncol=2)
>> cc<-matrix(c(12,13,17,16),nrow=2,ncol=2)
>> A<-list(raster(aa),raster(t(aa)))
>> B<-list(raster(bb),raster(t(bb)))
>> C<-list(raster(cc),raster(t(cc)))
>>
>>
>> f1 <- function(a,b,c) 1/(1-exp(-exp(-((a - b)/c))))
>
> # this function does not work, I suppose the "-" function, as shortcut for
> "-1 *" is not implemented for Raster objects
>> f1(A[[1]], B[[1]], C[[1]])
> Error in -((a - b)/c) : invalid argument to unary operator
>
> # You can make it work like this, because "-1 *" is understood
>> f2 <- function(a,b,c) 1/(1-exp(-1*exp(-1*((a - b)/c))))
>> f2(A[[1]], B[[1]], C[[1]])
> class       : RasterLayer
> dimensions  : 2, 2, 4  (nrow, ncol, ncell)
> resolution  : 0.5, 0.5  (x, y)
> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> data source : in memory
> names       : layer
> values      : 1.257289, 1.529642  (min, max)
>
> # now on to mapply
> # simple fucnction
>> f3 <- function(a,b,c) (a - b)/c
>
> # works
>> D <- mapply(f3, A,B,C)
>
> # original function does not work
>> D <- mapply(f2, A,B,C)
> Error in as.character(sys.call(sys.parent())[[1]]) :
>   cannot coerce type 'closure' to vector of type 'character'
> # probably because "exp" is part of the Math group generic and the
> # method for Raster objects gets lost along the way.
> # I need to investigate that. Ideas anyone?
>
>
> # However, I think the good news is that a more rasteresque solution works:
>> f2(stack(A), stack(B), stack(C))
> class       : RasterBrick
> dimensions  : 2, 2, 4, 2  (nrow, ncol, ncell, nlayers)
> resolution  : 0.5, 0.5  (x, y)
> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> data source : in memory
> names       :  layer.1,  layer.2
> min values  : 1.257289, 1.257289
> max values  : 1.529642, 1.529642
>
> # OR
> overlay(stack(A), stack(B), stack(C), fun=f2)
> class       : RasterBrick
> dimensions  : 2, 2, 4, 2  (nrow, ncol, ncell, nlayers)
> resolution  : 0.5, 0.5  (x, y)
> extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
> coord. ref. : NA
> data source : in memory
> names       :  layer.1,  layer.2
> min values  : 1.257289, 1.257289
> max values  : 1.529642, 1.529642
>
>>
>
>
>
> On Mon, Feb 11, 2013 at 5:31 AM, Roman Luštrik <roman.lustrik at gmail.com>
> wrote:
>>
>> There may be a problem.
>>
>> mapply(function(a, b, c) exp(a), A, B, C)
>>
>>
>> alas
>>
>> > exp(A[[1]])
>> > class       : RasterLayer
>> > dimensions  : 2, 2, 4  (nrow, ncol, ncell)
>> > resolution  : 0.5, 0.5  (x, y)
>> > extent      : 0, 1, 0, 1  (xmin, xmax, ymin, ymax)
>> > coord. ref. : NA
>> > data source : in memory
>> > names       : layer
>> > values      : 2.718282, 54.59815  (min, max)
>>
>>
>> Looks like exp() doesn't work from inside `mapply`. Robert?
>>
>> # I changed c to d, because it's hell to debug `c`
>>
>> D <- mapply(function(a, b, d) {
>> >   browser()
>>
>>   1/(1 - exp(-exp(-((a - b)/d))))
>> > }, A, B, C)
>>
>>
>>
>> Browse[1]> class(a)
>> > [1] "RasterLayer"
>> > attr(,"package")
>> > [1] "raster"
>> > Browse[1]> exp(a)
>> > Error in as.character(sys.call(sys.parent())[[1]]) :
>> >   cannot coerce type 'closure' to vector of type 'character'
>>
>>
>>
>>
>> Cheers,
>> Roman
>>
>>
>>
>> On Mon, Feb 11, 2013 at 2:14 PM, Lorenzo Alfieri
>> <alfios17 at hotmail.com>wrote:
>>
>> >
>> > Dear list,
>> > I'm trying to apply a function to lists of rasters, but without success.
>> > I'm not sure if it feasible of if I should use a for loop.
>> > Let's say a, b and c are three lists of argument to my function
>> >  1/(1-exp(-exp(-((a - b)/c))))  (2 rasters each)
>> > I'd like this function to be applied to these arguments, so that I get
>> > as
>> > result a list D of 2 rasters, where for example
>> > D[[1]] = 1/(1-exp(-exp(-((A[[1]] - B[[1]])/C[[1]])))) and the same for
>> > D[[2]]
>> >
>> > I tryed with this but didn't work:
>> >
>> > library(raster)
>> >
>> > aa<-matrix(c(2,1,3,4),nrow=2,ncol=2)
>> > bb<-matrix(c(5,7,4,6),nrow=2,ncol=2)
>> > cc<-matrix(c(12,13,17,16),nrow=2,ncol=2)
>> > A<-list(raster(aa),raster(t(aa)))
>> > B<-list(raster(bb),raster(t(bb)))
>> > C<-list(raster(cc),raster(t(cc)))
>> >
>> > D<-mapply(function(a,b,c) 1/(1-exp(-exp(-((a - b)/c)))), A,B,C)
>> >
>> > Any suggestions?
>> >
>> > Lorenzo
>> >
>> >
>> >
>> >         [[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
>> >
>>
>>
>>
>> --
>> In God we trust, all others bring data.
>>
>>         [[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
>
>