I am trying to create a raster of the residuals of a regression between two rasters. i.e. I would like to carry out a regression of one raster against another of the same extent and plot the residual for each cell of the raster at the same resolution and extent as the original rasters.

It seems the GRASS 7.0 has a function called r.regression.multi which would let me carry this out but I have zero experience in GRASS and I can't even get GRASS 7.0 to work (GRASS 6.4 works fine but doesn't have r.regression.multi so is useless!)

I currently have the files in almost every format possible (Idrisi raster, TIFF, ASCII .grd in R) any help would be much appreciated!



Using R, once you have read in or created rasters x and y, the regression is performed and a raster of residuals r can be produced with three commands:

model <- lm(getValues(y) ~ getValues(x), na.action=na.exclude)
r <- y
r[] <- residuals.lm(model)

(To verify everything is working correctly, printing and plotting the results is a good idea.) The workflow generalizes readily to multiple independent rasters and to other models such as generalized linear models, nonlinear fits, models with spatial correlation, etc.

Worked example

Let's create some sample data:

# Create raster objects.
nRows <- 7
nCols <- 10
x <- raster(nrow=nRows, ncol=nCols)
y <- x
e <- x
# Simulate a bivariate relation y ~ x with error e.
e[] <- rnorm(nRows*nCols)
x[] <- 1:nCols
y <- x * 0.5 + 5.0 + e

Thus, the independent raster is x and the dependent is y; the errors in the linear relation are e. (In practice the errors are never known--they have to be estimated by the residuals--but with a simulated dataset we can hold on to the error values for checking things later.)

Do the regression:

model <- lm(getValues(y) ~ getValues(x), na.action=na.exclude)

The output lets us check that everything is working as expected. Here's part of it:

             Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4.92223    0.25842   19.05   <2e-16 ***
getValues(x)  0.53734    0.04165   12.90   <2e-16 ***
Residual standard error: 1.001 on 68 degrees of freedom

The estimates are highly significant and quite close to the actual values (4.92 vs. 5.00 and 0.54 vs. 0.05). The residual standard error of 1.001 closely estimates the error built into the simulation, 1.0. Finally, 68 = 7*10 - 2 counts the cells minus the two estimated parameters. Everything looks good.

Let's create your grid of residuals and plot it:

r <- y # Create a grid object for the residuals
r[] <- residuals.lm(model)

Residual plot

As a double-check (because we have been copying values out of and then back into raster format), let's plot the residuals against the true errors. The points should closely follow the line y=x, with just a tiny bit of scatter:

plot(getValues(e), getValues(r), xlab="Error grid", ylab="Residual grid")


Everything is fine.

  • You make my efforts below look so...quaint. Nice job. The bit I included about managing the NA values one usually gets in raster may be helpful though. If, as I expect, your response gets booted to the top, some folks may want to look below for that bit of code. – csfowler Feb 10 '12 at 17:36
  • 1
    Thanks for pointing out the handling of NAs. That's useful to know. I made appropriate changes in my reply (and tested them by setting some values of y to NA before doing the regression, as in y[c(13,14,58)]<-NA). These occur at two places: in the na.action option of the lm() call and in using the residuals.lm method to extract the residuals (rather than accessing them directly). BTW, you control the sequence of answers: look for the "Active Oldest Votes" links at the bottom right of the question. – whuber Feb 10 '12 at 18:53
  • thanks very much for this, very helpful! Although I have one issue in that the dependent variable raster does not have any NAs in it whereas the independent variable raster has the NAs that I require - is it possible to copy the NAs across? – M Hudson Feb 14 '12 at 10:29
  • I don't think you need to (in my testing, I put NAs into the DV raster only), but if you wish you can do something like y[is.na(x)]<-NA. This can be extended to multiple DVs via y[is.na(x1) | is.na(x2)]<-NA. – whuber Feb 14 '12 at 15:51
  • 1
    The raster package can use GDAL, so it can read and write many raster formats. Even without GDAL it can write several basic formats, including the ESRI ASCII export format. See the help for writeRaster. – whuber Feb 15 '12 at 15:05

If you already have it in R's .grd format why not do it in R with the raster package

r <- raster(system.file("external/test.grd", package="raster"))
plot(r) #get an initial view
pt.1<-getValues(r) #extract values from raster into a vector
summary(pt.1) #what do our values look like
pt.2<-runif(9200,min=min(128),max=max(1806)) #generate something to compare with
pt.2[which(is.na(pt.1))]<-NA #presumably your two rasters would already have the
                             #the same NA values
df<-data.frame(pt.1,pt.2)    #combine in data frame

model<-lm(pt.1~pt.2,data=df,na.action=na.exclude) #model, na.exclude padds the
                             #results so we retain the position of the residuals
                             #in our original vector
summary(model) #sadly my model has lousy fit, such is the life of random numbers

residuals<-pt.2 #create a vector to hold our residuals with NA's in place

#this is the cool part. na.exclude padded our residuals, we can access the
#original slot through the 'names' function and insert our residuals where they 
#belong leaving the NA's in place
resid<-r  #create a new raster to hold our results
values(resid)<-residuals #replace our values in this new raster
plot(resid) #random error never looked so good.
  • +1 for pointing out that missing data (NA) may need special handling. – whuber Feb 10 '12 at 19:07

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.