Hi Jacob,

You are right it does produce a stack. I did not know that calc can 
handle a computation like that for stack, impressive.

So I have tried your suggestions and hate to keep bothering but it still 
is not completely working, but we are making progress:

func <- function(X) {x <- movingFun(c(X), n=3, sum, "around", 
circular=TRUE);return(x)}
wetness_quarters <- calc(mon_stac, func)
# It worked!
summary(wetness_quarters)
#Error in cbind(sumobj@values, as.matrix(sm)) :
#  number of rows of matrices must match (see arg 2)

str(wetness_quarters) # looks ok

wettest_quarter <- calc(wetness_quarters, function(x) which(x == max(x)))
#Error in setValues(out, x) : values must be numeric, integer or logical.

wettest_quarter <- calc(wetness_quarters, which.max)
#Error in .local(x, values) : values must be numeric, integer or logical.

Could it be that I have way too many places that have only NA values? 
i.e. NA values over the ocean?
?which.max tells me that NA's are discarded.

I know it is UGLY. I just didn't have a better idea, and I would have 
never dreamed up your proposed function. I am a little green, but I 
really appreciate the help.

Cheers,
Brian






On 11/8/2010 4:58 PM, Jacob van Etten wrote:
> Hi Brian,
>
> Something strange is going on. Normally, this should produce a stack.
>
> This works, for instance:
>
> fn <- system.file("external/test.grd", package="raster")
> s <- stack(fn, fn)
> t <- calc(s, function(x) x*10)
> plot(t)
>
> I think something strange is going on with your stack. Try some basic 
> operations on it like
>
> plot(mon_stac * 10)
>
> to see what´s going on. Also, check if you have the latest version of 
> raster, etc.
>
> Regarding your question what to do with the output you get with my 
> suggestion: well, it can be used to index the temperature raster to 
> retrieve the values you want. You can append the index raster to the 
> temperature raster
>
> s <- stack(mon_temp, wettest_quarter)
>
> and then use calc() with a function like this (assuming you took the 
> middle month, type="around", not type="from"):
>
> func <- function(x)
> {
>   from <- c(12,1:11)[x[13]]
>   to <- c(2:12,1)[x[13]]
>   return(mean(x[from:to]))
> }
>
> overlay() should also work, but you get these ugly functions with 12 
> arguments!
>
> Best,
>
> Jacob.
>
> --- On *Mon, 8/11/10, Brian Oney /<zenlines@gmail.com>/* wrote:
>
>
>     From: Brian Oney <zenlines@gmail.com>
>     Subject: Re: [R-sig-Geo] Bioclimatic variables - wettest quarter
>     To: "Jacob van Etten" <jacobvanetten@yahoo.com>
>     Cc: r-sig-geo@stat.math.ethz.ch
>     Date: Monday, 8 November, 2010, 14:48
>
>     Hi Jacob,
>     Thanks for the suggestion. So...
>
>     func <- function(X) {x <- movingFun(c(X), n=3, sum, 'from',
>     circular=TRUE);return(x)}
>     mon_stac <- stack(paste(pc_clim,"",y,"_prec",1:12,".asc", sep=""))
>     # then
>     wetness_quarters <- calc(mon_stac, func)
>     Error in .local(x, fun, ...) : function 'fun' returns more than
>     one value
>
>     # This make sense that this would do this..
>     # Now this to get the first month of the max quarter by returning
>     first month number:
>     func <- function(X) {x <- which.max(movingFun(c(X), n=3, sum,
>     'from', circular=TRUE));return(x)}
>
>     # Then
>     wetness_quarters <- calc(mon_stac, func)
>     Error in v[, tr$row[i]:(tr$row[i] + tr$nrows[i] - 1)] <-
>     matrix(sv, nrow = ncol(outraster)) :
>       incorrect number of subscripts on matrix
>     # This error being returned after a long time busying a quadcore
>     computer...
>
>     Even if I can get your suggestion to work, I am still faced with
>     the problem of retrieving the sum of values of the quarter marked
>     my the month number,
>     which should look something like this:
>
>     quart_value <- function(a,b,d,e,f,g,i,j,k,l,m,n,o){
>                         ifelse(o==1, X <- mean(c(a,b,d)),
>                         ifelse(o==2, X <- mean(c(b,d,e)),
>                         ifelse(o==3, X <- mean(c(d,e,f)),
>                         ifelse(o==4, X <- mean(c(e,f,g)),
>                         ifelse(o==5, X <- mean(c(f,g,i)),
>                         ifelse(o==6, X <- mean(c(g,i,j)),
>                         ifelse(o==7, X <- mean(c(i,j,k)),
>                         ifelse(o==8, X <- mean(c(j,k,l)),
>                         ifelse(o==9, X <- mean(c(k,l,m)),
>                         ifelse(o==10, X <- mean(c(l,m,n)),
>                         ifelse(o==11, X <- mean(c(m,n,a)),
>                         ifelse(o==12, X <- mean(c(n,a,b)),NA))))))))))))
>                         return(X)}
>
>     #For BIO8
>     mon_stac.wq <- stack(paste(pc_clim,"",y,"_tmean",1:12,".asc",
>     sep=""), wetness_quarters)
>     tmean.wq <- overlay(mon_stac.wq, fun=quart_value)
>
>     With "o" being the product of the "func"or "maxmonth" function
>     I think that I NEED to get the "overlay" function to work, but am
>     of course open to suggestions.
>     I appreciate the help.
>
>     Regards,
>     Brian
>
>
>     On 11/8/2010 10:59 AM, Jacob van Etten wrote:
>>     I don´t know what the problem is with overlay(), but you could
>>     use an alternative approach, using calc() instead.
>>
>>     Something like this (not tried)
>>
>>     library(raster)
>>     library(dismo)
>>     func <- function(x, ...) movingFun(x, n=3, sum, circular=TRUE)
>>     wetness_quarters <- calc(mon_stac, func)
>>     wettest_quarter <- calc(wetness_quarters, function(x) which(x ==
>>     max(x))
>>
>>     The output input should indicate the month in the centre of the
>>     wettest quarter.
>>
>>     Jacob.
>>
>>     --- On *Mon, 8/11/10, Brian Oney /<zenlines@gmail.com>
>>     </mc/compose?to=zenlines@gmail.com>/* wrote:
>>
>>
>>         From: Brian Oney <zenlines@gmail.com>
>>         </mc/compose?to=zenlines@gmail.com>
>>         Subject: [R-sig-Geo] Bioclimatic variables - wettest quarter
>>         To: r-sig-geo@stat.math.ethz.ch
>>         </mc/compose?to=r-sig-geo@stat.math.ethz.ch>
>>         Date: Monday, 8 November, 2010, 1:40
>>
>>         Hello All,
>>
>>         I am attempting to reproduce the bioclimatic (Worldclim,
>>         Hijmans et al
>>         2005) variables with the raster package and I am stuck on the
>>         group
>>         "temp of driest quarter" etc. variables.
>>         I am attempting it with the function "overlay" ("raster"
>>         package).
>>
>>         Here is how it looks:
>>         # To know the first month of the wettest quarter
>>         maxmonth<- function(a,b,d,e,f,g,i,j,k,l,m,n){
>>                             o<- sum(c(a,b,d))
>>                             p<- sum(c(b,d,e))
>>                             r<- sum(c(d,e,f))
>>                             u<- sum(c(e,f,g))
>>                             v<- sum(c(f,g,i))
>>                             w<- sum(c(g,i,j))
>>                             x<- sum(c(i,j,k))
>>                             y<- sum(c(j,k,l))
>>                             z<- sum(c(k,l,m))
>>                             aa<- sum(c(l,m,n))
>>                             ab<- sum(c(m,n,a))
>>                             ac<- sum(c(n,a,b))
>>                             ad<- which.max(c(o,p,r,u,v,w,x,y,z,aa,ab,ac))
>>                             return(ad)}
>>         > 
>>         maxmonth(a=13,b=16,d=41,e=61,f=41,g=16,i=15,j=14,k=13,l=11,m=31,n=11)
>>         # works: 2nd element (month) is the beginning of the wettest
>>         quarter
>>         [1] 2
>>         mon_stac<- stack(paste(pc_clim,"",y,"_prec",1:12,".asc", sep=""))
>>         # Now attempt it
>>         avg<- overlay(mon_stac, fun=maxmonth)
>>         Error in .overlayList(rasters, fun = fun, filename =
>>         filename, ...) :
>>           cannot use this formula; lenghts do not match
>>
>>         # How about manually?
>>         a<- raster(paste(pc_clim,"",y,"_prec",1,".asc", sep=""))
>>         b<- raster(paste(pc_clim,"",y,"_prec",2,".asc", sep=""))
>>         d<- raster(paste(pc_clim,"",y,"_prec",3,".asc", sep=""))
>>         e<- raster(paste(pc_clim,"",y,"_prec",4,".asc", sep=""))
>>         f<- raster(paste(pc_clim,"",y,"_prec",5,".asc", sep=""))
>>         g<- raster(paste(pc_clim,"",y,"_prec",6,".asc", sep=""))
>>         i<- raster(paste(pc_clim,"",y,"_prec",7,".asc", sep=""))
>>         j<- raster(paste(pc_clim,"",y,"_prec",8,".asc", sep=""))
>>         k<- raster(paste(pc_clim,"",y,"_prec",9,".asc", sep=""))
>>         l<- raster(paste(pc_clim,"",y,"_prec",10,".asc", sep=""))
>>         m<- raster(paste(pc_clim,"",y,"_prec",11,".asc", sep=""))
>>         n<- raster(paste(pc_clim,"",y,"_prec",12,".asc", sep=""))
>>
>>         avg<- overlay(a,b,d,e,f,g,i,j,k,l,m,n, fun=maxmonth)
>>         Error in .overlayList(rasters, fun = fun, filename =
>>         filename, datatype
>>         = datatype,  :
>>           cannot use this formula; lenghts do not match
>>
>>         Also doesn't work...
>>         >  length(mon_stac[1])
>>         [1] 12
>>         I have 12 rasters being fed into overlay, but it won't buy it.
>>         I think that my first step is to know the first month of the
>>         wettest
>>         month. Once I can do that I can do it for temperature as well.
>>
>>         Thanks and Regards,
>>         Brian
>>
>>         _______________________________________________
>>         R-sig-Geo mailing list
>>         R-sig-Geo@stat.math.ethz.ch
>>         https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>

	[[alternative HTML version deleted]]

