在光栅砖中使用 PCA 降维答案

【问题标题】：dimensionality reduction with PCA in raster brick在光栅砖中使用 PCA 降维
【发布时间】：2019-03-18 10:30:36
【问题描述】：

基于此处的示例：

How to perform dimensionality reduction with PCA in R

和

How to reverse PCA and reconstruct original variables from several principal components?

我正在尝试在光栅砖（69 层）上执行 PCA，然后获取领先的 PC，最后仅使用累积比例约为 95% 的 PC 重建原始变量。

library(raster)
library(ncdf4)

ln <- "https://www.dropbox.com/s/d88iuvp9oio14zk/test.nc?dl=1" # ~400 kb size

### DOWNLOAD THE FILE
download.file(ln,
              destfile="test.nc",
              method="auto")

st <- brick("test.nc")
nlayers(st)

### DO THE PCA
pca <- prcomp(st[]) 
# to visualize pcs as rasters
x <- predict(st, pca, index=1:4)
spplot(x) # there are the first 4 PCs explaining most of the data.

然后我尝试从前 4 台 PC 中重建原始变量，因为我对这些的空间分布感兴趣：

### PCA DETAILS
summary(pca) # importance of components
plot (pca) # scree plot
loadings(pca) #eigens 

mu <- colMeans(as.matrix(st)) # get the column means to use after

#### REDUCTION
nComp <- 4
Xhat <- pca$x[,1:nComp] %*% t(pca$rotation[,1:nComp])
Xhat <- scale(Xhat, center = -mu, scale = FALSE)

在这里，我以为我只会得到前 4 台 PC。但是，我还是像以前一样以 69 结束：

### CHECK THE DIMENSIONS
dim(Xhat)

### THEN CREATE THE RASTER WITH THE PCs
coords <- coordinates(st[[1]]) # get the lon/lat

rst <- cbind(coords, Xhat) # bind the coordinates 
rst <- rasterFromXYZ(rst) # create the raster
plot(rst)

我在这里错过了什么？我不是 PCA 方面的专家，但最初的想法是让更少的层能够解释原始数据中的模式。谢谢！

【问题讨论】：

标签： r pca r-raster svd

【解决方案1】：

在这里提问时，请不要指向 Dropbox 上的文件，而是包括一些示例数据，如下所示：

library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
s <- stack(b, flip(b, "y"), setValues(raster(b), runif(ncell(b))))
names(s) <- paste0("var", 1:nlayers(s))    

pca <- prcomp(values(s)) 
x <- predict(s, pca, index=1:4)

您通过子集 PC 创建 Xhat，但 pca$rotation 包含所有变量

round(pca$rotation[,1:nComp],1)
      PC1  PC2  PC3  PC4
var1 -0.4  0.4 -0.4  0.4
var2 -0.4  0.4 -0.2  0.2
var3 -0.4  0.4  0.6 -0.6
var4 -0.4 -0.4 -0.4 -0.4
var5 -0.4 -0.4 -0.2 -0.2
var6 -0.4 -0.4  0.6  0.6
var7  0.0  0.0  0.0  0.0

这是有道理的，因为您说您想要“从前 4 台 PC 重建原始变量，因为我对这些的空间分布感兴趣”。所有变量都对 PC 有贡献。

你可能真正想要的是plot(x)？

您使用糟糕的代码从 Xhat 创建 RasterBrick。你可以这样做：

Xhat <- pca$x[,1:nComp] %*% t(pca$rotation[,1:nComp])
Xhat <- scale(Xhat, scale = FALSE)

b <- brick(s, values=FALSE)
b <- setValues(b, Xhat)
b
#class      : RasterBrick 
#dimensions : 77, 101, 7777, 7  (nrow, ncol, ncell, nlayers)
#resolution : 1, 1  (x, y)
#extent     : 0, 101, 0, 77  (xmin, xmax, ymin, ymax)
#crs        : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
#source     : memory
#names      :          var1,          var2,          var3,          var4,          var5,          var6,          var7 
#min values : -184.32344039, -184.48714657, -193.05823803, -184.32341010, -184.48718831, -193.05827512,   -0.01466663 
#max values :   73.33872354,   70.38724578,   63.48912986,   73.33875039,   70.38723822,   63.48906605,    0.01193009

比较 b 和 s

m <- cellStats(s, mean)
bb <- b + m
plot(bb, s)

【讨论】：