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:

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!

score 1 · Accepted Answer · answered May 29 '14 at 20:28

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

answered May 29 '14 at 20:28

MrFlick

195,160
17
277
295

Thanks MrFlick, I tried the modified code above, and go the following error, it's the same as before: `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(s, funcal) 1 calc(s, funcal)` – lamochila May 30 '14 at 07:02
@H.Abdi Looks like `cor.test` may need at least three observations. Try changing it to `nrow(xy)<3`. You might need to keep increasing that till the error goes away. I don't know what the minimum number of finite observations is. – MrFlick May 30 '14 at 07:10
Thanks, I'll try different values and report back. – lamochila May 30 '14 at 07:13
I've tried values 2 to 6. At values 2, I get the "not enough finite observations" response. From values 3 and up, I get a black 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)` Not sure what else to try? – lamochila May 30 '14 at 07:50
@H.Abdi Isn't that the same object that you get when you run the above code with `cor`? – MrFlick May 30 '14 at 14:27
Thanks @MrFlick, your suggestion was very useful in helping figure out the solution. – lamochila Oct 22 '19 at 18:03

score 1 · Answer 2 · answered Aug 17 '21 at 14:14

1

Use the function corLocal from the raster package.

corLocal(stack1, stack2, test = T)

the returned results are rasterBrick with two layers, the first one 'r' value and the second one p.value

answered Aug 17 '21 at 14:14

Feng Tian

41
1
3

Raster correlation and p-values from cor.test

2 Answers2