Raster correlation and p-values from cor.test

Question

I am trying to get pixel-wise correlations and significance (p-value) between two raster bricks using cor and cor.test. My data are here:

• Brick 1
• Brick 2

They're both fairly small, less than 2MB altogether.

I found the following two codes (both from Robert Hijmans) from a previous discussions on StackOverflow and r-sig-geo:

#load data
require(raster)
sa <- brick("rain.grd")
sb <- brick("pet.grd")

# 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 a previous r-sig-geo discussion, a user had asked about this error but it was left unanswered. So I asked again and received one response to my inquiry where a user pointed out that cor is able to input matrices and cor.test cannot, but even after converting 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]
  #  if (ncol(xy) < 2) {
   # NA
  #} else {
    cor.test(x[1,], y[1,])$p.value
  #}
}
calc(s, cor.pval)

I am faced with the following error:

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

I was wondering if anyone can help with 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    

Thanks!

Answer 1

The difference here is that they behave differently with empty vectors.

cor(numeric(0), numeric(0))
# -> NA
cor.test(numeric(0), numeric(0))
#->   Error in cor.test.default(numeric(0), numeric(0)) : 
#       not enough finite observations

It seems that your na.omit can remove all the records from certain combinations. Right now you only check the number of columns, when you should also check to see if there are any rows.

This

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
    }
}