Thanks for pointing that out, Pascal.

By the way, if anyone wants sample data to try out these codes, I have them
here:

Brick 1: https://www.dropbox.com/s/bzgz2w7ao8kld27/pet.zip
Brick 2: https://www.dropbox.com/s/jyuc0u0pt3bozn4/rain.zip

I tried to convert the data into a numeric vector:

cor.pval <- function(xy) { # Pearson product moment correlation
>   xy <- na.omit(as.matrix(xy, ncol=2))
>   x <- xy[,1:11]
>   y <- xy[,12:22]
>     cor.test(x[1,], y[1,])$p.value
> }
> calc(s, cor.pval)


But I got the following error:

Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
>   cannot use this function


So, I then tried to check the number of rows as well as suggested by
another user and modified the first code:

funcal <- function(xy) {
>     xy <- na.omit(matrix(xy, ncol=2))
>     if (ncol(xy) < 2 | nrow(xy) < 1 ) {
>         NA
>     } else {
>         cor.test(xy[, 1], xy[, 2], method="pearson")$p.value
>     }
> }


I varied the condition of "nrow(xy) < 1" between 2 and 6 to see what I
output I would get. At value 2, I get the aforementioned "not enough finite
observations" response. From values 3 and up, I get a blank raster as the
output, as shown below:

class       : RasterLayer
> dimensions  : 273, 146, 39858  (nrow, ncol, ncell)
> resolution  : 0.05, 0.05  (x, y)
> extent      : 43.2, 50.5, -25.6, -11.95  (xmin, xmax, ymin, ymax)
> coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
> +towgs84=0,0,0
> data source : in memory
> names       : layer
> values      : NA, NA  (min, max)


Does anyone have any ideas on how to fix this?

My sessionInfo() is below:

R version 3.0.1 (2013-05-16)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
> locale:
> [1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United
> States.1252
> [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
>
> [5] LC_TIME=English_United States.1252
> attached base packages:
> [1] grid      stats     graphics  grDevices utils     datasets  methods
> base
> other attached packages:
>  [1] car_2.0-18          nnet_7.3-7          MASS_7.3-29
> rgdal_0.8-11
>  [5] plyr_1.8            rasterVis_0.21      hexbin_1.26.2
> latticeExtra_0.6-26
>  [9] RColorBrewer_1.0-5  lattice_0.20-23     raster_2.1-49
> maptools_0.8-27
> [13] sp_1.0-13
> loaded via a namespace (and not attached):
> [1] foreign_0.8-55 tools_3.0.1    zoo_1.7-10




On Thu, May 29, 2014 at 10:56 AM, Pascal Oettli <kridox@ymail.com> wrote:

> Hello Hakim,
>
> The function "cor" can take matrix as input. Not the function "cor.test".
>
> Hope this helps,
> Pascal
>
>
> On Thu, May 29, 2014 at 5:45 PM, Hakim Abdi <hakim.abdi@nateko.lu.se>
> wrote:
> > Dear list:
> >
> > I am trying to get pixel-wise correlations and significance between two
> > raster bricks; one representing monthly NDVI means and another
> representing
> > monthly rainfall means.
> >
> > I found the following two codes (both from Robert Hijmans) from a
> previous
> > discussions on this list:
> >
> http://r-sig-geo.2731867.n2.nabble.com/questions-on-RasterStack-Brick-td5553580.html
> > and
> > from StackOverflow:
> >
> http://stackoverflow.com/questions/16698048/r-calculating-pearson-correlation-coefficient-in-each-cell-along-time-line
> >
> >
> > # code 1
> >> funcal <- function(xy) {
> >>     xy <- na.omit(matrix(xy, ncol=2))
> >>     if (ncol(xy) < 2) {
> >>         NA
> >>     } else {
> >>         cor(xy[, 1], xy[, 2])
> >>     }
> >> }
> >> s <- stack(sa, sb)
> >> calc(s, funcal)
> >> # code 2
> >> stackcor <- function(s1, s2, method='spearman') {
> >>         mycor <- function(v) {
> >>                 x <- v[1:split]
> >>                 y <- v[(split+1):(2*split)]
> >>                 cor(x, y, method=method)
> >>         }
> >>         s <- stack(s1, s2)
> >>         split <- nlayers(s)/2
> >>         calc(s, fun=mycor )
> >> }
> >
> >
> > Both codes work as expected with the "cor" function, producing
> correlation
> > grids. However, I tried substituting the cor with cor.test in order to
> > extract the p-values:
> >
> > # using the first code
> >> funcal <- function(xy) {
> >>     xy <- na.omit(matrix(xy, ncol=2))
> >>     if (ncol(xy) < 2) {
> >>         NA
> >>     } else {
> >>         cor.test(xy[, 1], xy[, 2],method="pearson")$p.value
> >>     }
> >> }
> >> s <- stack(sa, sb)
> >> calc(s, funcal)
> >
> >
> > I am met with the following error (with trackback in RStudio):
> >
> > Error in cor.test.default(xy[, 1], xy[, 2], method = "pearson") :
> >>   not enough finite observations
> >> 8 stop("not enough finite observations")
> >> 7 cor.test.default(xy[, 1], xy[, 2], method = "pearson")
> >> 6 cor.test(xy[, 1], xy[, 2], method = "pearson")
> >> 5 FUN(newX[, i], ...)
> >> 4 apply(v, 1, fun)
> >> 3 .local(x, fun, ...)
> >> 2 calc(meistack, brick.cor.pval)
> >> 1 calc(meistack, brick.cor.pval)
> >
> >
> > In the previous r-sig-geo, a user had asked about this error but it was
> > left unanswered. I was wondering if anyone can help with this?
> >
> > Thanks!
> >
> > Hakim
> > ____________________________________________________
> >
> > Hakim Abdi |  PhD Candidate
> > Geobiosphere Science Center
> >
> > Department of Physical Geography and Ecosystem Science
> > Lund University, Sölvegatan 12, SE-223 62 Lund, Sweden
> >
> > Email: hakim.abdi@nateko.lu.se
> >
> >         [[alternative HTML version deleted]]
> >
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo@r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>
>
> --
> Pascal Oettli
> Project Scientist
> JAMSTEC
> Yokohama, Japan
>

	[[alternative HTML version deleted]]

