我有几个想要在 PCA 中处理的大型栅格（以生成汇总栅格）。 我见过几个例子，人们似乎只是简单地调用 prcomp 或 princomp。但是，当我这样做时，我收到以下错误消息：
Error in as.vector(data): no method for coercing this S4 class to a vector
示例代码：
files<-list.files() # a set of rasters
layers<-stack(files) # using the raster package
pca<-prcomp(layers)
我尝试过使用光栅砖而不是堆栈，但这似乎不是问题。我需要提供什么方法才能将栅格数据转换为矢量格式？我知道有一些方法可以对栅格进行采样并从中运行 PCA，但我真的很想了解为什么上述方法不起作用。
谢谢！
【问题讨论】：
上面的方法不起作用仅仅是因为 prcomp 不知道如何处理光栅对象。它只知道如何处理向量，并且强制到向量不起作用，因此出错。
您需要做的是将每个文件读入一个向量，然后将每个栅格放入矩阵的一列中。然后，每一行将是单个空间位置的时间序列值，每一列将是某个时间步长的所有像素。请注意，这种方法不需要精确的空间坐标。这个矩阵作为prcomp的输入。
可以使用readGDAL 读取文件，并使用as.data.frame 将空间数据转换为data.frame。
【讨论】：
这是一个可行的解决方案：
library(raster)
filename <- system.file("external/rlogo.grd", package="raster")
r1 <- stack(filename)
pca<-princomp(r1[], cor=T)
res<-predict(pca,r1[])
显示结果：
r2 <- raster(filename)
r2[]<-res[,1]
plot(r2)
【讨论】：
回答我自己的问题：我最终做了一些稍微不同的事情：我没有使用每个栅格单元作为输入（非常大的数据集），而是取了一个点样本，运行 PCA，然后保存输出模型，以便我可以对每个网格单元进行预测……也许不是最好的解决方案，但它确实有效：
rasters <- stack(myRasters)
sr <- sampleRandom(rasters, 5000) # sample 5000 random grid cells
# run PCA on random sample with correlation matrix
# retx=FALSE means don't save PCA scores
pca <- prcomp(sr, scale=TRUE, retx=FALSE)
# write PCA model to file
dput(pca, file=paste("./climate/", name, "/", name, "_pca.csv", sep=""))
x <- predict(rasters, pca, index=1:6) # create new rasters based on PCA predictions
【讨论】：
RStoolbox包http://bleutner.github.io/RStoolbox/rstbx-docu/rasterPCA.html中有rasterPCA函数
例如：
library('raster')
library('RStoolbox')
rasters <- stack(myRasters)
pca1 <- rasterPCA(rasters)
pca2 <- rasterPCA(rasters, nSamples = 5000) # sample 5000 random grid cells
pca3 <- rasterPCA(rasters, norm = FALSE) # without normalization
【讨论】：
另一种选择是从光栅堆栈中提取值，即：
rasters <- stack(my_rasters)
values <- getValues(rasters)
pca <- prcomp(values, scale = TRUE)
【讨论】：
这是另一种扩展 @Daniel 提出的 getValues 方法的方法。结果是一个栅格堆栈。索引 (idx) 引用非 NA 位置，以便考虑 NA 值。
library(raster)
r <- stack(system.file("external/rlogo.grd", package="raster"))
r.val <- getValues(r)
idx <- which(!is.na(r.val))
pca <- princomp(r.val, cor=T)
ncomp <- 2 # first two principle components
r.pca <- r[[1:ncomp]]
for(i in 1:ncomp) { r.pca[[i]][idx] <- pca$scores[,i] }
plot(r.pca)
【讨论】：