在 R 中的 RasterStack 上运行 PCA

Question

我正在尝试对 Bioclim 中可用的一组变量运行 PCA。我注意到 rasterPCA() 在 R 中不再可用，因为它的软件包已停止使用，看起来这是因为该软件包的一些问题尚未解决。下面，我尝试直接在 RasterStack 上运行 PCA（在本例中这是必要的，因为如果我转换为数据帧，我将丢失每行的空间信息，我的目标是创建一个新的 rasterstack，其中包含所有此 PCA 将生成的 PC）。

#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'

有谁知道如何使这种直接 PCA 分析工作，类似于 rasterPCA 过去的工作方式？或者在运行 PCA 时保留空间数据的方法？

Answer 1

我正在使用“terra”，替代“raster”。该示例来自

?terra::predict

，但您可以使用栅格执行类似的操作。

示例数据（在提出 R 问题时始终包含一些数据）

library(terra)
logo <- rast(system.file("ex/logo.tif", package="terra"))   

解决方案：

pca <- prcomp(logo)

pca
#Standard deviations (1, .., p=3):
#[1] 124.814772  17.084151   1.456423
#
#Rotation (n x k) = (3 x 3):
#             PC1        PC2        PC3
#red   -0.5935507  0.5073220 -0.6247575
#green -0.5846008  0.2617406  0.7679412
#blue  -0.5531179 -0.8210458 -0.1412245

您可以使用

predict

来映射主成分：

x <- predict(logo, pca)
plot(x)

如果栅格非常大，您可以使用采样来避免内存不足

sr <- spatSample(logo, 50000, "regular")
pca <- prcomp(sr)
# etc

要从 RasterStack 中创建 SpatRaster

s

，您可以执行

r <- rast(s)

并使用您可以执行的其他方式

ss <- stack(r)

Answer 2

crop_bioclim

可能不是

prcomp()

函数可接受的格式。尝试先将其转换为

data frame object

并执行

pca

:

...
crop_bioclim <- crop(x = bioclim_all, y = bounding_box)

# convert to data frame
crop_bioclim_df <- as.data.frame(crop_bioclim)

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

Answer 3

我也遇到了同样的问题。解决方案分为两部分。一，我找到了 rasterPCA() 的替代品，形式为 raster_pca() (package 'synoptReg')。但是，它仍然会给我以下错误：

# Error in (function (classes, fdef, mtable)  : 
# unable to find an inherited method for function ‘nlayers’ for signature ‘"SpatRaster"’

看起来 PCA 函数不能与 SpatRaster 对象一起使用，但它们可以与 RasterStack 对象一起使用（也许其他人可以解释为什么会出现这种情况）。

使用 worldclim_global() 获取 Worldclim 数据将创建 SpatRaster，rast() 也将创建。

制作一个普通的旧 RasterStack 需要一些额外的工作，但这段代码对我来说很实用：

library(terra)

bio1<-raster('C:/.../wc2.1_2.5m/wc2.1_2.5m_bio_1.tif')
bio2<-raster('C:/.../wc2.1_2.5m/wc2.1_2.5m_bio_2.tif')
bio3<-raster('C:/.../wc2.1_2.5m/wc2.1_2.5m_bio_3.tif')
bio4<-raster('C:/.../wc2.1_2.5m/wc2.1_2.5m_bio_4.tif')

biostack<-stack(bio1,bio2,bio3,bio4)

# pca1<-prcomp(biostack)
## Gives: Error in svd(x, nu = 0, nv = k) : infinite or missing values in 'x'
## clearly base R stats is just useless for this

#install.packages('synoptReg')
library(synoptReg)# raster_pca function

pca1<-raster_pca(biostack,aggregate=0,focal=0)

pca1

# $rasterPCA
# class      : RasterStack 
# dimensions : 4320, 8640, 37324800, 4  (nrow, ncol, ncell, nlayers)
# resolution : 0.04166667, 0.04166667  (x, y)
# extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# crs        : +proj=longlat +datum=WGS84 +no_defs 
# names      :       PC1,       PC2,       PC3,       PC4 
# min values : -2.656944, -4.054211, -2.162238, -2.014184 
# max values :  3.999249,  3.174859,  2.088713,  0.939421 
# $summary
#                  Comp.1    Comp.2     Comp.3     Comp.4
# sdev          1.6386658 0.9594545 0.54523932 0.31134434
# prop.variance 0.6713064 0.2301383 0.07432148 0.02423382
# cum.variance  0.6713064 0.9014447 0.97576618 1.00000000