0
$\begingroup$

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')
$\endgroup$

1 Answer 1

0
$\begingroup$

According to here: "The names in the Raster object should exactly match those expected by the model."

This means that even if I used coarse resolution data to perform the LR, I could use the predict function to apply the model parameters to a finer spatial scale, as long my SpatRaster has the same name as my input data.

library(terra)

wd = "path/"

ntl = rast(paste0(wd, "ntl.tif"))

tirs = rast(paste0(wd, "tirs.tif"))
tirs_res <- resample(tirs, 
                     ntl, 
                     method = "bilinear")

s = c(ntl, tirs_res)

names(s)

model <- lm(formula = ntl ~ B10_median, data = s, na.action = na.omit) # B10_median is how the tirs band was named when downloaded using Google Earth Engine
summary(model)

p = predict(tirs, model)
res(p)
$\endgroup$
1
  • 1
    $\begingroup$ Your answer could be improved with additional supporting information. Please edit to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center. $\endgroup$
    – Community Bot
    Commented Dec 11, 2022 at 13:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.