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 <zenli...@gmail.com> wrote:
From: Brian Oney <zenli...@gmail.com>
Subject: Re: [R-sig-Geo] Bioclimatic variables - wettest quarter
To: "Jacob van Etten" <jacobvanet...@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 <zenli...@gmail.com>
wrote:
From: Brian Oney <zenli...@gmail.com>
Subject: [R-sig-Geo] Bioclimatic variables - wettest
quarter
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]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
