Yes, the original example is also not correct too. But in this case it would be
CRPS <- overlay(WB_LT, EUE, fun=function(x,y) crpsDecomposition(x,y)$CRPS)
WB_LT[1] <- NA
EUE[1] <- NA
CRPS2 <- overlay(WB_LT, EUE, fun=function(x,y) crpsDecomposition(x,y)$CRPS)
par(mfrow = c(2,1))
plot(CRPS);plot(CRPS2)
which produce:
Error in function (x, fun, filename = "", recycle = TRUE, ...) :
cannot use this formula, probably because it is not vectorized
(i.e., crpsDecomposition(x,y)$CRPS instead of crps(x,y)$crps)
Could it be that this formula cannot be used with overlay?
I'm trying all combinations but doesn't seem to work...
Lorenzo
Date: Thu, 5 Sep 2013 09:43:48 -0700
Subject: Re: [R-sig-Geo] calculate CRPS on rasters
From: vitale232@gmail.com
To: alfios17@hotmail.com
CC: r-sig-geo@r-project.org
The original example you gave also produces a raster of 1 value throughout.
If you use the suggestion Robert made to use overlay and achieve unique values for each cell, it should still run with NA values.
CRPS <- overlay(WB_LT, EUE, fun=function(x,y) crps(x,y)$crps)WB_LT[1] <- NAEUE[1] <- NACRPS2 <- overlay(WB_LT, EUE, fun=function(x,y) crps(x,y)$crps)
par(mfrow = c(2,1))plot(CRPS);plot(CRPS2)
On Thu, Sep 5, 2013 at 9:33 AM, Lorenzo Alfieri wrote:
That's true, but I get as result one value which is the same for all the raster points.
I thought I could apply this function to each point by using the [] but it doesn't work...
for example, if after
library(verification)
library(raster)
set.seed(1)
r <- raster(ncol=10, nrow=10)
pred1 <- init(r, fun=runif)
pred2 <- init(r, fun=runif)
pred3 <- init(r, fun=runif)
EUE <- stack(pred1, pred2, pred3)
WB_LT <- init(r, fun=runif)
EUE[1]<-NA
WB_LT[1]<-NA
# you run:
CRPS2 <- raster(WB_LT)
CRPS2[]<-crpsDecomposition(na.omit(WB_LT[]),na.omit(EUE[]))$CRPS
CRPS2[2]
# you get a different result from
crpsDecomposition(na.omit(WB_LT[2]),na.omit(EUE[2]))$CRPS
#The correct result could be obtained with
for (i in 1:100){
CRPS2[i]<-crpsDecomposition(na.omit(WB_LT[i]),na.omit(EUE[i]))$CRPS
}
#but I'd like to avoid the for loop as not efficient
Also, the na.omit doesn't overcome the problem in this case, so it crashes for points containing NA. In the example above it only runs if you start from i=2 (as i=1 has NA)
for (i in 2:100){
CRPS2[i]<-crpsDecomposition(na.omit(WB_LT[i]),na.omit(EUE[i]))$CRPS
}
Any suggestion?
Lorenzo
Date: Thu, 5 Sep 2013 08:50:15 -0700
Subject: Re: [R-sig-Geo] calculate CRPS on rasters
From: vitale232@gmail.com
To: alfios17@hotmail.com
CC: r-sig-geo@r-project.org
This seems to run:CRPS2[]<-crpsDecomposition(na.omit(WB_LT[]),na.omit(EUE[]))$CRPS
On Thu, Sep 5, 2013 at 6:43 AM, Lorenzo Alfieri wrote:
Dear list,
following this thread, I'm having another problem.
I'm now using the command crpsDecomposition, of the verification package.
It's quite similar to the crps.
library(verification)
library(raster)
set.seed(1)
r <- raster(ncol=10, nrow=10)
pred1 <- init(r, fun=runif)
pred2 <- init(r, fun=runif)
pred3 <- init(r, fun=runif)
EUE <- stack(pred1, pred2, pred3)
WB_LT <- init(r, fun=runif)
CRPS <- raster(WB_LT)
CRPS[]<-crpsDecomposition(WB_LT[],EUE[])$CRPS
Up to here it works fine.
Problems start when there are NA values in the input.
For example:
EUE[1]<-NA
WB_LT[1]<-NA
CRPS[]<-crpsDecomposition(WB_LT[],EUE[])$CRPS
In this case I get as output:
Error in prev[index, 1] - obs[index] :
non-numeric argument to binary operator
In fact if I run:
crpsDecomposition(WB_LT[1],EUE[1])$CRPS
I get as output:
Error in which(obs < prev[, 1]) : subscript out of bounds
One option is replacing all NAs with zeros
EUE[is.na(EUE)]<-0
before running
CRPS[]<-crpsDecomposition(WB_LT[],EUE[])$CRPS
but it's not very correct and in addition it takes ages, as the EUE is a big variable in my case.
Anyone can help me address this issue?
Regards
Lorenzo
> From: alfios17@hotmail.com
> To: r.hijmans@gmail.com
> CC: r-sig-geo@r-project.org
> Subject: RE: [R-sig-Geo] calculate CRPS on rasters
> Date: Fri, 30 Nov 2012 09:49:00 +0100
>
>
> Robert,
> you're right, it works.
> I must have done something wrong yesterday
> I'm trying to optimize the computation time
> I've tested the two alternatives on a predictor of 51 layers of 680x810 pixels
>
> CRPS[] <- crps(obs[],pred[])$crps took 0.98 second
> CRPS <- overlay(obs, pred, fun=function(x,y) crps(x,y)$crps) took 1.45 seconds
>
> anyway, thank you all!
>
> Lorenzo
>
>
> ________________________________
> > Date: Thu, 29 Nov 2012 10:34:36 -0800
> > Subject: Re: [R-sig-Geo] calculate CRPS on rasters
> > From: r.hijmans@gmail.com
> > To: alfios17@hotmail.com
> > CC: r-sig-geo@r-project.org
> >
> > Lorenzo,
> >
> > The below (your example) works for me, using the overlay as I suggested.
> > Perhaps you did not load the verification package?
> >
> > Robert
> >
> > library(verification) #includes function crps
> > library(raster)
> > set.seed(1)
> > r <- raster(ncol=10, nrow=10)
> > pred1 <- init(r, fun=runif)
> > pred2 <- init(r, fun=runif)
> > pred3 <- init(r, fun=runif)
> > pred <- stack(pred1, pred2, pred3)
> > obs <- init(r, fun=runif)
> >
> > CRPS <- overlay(obs, pred, fun=function(x,y) crps(x,y)$crps)
> >
> >
> > > CRPS
> > class : RasterLayer
> > dimensions : 10, 10, 100 (nrow, ncol, ncell)
> > resolution : 36, 18 (x, y)
> > extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
> > coord. ref. : +proj=longlat +datum=WGS84
> > data source : in memory
> > names : layer
> > values : 0.003277372, 0.7726679 (min, max)
> >
> >
> > # More fancy, to get all the crps output:
> >
> > x <- overlay(obs, pred, fun=function(x,y) as.matrix(data.frame(crps(x,y))))
> >
> >
> >
> > On Thu, Nov 29, 2012 at 1:15 AM, Lorenzo Alfieri
> > > wrote:
> >
> > Hi Robert,
> > I was hoping to use the overlay function too, but it seems it doesn't
> > work in this case. See the result:
> >
> > CRPS <- overlay(obs, pred, fun=function(x,y) crps(x,y)$crps)
> >
> > Error in function (x, fun, filename = "", recycle = TRUE, ...) :
> > cannot use this formula, probably because it is not vectorized
> >
> > Lorenzo
> >
> >
> >
> > Date: Wed, 28 Nov 2012 11:25:36 -0800
> > Subject: Re: [R-sig-Geo] calculate CRPS on rasters
> > From: r.hijmans@gmail.com
> > To: alfios17@hotmail.com
> > CC: etiennebr@gmail.com;
> > r-sig-geo@r-project.org
> >
> > Lorenzo,
> > Given that solution, I think you can also express this (in a
> > memory-safe fashion) like this :
> > CRPS <- overlay(obs, pred, fun=function(x,y) crps(x, y)$crps)
> >
> > Robert
> >
> > On Tue, Nov 27, 2012 at 3:11 AM, Lorenzo Alfieri
> > > wrote:
> >
> >
> >
> >
> > Etienne,
> >
> > thank you for the tip
> >
> > Now it runs trough, by using
> >
> >
> >
> > CRPS <- raster(ncol=10, nrow=10)
> >
> > CRPS[] <- crps(obs[],pred[])$crps
> >
> >
> >
> > Lorenzo
> >
> >
> >
> >
> >
> > Date: Mon, 26 Nov 2012 21:37:13 -0500
> >
> > Subject: Re: [R-sig-Geo] calculate CRPS on rasters
> >
> > From: etiennebr@gmail.com
> >
> > To: alfios17@hotmail.com
> >
> > CC: r-sig-geo@r-project.org
> >
> >
> >
> > Lorenzo,
> >
> >
> >
> > I don't know about your specific function, but you can access matrix
> > using the []'s. So maybe something like :
> >
> > CRPS <- raster(ncol=10, nrow=10)
> >
> > CRPS[] <- crps(obs[],pred[])
> >
> >
> >
> >
> >
> > Etienne
> >
> >
> >
> > 2012/11/26 Lorenzo Alfieri
> > >
> >
> >
> >
> > CRPS <- raster(ncol=10, nrow=10)
> >
> >
> >
> > for (i in 1:length(obs)){
> >
> >
> >
> > CRPS[i] <- crps(obs[i],pred[i])$CRPS
> >
> >
> >
> > [[alternative HTML version deleted]]
> >
> >
> >
> > _______________________________________________
> >
> > R-sig-Geo mailing list
> >
> > R-sig-Geo@r-project.org
> >
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> >
> >
> >
>
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
--
Andrew P. VitaleMasters StudentDepartment of GeographyUniversity of Nevada, Reno(412) 915-3632
vitale232@gmail.com
--
Andrew P. VitaleMasters StudentDepartment of GeographyUniversity of Nevada, Reno(412) 915-3632
vitale232@gmail.com
[[alternative HTML version deleted]]