[R-sig-Geo] Error running Mann-Kendall trend test on a raster stack

[Thiago V. dos Santos]

Hi Tim and Loïc,

Thanks for your suggestions.

I was able to run the analysis using the following function (which also accounts for NA values in the raster object):

----------------
library(raster)
library(rkt)

# Create the date sequence
idx <- seq.Date(as.Date("2011-01-01"), as.Date("2099-01-01"), by='year')

# Create raster stack and apply the date
r <- raster(ncol=23, nrow=19)
s <- stack(lapply(1:length(idx), function(x) setValues(r, runif(ncell(r), 0, 1200))))
s <- setZ(s, idx)
s

year <- as.numeric(substr(getZ(s), 1, 4))
# Now I define a function
tsfun <- function(x) {
if(all(is.na(x))){
c(NA, NA)
} else {
r.kenn <- rkt(year, x)
a <- r.kenn$B
b <- r.kenn$sl
return(cbind(a, b))
}

# I then apply the function - takes a while to run!
raster.kenn <- calc(s, fun=tsfun)
raster.kenn
----------------
Greetings,
 -- Thiago V. dos Santos

PhD student
Land and Atmospheric Science
University of Minnesota



On Saturday, February 27, 2016 5:19 AM, Loïc Dutrieux <loic.dutrieux at wur.nl> wrote:
Hi Thiago,

A few calc tips inline that you can use together with Tim's suggestion.

On 02/27/2016 07:47 AM, Tim Appelhans wrote:
> Dear Thiago,
> have a look at the gimms package, available on CRAN. The function of
> interest is significantTau().
>
> https://github.com/environmentalinformatics-marburg/gimms/blob/master/R/significantTau.R
>
>
> Hope this helps,
> Tim
>
> On 27.02.2016 09:49, Thiago V. dos Santos wrote:
>> Dear colleagues,
>>
>> I have a raster stack with 89 layers, each layer representing yearly
>> precipitation.
>>
>> I am trying to use the function rkt (from package "rkt") to detect a
>> possible trend in my precipitation time series.
>>
>> Its usage is pretty simple, and this is how it runs on a data frame:
>>
>> --------------------------------
>> library(rkt)
>>
>> # generate some random data
>> myDF <- data.frame(Date = seq.Date(as.Date("2011-01-01"),
>> as.Date("2099-01-01"), by='year'),
>> value = runif(89,0,1200))
>>
>> # extract year from date, as numeric
>> myDF$year <- as.numeric(format(myDF$Date, "%Y"))
>>
>> # run mann-kendall test
>> kenn <- rkt(myDF$year, myDF$value)
>> kenn
>>
>> --------------------------------
>>
>> Now, this is my attempt to apply this test on a raster stack, using a
>> calc function:
>>
>> --------------------------------
>> library(raster)
>>
>> # Create the date sequence
>> idx <- myDF$Date
>>
>> # Create raster stack and apply the date
>> r <- raster(ncol=50, nrow=50)
>> s <- stack(lapply(1:length(idx), function(x) setValues(r,
>> runif(ncell(r), 0, 1200))))
>> s <- setZ(s, idx)
>> s
>>
>> # Now I define a function
>> tsfun <- function(x) {
>>      year <- as.numeric(substr(getZ(x), 1, 4))
>>      r.kenn <- rkt(year, x)
>>      return(r.kenn)
>> }

You can't call getZ() on x. Within the calc function environment x is a 
numeric vector not a RasterBrick.
If you do it outside of the function and refer to it within the function 
it should work.
dates <- getZ(s) # then use the dates variable in the function.

You can only store numerics in raster cells, and r.kenn is a full 
object. You need to choose which numerics contained in that object you 
want the function to return and subset the object accordingly.

Functions passed to calc should work on numeric vectors and ideally they 
should also be vectorized. It's not always easy to vectorize a function 
though; so in such case you can set the forceapply= argument to TRUE. I 
believe that may cost a bit of speed though.

Cheers,
Loïc

>>
>> # I apply the function
>> raster.kenn <- calc(s, fun=tsfun)
>>
>> --------------------------------
>>
>>
>> And the error message I get is:
>>
>> Error in .calcTest(x[1:5], fun, na.rm, forcefun, forceapply) :
>> cannot use this function
Always loved that very informative error message from calc :)

>>
>> What am I doing wrong here?
>> Greetings,
>>   -- Thiago V. dos Santos
>>
>> PhD student
>> Land and Atmospheric Science
>> University of Minnesota
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo