R:线性模型的光栅预测问题

Question

我使用函数raster::predict将线性模型的预测部分提取为光栅，但我得到了以下错误：

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 622 rows but variables found have 91 rows

我的数据集是两个卫星图像的RasterStack(相同的CRS和数据类型(。我找到了这个问题，但我无法解决我的问题。这是代码和数据：

library(raster)
ntl = raster ("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 = raster ("path/ebbi.tif")
ebbi <- resample(ebbi, ntl, method = "bilinear")
vals_ebbi <- as.data.frame(values(ebbi))
s = stack(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
r1 <- raster::predict(s, model, progress = 'text', na.rm = T)
plot(r1)
writeRaster(r1, filename = "path/lm_predict.tif")

数据可以从这里下载(我不知道共享一个较小的数据集是否仍然存在问题，所以我决定共享完整的数据集，当使用dput命令复制粘贴时，这个数据集相当大(

Answer 1

dput通常对空间数据不太有用，这是正确的；并且您应该避免使用它。但是，在大多数情况下，没有必要共享数据，因为您可以使用代码或附带R的数据创建示例数据，就像本网站上帮助文件和问题中的大多数示例一样。说"；我不知道通过共享一个较小的数据集，这个问题是否仍然存在；建议你应该做的第一件事就是找出答案。

如果你有一个想要复制的SpatRasterx，你可以从as.character(x)开始，这就是我所做的。

library(terra)
ntl <- rast(ncols=48, nrows=91, nlyrs=1, xmin=582360, xmax=604440, ymin=1005560, ymax=1047420, names=c('avg_rad'), crs='EPSG:7767')
ebbi <- rast(ncols=48, nrows=91, nlyrs=1, xmin=582360, xmax=604440, ymin=1005560, ymax=1047420, names=c('B6_median'), crs='EPSG:7767')
values(ntl) <- sample(100, ncell(ntl), replace=TRUE)
values(ebbi) <- runif(ncell(ebbi))

组合，设置名称，并将值放入data.frame中。对于较大的数据集，您可以使用spatSample(x, type="regular")进行采样。

x <- c(ntl, ebbi)
names(x) <- c("ntl", "ebbi")  

安装模型。可以分两步完成

v <- as.data.frame(x, na.rm=TRUE)
model <- lm(ntl ~ ebbi, data=v)    

或者一步

model <- lm(ntl ~ ebbi, data=x)    

现在预测(如果要将光栅保存到磁盘，请设置文件名(。

p <- predict(x$ebbi, model, filename="")

predict的第一个(SpatRraster(参数的名称必须与模型中的名称相匹配，这一点非常重要。所以在这种情况下，你可以使用x$ebbi或x[[2]]，但如果你使用ebbi，你会得到一个神秘的错误消息

p <- predict(ebbi, model)
#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 48 rows but variables found have 91 rows 

除非你先做

names(ebbi) <- "ebbi"
p <- predict(ebbi, model)

Answer 2

备选方案，使用raster包，解决方案为：

library(raster)
ntl = raster ("path/ntl.tif")
ebbi = raster ("path/ebbi.tif")
ebbi <- resample(ebbi, ntl, method = "bilinear")
s = stack(ntl, ebbi)
names(s) = c('ntl', 'ebbi') # important step in order to run the predict function successfully
block.data = data.frame(na.omit(values(s)))
names(block.data) <- c('ntl', 'ebbi')
model <- lm(formula = ntl ~ ebbi, data = block.data)
#predict to a raster
r1 <- raster::predict(s, model, progress = 'text', na.rm = T)
plot(r1)
writeRaster(r1, filename = "path/lm_predict.tif")

我根据这篇帖子找到了答案。