r-在栅格单元级别跨光栅文件的层应用函数

Question

我有一个具有多层的netcdf光栅对象(时间：30年*12个月=360步(。

光栅由多个网格单元(例如1000个(组成。我想分别对每个月和每个网格单元进行线性去预测，以获得与输入具有相同尺寸的残差光栅对象。

对于列[x，y，value，month，year]的长格式数据帧，我会这样做：

library(tidyverse)

data %>%
group_by(x, y, month) %>%
nest() %>%
mutate(model_trend = map(data,  ~ lm(value ~ year, data = .))) %>%
mutate(augment_trend  = map(model_trend,  ~ broom::augment(.))) %>%
unnest(cols = c(augment_trend)) %>%
ungroup() %>%
dplyr::select(x, y, year, month, .resid)

我相信，使用空间包(如terra, stars等(，有更高效、同样简单的方法可以完成此操作。你能给我一些建议吗？

感谢

Answer 1

我不太确定你在做什么，因为你没有提供数据，也没有使用标准R。但以下是如何使用线性回归来去除SpatRast中的细胞。

library(terra)
s <- rast(system.file("ex/logo.tif", package="terra"))   
n <- 1:nlyr(s)
f <- function(y) {
residuals( lm(y ~ n) )
}
a <- app(s, f)

还有一种更快的线性代数方法(可能不适用于您的数据(

X <- cbind(1, n)
invXtX <- solve(t(X) %*% X) %*% t(X)
qfun <- function(y) {
p <- invXtX %*% y
y - (p[1] + p[2] * n)
}
app(s, qfun) 
#class       : SpatRaster 
#dimensions  : 77, 101, 3  (nrow, ncol, nlyr)
#resolution  : 1, 1  (x, y)
#extent      : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source      : memory 
#names       :       red,      green,      blue 
#min values  : -3.833333, -17.000000, -3.833333 
#max values  :  8.500000,   7.666667,  8.500000 

通过使用tapp而不是app，您可以在所有月份使用以上内容