“时间模式”解释了时间序列在所有网格中占主导地位的时间变化，它由PCA的主成分（PCs，多个时间序列）表示。在 R 中，最重要的 PC PC1 是prcomp(data)$x[,'PC1']。
“空间模式”解释了 PC 对某些变量（在您的例子中是地理）的依赖程度，它由 每个主成分的负载 表示。例如，对于 PC1，它是prcomp(data)$rotation[,'PC1']。
这是一个使用您的数据为 R 中的时空数据构建 PCA 并显示时间变化和空间异质性的示例。
首先，必须将数据转换为带有变量（空间网格）和观测值（yyyy-mm）的 data.frame。
加载和转换数据：
load('spei03_df.rdata')
str(spei03_df) # the time dimension is saved as names (in yyyy-mm format) in the list
lat <- spei03_df$lat # latitude of each values of data
lon <- spei03_df$lon # longitude
rainfall <- spei03_df
rainfall$lat <- NULL
rainfall$lon <- NULL
date <- names(rainfall)
rainfall <- t(as.data.frame(rainfall)) # columns are where the values belong, rows are the times
要了解数据，请在地图上绘制 1950 年 1 月的数据：
library(mapdata)
library(ggplot2) # for map drawing
drawing <- function(data, map, lonlim = c(-180,180), latlim = c(-90,90)) {
major.label.x = c("180", "150W", "120W", "90W", "60W", "30W", "0",
"30E", "60E", "90E", "120E", "150E", "180")
major.breaks.x <- seq(-180,180,by = 30)
minor.breaks.x <- seq(-180,180,by = 10)
major.label.y = c("90S","60S","30S","0","30N","60N","90N")
major.breaks.y <- seq(-90,90,by = 30)
minor.breaks.y <- seq(-90,90,by = 10)
panel.expand <- c(0,0)
drawing <- ggplot() +
geom_path(aes(x = long, y = lat, group = group), data = map) +
geom_tile(data = data, aes(x = lon, y = lat, fill = val), alpha = 0.3, height = 2) +
scale_fill_gradient(low = 'white', high = 'red') +
scale_x_continuous(breaks = major.breaks.x, minor_breaks = minor.breaks.x, labels = major.label.x,
expand = panel.expand,limits = lonlim) +
scale_y_continuous(breaks = major.breaks.y, minor_breaks = minor.breaks.y, labels = major.label.y,
expand = panel.expand, limits = latlim) +
theme(panel.grid = element_blank(), panel.background = element_blank(),
panel.border = element_rect(fill = NA, color = 'black'),
axis.ticks.length = unit(3,"mm"),
axis.title = element_text(size = 0),
legend.key.height = unit(1.5,"cm"))
return(drawing)
}
map.global <- fortify(map(fill=TRUE, plot=FALSE))
dat <- data.frame(lon = lon, lat = lat, val = rainfall["1950-01",])
sample_plot <- drawing(dat, map.global, lonlim = c(-180,180), c(-90,90))
ggsave("sample_plot.png", sample_plot,width = 6,height=4,units = "in",dpi = 600)
如上所示，所提供的链接给出的网格数据包括代表加拿大降雨量（某些类型的指数？）的值。
主成分分析：
PCArainfall <- prcomp(rainfall, scale = TRUE)
summaryPCArainfall <- summary(PCArainfall)
summaryPCArainfall$importance[,c(1,2)]
这表明前两台 PC 解释了降雨数据中 10.5% 和 9.2% 的方差。
我提取了前两台 PC 的负载和 PC 时间序列本身：
“空间格局”（载荷），显示趋势（PC1和PC2）强度的空间异质性。
loading.PC1 <- data.frame(lon=lon,lat=lat,val=PCArainfall$rotation[,'PC1'])
loading.PC2 <- data.frame(lon=lon,lat=lat,val=PCArainfall$rotation[,'PC2'])
drawing.loadingPC1 <- drawing(loading.PC1,map.global, lonlim = c(-180,-30), latlim = c(40,90)) + ggtitle("PC1")
drawing.loadingPC2 <- drawing(loading.PC2,map.global, lonlim = c(-180,-30), latlim = c(40,90)) + ggtitle("PC2")
ggsave("loading_PC1.png",drawing.loadingPC1,width = 6,height=4,units = "in",dpi = 600)
ggsave("loading_PC2.png",drawing.loadingPC2,width = 6,height=4,units = "in",dpi = 600)
“时间模式”，前两个 PC 时间序列，显示数据的主要时间趋势
library(xts)
PC1 <- ts(PCArainfall$x[,'PC1'],start=c(1950,1),end=c(2014,12),frequency = 12)
PC2 <- ts(PCArainfall$x[,'PC2'],start=c(1950,1),end=c(2014,12),frequency = 12)
png("PC-ts.png",width = 6,height = 4,res = 600,units = "in")
plot(as.xts(PC1),major.format = "%Y-%b", type = 'l', ylim = c(-100, 100), main = "PC") # the black one is PC1
lines(as.xts(PC2),col='blue',type="l") # the blue one is PC2
dev.off()
但是，此示例绝不是您数据的最佳 PCA，因为 PC1 和 PC2 存在严重的季节性和年度变化（当然，夏季下雨更多，并查看 PC 的弱尾部)。
您可以通过对数据进行去季节化处理或通过回归去除年度趋势来改进 PCA，正如文献中所建议的那样。但这已经超出了我们的主题。