3

I have a RasterBrick object that represents a raster time series (say different images of the same scene over a thirty year period) and I'm trying to perform a Theil-Sen regression on this time series.

I'm attempting to do this through the zyp.sen() and calc() functions but I've run into a roadblock.

Here's what I've tried running:

 temp <- list()
 for(i in 1:10) {
     temp[[i]] <- runif(36,-1, 1)
     temp[[i]] <- matrix(temp[[i]], 6, 6)
     temp[[i]] <- raster(temp[[i]])}
 temp <- brick(temp)

 time <- 1:nlayers(temp)

 fun=function(x) { if (is.na(x[1])){ NA } else { m = zyp.sen(x~time)}}

 calc(temp,fun)

However, I keep running into this error:

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

And I'm not really sure why, or what I'm doing wrong.

6

The problem is not entirely yours. zyp.sen is not well written: it fails unless you guess correctly concerning the syntax it is hoping for. (Whether that is its only problem is a question I have not pursued, but you might want to.)

One way to verify this conclusion is to test your function--which is a wise idea in any case, even when you are confident in the software you are using.

x <- rnorm(length(time))
fun(x)

This fails. The error is

Error in yy[1:(n - i)] : only 0's may be mixed with negative subscripts

A quick debugging session indicates the problem is in zyp.sen (a part of the zyp package), which attempts to index an array with the value 0-1.

This behavior would be enough to make me look for a different package (or even write my own Theil-Sen regression, since it's easy to do). But if you must use this particular one, here's a workaround. It coddles zyp.sen to match precisely what it seems to be looking for: an explicit data frame argument and explicit naming of all variables in the formula. I have had to modify your function, though, because it was doomed to failure anyway: you were attempting to return the entire Theil-Sen model and store it in the brick. That would produce an error in any event. You need to choose which part(s) of the model you want: the coefficients? Just the slope? Something else? This example returns just the slope:

library(raster)
library(zyp)
fun <- function(y) { 
  if (any(is.na(y))){ 
    NA_real_ 
  } else { 
    X <- data.frame(y=y, time=time)
    model <- zyp.sen(y ~ time, X)
    beta.hat <- coef(model)
    beta.hat["time"]
  }
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.