在R中对RasterStack执行PCA遇报错，求rasterPCA替代方案

问题描述 #

尝试对Bioclim变量集运行PCA，因rasterPCA所属包停用无法使用，直接对RasterStack调用prcomp时出现报错，目标是生成包含所有主成分的新RasterStack且不丢失空间信息：

用户代码：

# Libraries:
library(geodata)
library(raster)

# Downloading the data:
bioclim_all <- worldclim_global(var = "bio", 
                                res = 0.5, 
                                path = "/data")

# Creating a bounding box:
bounding_box <- extent(x = c(-118.2724, -86.4236, 14.3237, 32.4306))

# Cropping to a smaller resolution:
crop_bioclim <- crop(x = bioclim_all, y = bounding_box)

# Conduct a PCA with standardization:
pca <- prcomp(crop_bioclim, center = TRUE, scale = TRUE)

报错信息：

Error in svd(x, nu = 0, nv = k) : infinite or missing values in 'x'

解决方案 #

1. 处理栅格缺失值 #

报错根源是栅格数据中存在缺失值（NA），需先清理数据：

# 提取所有无缺失值的栅格单元格数据
valid_values <- na.omit(getValues(crop_bioclim))

2. 执行PCA分析 #

基于清理后的数据集运行标准化PCA：

pca <- prcomp(valid_values, center = TRUE, scale. = TRUE)

3. 将PCA结果映射回栅格（保留空间信息） #

创建与原栅格空间属性一致的空栅格栈，将PCA结果填充到对应位置，同时保留原缺失值的位置：

# 创建空栅格栈，层数等于主成分数量
pc_stack <- raster(crop_bioclim, nlayers = ncol(pca$x))
names(pc_stack) <- paste0("PC", 1:ncol(pca$x))

# 获取原栅格中有效单元格的索引
valid_indices <- which(!is.na(getValues(crop_bioclim[[1]])))

# 逐个主成分填充数据
for(i in seq_len(ncol(pca$x))){
  pc_stack[[i]][valid_indices] <- pca$x[, i]
}

更简便的替代方案：使用RStoolbox的rasterPCA #

若可以安装RStoolbox包，其内置的rasterPCA功能可直接处理RasterStack，自动处理缺失值并返回带空间信息的主成分栅格栈：

# 安装并加载RStoolbox
install.packages("RStoolbox")
library(RStoolbox)

# 运行PCA，scale参数控制是否标准化
pca_result <- rasterPCA(crop_bioclim, scale = TRUE)

# 提取主成分栅格栈
pc_stack <- pca_result$map

内容的提问来源于stack exchange，提问作者Cameron