Home > Net >  R: apply regression's model parameters at another spatial scale
R: apply regression's model parameters at another spatial scale

Time:12-03

I have three raster layers, two coarse resolution and one fine resolution. I perform simple linear regression using the coarse resolution rasters and then I extract the coefficients of the regression (slope and intercept) and I use them with my fine resolution raster. So far, I am doing this manually by plotting the coefficients in the console. Is there any way I can do this automatically (somehow to use the coefficients without having to print them and copy-paste them, because I have 10-20 rasters to do this).

Here is what I am doing so far:

library(terra)
library(sp)

ntl = rast("path/ntl.tif")
vals_ntl <- as.data.frame(values(ntl))
ntl_coords = as.data.frame(xyFromCell(ntl, 1:ncell(ntl)))
combine <- as.data.frame(cbind(ntl_coords,vals_ntl))

ebbi = rast("path/tirs010.tif")
ebbi <- resample(ebbi, ntl, method="bilinear")
vals_ebbi <- as.data.frame(values(ebbi))

s = c(ntl, ebbi)
names(s) = c('ntl', 'ebbi')

block.data <- as.data.frame(cbind(combine, vals_ebbi))
names(block.data)[3] <- "ntl"
names(block.data)[4] <- "ebbi"

block.data <- na.omit(block.data)

model <- lm(formula = ntl ~ ebbi, data = block.data)

#predict to a raster
summary(model)
model$coefficients # here I plot the coefficients into the console and then manually I use them below

pop = rast("path/pop.tif")
lm_pred010 = 19.0540153   0.2797187 * pop
plot(lm_pred010)
res(lm_pred010)

writeRaster(lm_pred010, 
            filename = "path/lm_pred010.tif",
            overwrite = T)

The rasters I am using:

ntl = rast(ncols=101, nrows=84, nlyrs=1, xmin=509634.6325, xmax=550034.6325, ymin=161998.158, ymax=195598.158, names=c('ntl'), crs='EPSG:27700')

ebbi = rast(ncols=101, nrows=84, nlyrs=1, xmin=509634.6325, xmax=550034.6325, ymin=161998.158, ymax=195598.158, names=c('focal_sum'), crs='EPSG:27700')

pop = rast(ncols=407, nrows=343, nlyrs=1, xmin=509615.9028, xmax=550315.9028, ymin=161743.6035, ymax=196043.6035, names=c('pop'), crs='EPSG:27700')

CodePudding user response:

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)
  • Related