在R中基于栅格栈PCA绘制带相关性的论文同款主成分图

AI 小北

2026-5-29

Hey there! Let's work through how to create those standardized principal component plots with correlation annotations, just like Figure 1 in your target paper. I’ll build on your existing R code and break this down step by step:

1. Finalize PCA & Extract Key Outputs

First, let’s tweak your PCA code to match the paper’s 8 principal components, and ensure we’re capturing scores and loadings (the correlations between PCs and original monthly NDVI data):

library(psych)
library(raster)

# Replace your sample raster stack with your actual 36-month NDVI stack (1986-1988)
# For example: x <- stack(list.files(path = "your_ndvi_folder", pattern = ".tif", full.names = TRUE))

# Extract pixel values and transpose for PCA (rows = months, cols = pixels)
extract.value.from.raster.stack <- extract(x, 1:ncell(x))
pcaan <- t(extract.value.from.raster.stack)

# Run PCA for 8 components, keep scores, and standardize data (matches "standardized" in the paper)
pca3 <- principal(pcaan, nfactor = 8, rotate = "none", scores = TRUE, scale = TRUE)

2. Convert PCA Scores Back to Raster Layers

We need to turn the PCA scores (one per pixel per component) back into spatial raster layers, which will form the base of our plots:

# Extract the first 8 principal component scores
pc_scores <- pca3$scores[, 1:8]

# Create a raster stack for all 8 PCs, matching the original NDVI raster's extent/projection
pc_raster_stack <- stack()
for (i in 1:8) {
  pc_raster <- raster(x)  # Reuse the original raster's spatial properties
  values(pc_raster) <- pc_scores[, i]  # Assign scores to the raster
  pc_raster_stack <- addLayer(pc_raster_stack, pc_raster)
}

3. Calculate & Prepare Correlation Annotations

The paper’s plots likely include correlations between each PC and the original monthly NDVI data (these are the loadings from the PCA output). Let’s extract and format these for plotting:

# Extract the loading matrix (rows = months, cols = PCs)
loadings_matrix <- as.data.frame(unclass(pca3$loadings[, 1:8]))
rownames(loadings_matrix) <- c(paste0("1986-", 1:12), paste0("1987-", 1:12), paste0("1988-", 1:12))

# For each PC, get the top 3 months with the strongest correlations (positive or negative)
plot_titles <- lapply(1:8, function(pc_num) {
  top_cor <- head(sort(abs(loadings_matrix[, pc_num]), decreasing = TRUE), 3)
  title_text <- paste0("PC ", pc_num, "\nTop Correlations:\n",
                       paste(names(top_cor), round(top_cor, 2), sep = ": ", collapse = "\n"))
  return(title_text)
})

4. Plot the Standardized PC Images with Annotations

Use rasterVis for clean, grid-aligned plotting that includes our correlation annotations:

library(rasterVis)

# Plot 8 PCs in a 2x4 grid, with custom titles showing key correlations
levelplot(pc_raster_stack, 
          layout = c(4, 2),  # 4 columns, 2 rows
          main = "Standardized Principal Component Images (1986-1988 Monthly NDVI)",
          names.attr = plot_titles,  # Add our correlation titles to each subplot
          col.regions = colorRampPalette(c("darkblue", "white", "darkred"))(100))

This will generate a set of 8 spatial plots, each showing the spatial distribution of a standardized principal component, with annotations highlighting which months have the strongest correlation to that component—just like the figure in your paper!

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