I have two raster layers and I want to perform linear regression (LR). The first raster has 500m pixel size and the other has 100m. In order to perform LR I need to aggregate the fine resolution raster to match the spatial resolution of coarse one. 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
.
Another error occurs when I convert the fine resolution raster to a data.frame
:
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 92009 rows but variables found have 344 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')