Hi Brian,
My guess is that it goes wrong when this happens:
which.max(c(NA,NA))
Try something like this instead:
func <- function (x) {
if(all(is.na(x))) {result <- 0
} else {result <- which.max(x)}
return(result)
}
Best,
Jacob.
--- On Tue, 9/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: Tuesday, 9 November, 2010, 2:11
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(sum...@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 <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
