I try to make a prediction in raster using linear regression based on past rasters but I can not capture the lm() model in order to feed it to predict().

# calculate regression in a raster stack (time series)
# and predict.

# make the list of rasters
rasters <- list.files(pattern='*.tif$')

# make the raster stack
s <- stack(rasters)

# crop the stack to the extent of Sicily
sicily_ext <- c(12, 16, 36.5, 38.5)
sicily <- crop(s, sicily_ext)

# make a time variable (to be used in regression)
time <- 1:nlayers(s)

# run the regression
fun <- function(x) { if (is.na(x[1])){ NA } else {lm(x ~ time)$coefficients[2] }} 
x2 <- calc(sicily, fun)

#predict to a raster
predicted <- predict(sicily, x2, progress='text')


The exact output is:

Loading required package: raster
Loading required package: sp
Loading required package: rgdal
rgdal: version: 1.0-6, (SVN revision 555)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 1.11.2, released 2015/02/10
 Path to GDAL shared files: /usr/share/gdal/1.11
 Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.1-1 

0%Error in UseMethod("predict") : no applicable method for 'predict' applied to an object of class "data.frame"

and the resulting x2 raster is shown here enter image description here

  • 3
    Can't say exactly without seeing your output but you are returning coefficients in your lm call and not a lm model object. Because of this, you are passing predict an invalid lm model object. You will have to write a function that applies the coefficients in a regression equation. Nov 20 '15 at 19:12
  • But how can I return the lm() model object?
    – jim
    Nov 21 '15 at 15:34
  • x <- lm(...); now x is an lm model object. Nov 21 '15 at 20:30
  • does not work: x <- lm(sicily ~ time) returns Error in model.frame.default(formula = sicily ~ time, drop.unused.levels = TRUE) : object is not a matrix
    – jim
    Nov 21 '15 at 21:58
  • If you want to fit linear models to spatial (or spatio-temporal) data in R, you should look into more specialised packages for fitting GLM (and GAM, for nonlinear models). Both INLA and Stan have good tutorials on how to do this.
    – Matt
    Dec 11 '15 at 22:41

You cannot use predict here. predict would work if you had a global model and you wanted to compute local (grid cell) level predictions. But it appears that you want to fit local models (a different model for each grid cell) and make predictions with these.

With these steps

fun <- function(x) { if (is.na(x[1])){ NA } else {lm(x ~ time)$coefficients[2] }} 
x2 <- calc(sicily, fun)

You got the slope. You could get both the slope and the intercept and than manually compute the values like this:

# get all coefficients
bfun <- function(x) { if (is.na(x[1])){ c(NA, NA) } else {lm(x ~ time)$coefficients}} 
x <- calc(sicily, bfun)
p1 <- x[[1]] + x[[2]] * time

Or in one step

pfun <- function(x) { if (is.na(x[1])){ rep(NA, length(x)) } else { predict( lm(x ~ time))}} 
p2 <- calc(sicily, pfun)
  • What you propose is basically identical to what @fdetsch said i.e. manually use of lm() output to compute y=a + bx. I believe that if I have a time series of rasters say (R1, R2, R3) and want to predict the next year say R4 the function should be R4 = a + bR3. How would we compute the global model for lm()?
    – jim
    Nov 23 '15 at 13:24
  • In principle, there are many ways to create a single, global, model. And then you could do predict(sicily, model). But it does not seem to be what you are after, conceptually. But, for example, you could do something along these lines v <- values(sicily); m <- lm(x ~ y + z , data=data.frame(v)); predict(sicily, m). Nov 24 '15 at 18:52
  • I am interested in per pixel regression but also to see what are the possibilities. I tried v <- values(sicily); m <- lm(sicily ~ time , data=data.frame(v)) but I get the same error ('object is not a matrix'). time and sicily are loaded as shown in the original code above.
    – jim
    Nov 25 '15 at 6:59

As long as it's linear regression you are dealing with, you could simply create two raster layers - one for slope, one for intercept - and do the prediction manually (i.e., y = m*x + t). Here is some sample code using GIMMS NDVI data downloaded and processed using the gimms package. You should be able to replace it with whatever kind of data you are dealing with.

## load 'gimms' package

## acquire sample data
gimms_envi <- downloadGimms(as.Date("2000-01-01"), as.Date("2000-06-30"))
gimms_tifs <- rasterizeGimms(gimms_envi, remove_header = TRUE, 
                             filename = gimms_envi, format = "GTiff")

## create spatial subset (sicily)
sicily_ext <- extent(c(12, 16, 36.5, 38.5))
gimms_tifs <- crop(gimms_tifs, sicily_ext)

## create time information
timestamps <- 1:nlayers(gimms_tifs)

## slope and intercept
lm_intercept <- calc(gimms_tifs, fun = function(x) {
  if (all(is.na(x)))
    return(coef(lm(x ~ timestamps))[1])

lm_slope <- calc(gimms_tifs, fun = function(x) {
  if (all(is.na(x)))
    return(coef(lm(x ~ timestamps))[2])

## predict ndvi (y = m*x + t)
lm_slope * gimms_tifs + lm_intercept
  • I adapted your code to my data and it works for the linear case. I use sicily[[3]] because my time series contains at the moment 3 rasters (say R1, R2, R3) and the predicted is say R4. So R4 has to be based on R3 I believe i.e. R4 = intersect + slope*R3 (note R3 not the complete stack). Since the objective is to also use non linear models, I will keep looking for a way to make predict() work. Thanks.
    – jim
    Nov 22 '15 at 9:04

An old topic. But could be useful for someone. Basic example:


set.seed(321) # make this example fully reproducible

r <- raster()

res(r) <- c(5,5)

r_list <- list()

time <- seq(8)

for (i in time) {
  r_list[[i]] <- setValues(x = r, values = rnorm(n = ncell(r), mean = i, sd = 0.1))

sicily <- stack(r_list)

Is the same approach than @RobertH, but adding newdata argument in prediction:

# pixel-by-pixel approach

pfun <- function(x) { if (is.na(x[1])){ rep(NA, length(x)) } else { predict( lm(x ~ time), newdata = data.frame(time = c(9:12)))}} 
p2 <- calc(sicily, pfun)


enter image description here

And the aim of @jim, searching a global fit:

# global model


df <- as.data.frame(sicily)
names(df) <- time

df %>% 
  gather(key = 'time', value = 'x') %>%
  mutate(time = as.numeric(time)) %>%
  lm(x ~ time, data = .) %>% 
  predict(newdata = data.frame(time = c(9:12)))

##         1         2         3         4 
##  8.996832  9.996315 10.995799 11.995282 

But, the output is only one value per time step

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.