there are a number of posts (1, 2) showing how to conduct OLS regression on raster stacks in R - but I can't find anything about using GLS regression with a specified spatial correlation structure.

Some code below

library(nlme) library(raster) # example raster ras <- brick(nrow = 72, ncol = 144, nl = 5) values(ras) <- runif(ncell(ras) * nlayers(ras), -1, 1) time <- ts(seq(1:5)) ras <- setZ(ras, time, "Time") ras # define the function to calculate/extract lm slope calculateSlope <- function(x) {if (is.na(x[1])) {NA} else{ lm(x ~ time)$coefficients[2] }} # lm slope slopeRaster <- calc(ras, fun = calculateSlope) slopeRaster # define the function to calculate/extract gls slope calculateGLSSlope <- function(x) {if (is.na(x[1])) {NA} else{ nlme::gls(x ~ time)$coefficients[2] }} # GLS slope slopeRasterGLS <- calc(ras, fun = calculateGLSSlope) slopeRasterGLS # use all.equal to allow for small differences (sqrt(.Machine$double.eps)) # gls and ols are equal when no correlation structre is defined all.equal(values(slopeRaster), values(slopeRasterGLS)) [1] TRUE

If I then re-write the slope function to something like

calculateGLSSlope <- function(x) {if (is.na(x[1])) {NA} else{ nlme::gls(x ~ time, correlation = corSpher(form = ~ time)))$coefficients[2] }}

but I'm not 100 % sure if this accounts for the X and Y coords of the data, only the correlation between time. How is it possible to specify the x~y coords of each of the raster cells when calc works on a pixel-by-pixel basis? I thought something along the lines of

calculateSlope <- function(j) {if (is.na(j[1])) {NA} else{ coords <- xyFromCell(ras, j) nlme::gls(j ~ Year, correlation = corSpher(form = ~ coords$x + coords$y + time))$coefficients[2] }}

But in this case j is represented by the cell-value and there is no information about the cell location, so xyFromCell fails

Edit 11/05/2017

Following on from @geoCompR answer below, I've attempted to follow their steps to get the lat/long coordinates in a model. The additions to the code above are below.

# create rasters for latitude and longitude rLat <- raster(ncols = ncol(ras), nrows = nrow(ras)) values(rLat) <- values(init(ras, "y")) rLon <- raster(ncols = ncol(ras), nrows = nrow(ras)) values(rLon) <- values(init(ras, "x")) # Stack the original raster with the lon and lat rasters newStack <- stack(ras, rLon, rLat) newStack # define the function to calculate/extract lm slope with coords as var calculateSlopeWCoords <- function(x) {if (is.na(x[1])) {NA} else{ # layers x[1:5] is the vector of dependent variables, time is the time vector from above # x[6] is the longitude, x[7] is the latitude lm(x[1:5] ~ time + x[6] + x[7])$coefficients[2] }} slopeRaster2 <- calc(newStack, fun = calculateSlopeWCoords)

However, when attempting to run slopeRaster2 I receive the following error

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

How can I include the lat and long variables in a function that can be passed to calc?

If I try to create a raster stack where I have 5 layers of dependent variables, 5 layers for time, and then 5 each for latitude and longitude; I could then write a function like this

slopeFun = function(x) { # In case of any NA return four NA values for coefficient estimates # Four variables = four coefficients (1 intercept, 3 indep) if (any(is.na(x))) return(c(NA, NA, NA, NA)) model <- lm(x[1:5] ~ x[6:10] + x[11:15] + x[16:20], na.action = na.omit) coefficients(model)}

This can be passed to calc like so: slopeRas <- calc(rasStack, fun = slopeFun) , which returns a raster stack of model coefficients, but the coefficients for the Latitude and Longitude layers all come back as NA as each of the layers for latitude, and each of the layers for longitude, are perfectly correlated with themselves and therefore can't be modelled with lm . The code below shows this approach.

library(raster) library(nlme) # Generate an example raster brick set.seed(1) ras <- raster(res = 2.5) brRas <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5) values(brRas) <- rnorm(ncell(brRas) * nlayers(brRas), 0, 1) brRas # generate a brick to represent time timeRas <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5) for (r in 1:nlayers(timeRas)){ values(timeRas[[r]]) <- seq(1:5)[r] } timeRas time <- ts(seq(1:5)) # Generate Lat/Long raster bricks lat <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5) values(lat) <- values(init(brRas, "y")) long <- brick(ncols = ncol(ras), nrows = nrow(ras), nl = 5) values(long) <- values(init(brRas, "x")) # Set the layer names so they are easily identifiable names(brRas) <- paste0("Layer", seq(1:5)) names(timeRas) <- paste0("Time", seq(1:5)) names(lat) <- paste0("Lat", seq(1:5)) names(long)<- paste0("Long", seq(1:5)) # Stack it all together and check rasStack <- stack(brRas, timeRas, long, lat) rasStack names(rasStack) # Define the function for the coefficients slopeFun <- function(x) { # In case of any NA return four NA values for coefficient estimates # Four variables in the model = four coefficients (1 intercept, 3 indep) if (any(is.na(x))) return(c(NA, NA, NA, NA)) model <- lm(x[1:5] ~ x[6:10] + x[11:15] + x[16:20], na.action = na.omit) coefficients(model)} # Pass the function to calc slopeRas <- calc(rasStack, fun = slopeFun) slopeRas class : RasterBrick dimensions : 72, 144, 10368, 4 (nrow, ncol, ncell, nlayers) resolution : 2.5, 2.5 (x, y) extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax) coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 data source : in memory names : X.Intercept., x.6.10., x.11.15., x.16.20. min values : -4.236375, -1.066759, NA, NA max values : 4.112934, 1.228611, NA, NA

The only other thing I can think of doing, is creating a dataframe of my raster stack containing xy centroids, the value of the dependent variable, the z-index of the raster (time), and the CellID, and then fitting the models row-wise, and then summarising the results by CellID to avoid the perfect collinearity issue.