Search code examples
rlinear-regressionrastercoefficients

Spatial linear regression: apply regression's model parameters at another spatial scale


I have two raster layers and I want to perform linear regression (LR). The first raster has 500m pixel size (called ntl.tif) and the other has 100m (called tirs.tif). In order to perform LR I need to aggregate the tirs to match the spatial resolution of ntl. After that I can perform LR and I can use the predict function at the coarse spatial scale.

My question is this: How can I apply the model parameters (intercept and slope) to predict at the fine spatial scale? What I mean is, I want to use the function predict so I can create the lm_pred raster (see code below) and not by typing the coefficients manually like I do. I know I need to check that the names (band names) of the raster from which the model was fitted should be the same as the names of the raster to which the model will be applied, but I can't figure it out how I can do that: This is what I am doing so far:

library(terra)

ntl = rast("path/ntl.tif") # coarse resolution raster

tirs = rast("path/tirs.tif") # fine resolution raster

tirs_res <- resample(tirs, ntl, method="bilinear")

s = c(ntl, tirs_res)

names(s) = c("ntl", "tirs")

model <- lm(formula = ntl ~ tirs, data = s)

lm_pred = model$coefficients[1] + model$coefficients[2] * tirs

If I use the predict function:

p = predict(tirs, model)

I am getting this error: Error in model.frame.default(Terms, newdata, na.action = na.action, xlev = object$xlevels) : object is not a matrix. In addition: Warning message:'newdata' had 1377 rows but variables found have 1008 rows

Here's my dataset:

ntl = rast(ncols=272, nrows=200, nlyrs=1, xmin=12662503.7366, xmax=12798503.7366, ymin=3532049.3009, ymax=3632049.3009, names=c('ntl'), crs='EPSG:3857')

tirs = rast(ncols=1377, nrows=1008, nlyrs=1, xmin=12662000, xmax=12799700, ymin=3531700, ymax=3632500, names=c('B10_median'), crs='EPSG:3857')

Solution

  • You can use predict to apply that fitted model to any raster you want. You just need to check that the names (band names) of the raster from which the model was fitted are the same as the names of the raster to which the model will be applied.

    library(terra)
    lm_pred010 <- predict(pop, model)