[R-sig-Geo] raster stack/brick cor.test p-values

[Pascal Oettli]

Please reply to the list, not privately.


> Thanks for the tip. It's a bit strange because the same code, with a minor
> modification, ran well until all of sudden it decided to spit out these
> error. Here it is:
>
> brick.cor.pval <- function(xy) { # Pearson product moment correlation
> p-values
>   xy <- na.omit(matrix(xy, ncol=2))
>   if (ncol(xy) < 2) {
>     NA} else if
>   (is.na(xy[1])) {
>     NA} else
>   { cor.test(xy[, 1], xy[, 2],method="pearson")$p.value
>   }
> }
>
> I just added the "(is.na(xy[1]))" portion. Any ideas as to why this is?
>
>
>
> On Thu, May 29, 2014 at 10:56 AM, Pascal Oettli <kridox at 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 at 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 at nateko.lu.se
>> >
>> >         [[alternative HTML version deleted]]
>> >
>> >
>> > _______________________________________________
>> > R-sig-Geo mailing list
>> > R-sig-Geo at r-project.org
>> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>> >
>>
>>
>>
>> --
>> Pascal Oettli
>> Project Scientist
>> JAMSTEC
>> Yokohama, Japan
>
>



-- 
Pascal Oettli
Project Scientist
JAMSTEC
Yokohama, Japan