如何在Python（Jupyter Notebook）中完成Landsat 5波段PCA变化检测

Got it, let's walk through completing the PCA for your Landsat bands in Python—this is straightforward with scikit-learn, we just need to adjust the data shape first since your stacked array is in spatial (rows, cols, bands) format. Here's a step-by-step breakdown tailored to your setup:

Step 1: Reshape your 3D data to 2D #

Scikit-learn's PCA expects a 2D array where each row is a sample (one pixel) and each column is a feature (one Landsat band). Your current data is (353, 588, 5)—we need to flatten the spatial dimensions into a single "samples" axis:

import numpy as np
from sklearn.decomposition import PCA

# Get the shape of your centered data
n_rows, n_cols, n_bands = data.shape

# Reshape to (total_pixels, n_bands)
data_2d = data.reshape(n_rows * n_cols, n_bands)

Step 2: Fit and apply PCA #

Next, initialize the PCA model. You can specify how many components to keep (e.g., 3 for most variance explained) or use a variance threshold (like 0.95 to retain 95% of the total variance). Let's stick to 3 components for this example:

# Initialize PCA with 3 components (adjust n_components as needed)
pca = PCA(n_components=3)

# Fit the model to your data and transform the pixels into PCA space
pca_transformed = pca.fit_transform(data_2d)

# Check how much variance each component captures (super useful for interpreting results!)
print("Variance explained per component:", pca.explained_variance_ratio_)
print("Total variance explained:", np.sum(pca.explained_variance_ratio_))

Quick side note: Scikit-learn's PCA automatically centers the data by subtracting the mean, so your earlier data = stacked - np.mean(stacked, axis=0) step is redundant if you use the default parameters. It doesn't break anything to keep it, but you could skip it next time if you want!

Step 3: Reshape back to 3D for spatial analysis #

To work with the PCA results in their original spatial context (for visualization or change detection), reshape the 2D PCA output back to the original (rows, cols, components) shape:

# Reshape to (n_rows, n_cols, n_components)
pca_spatial = pca_transformed.reshape(n_rows, n_cols, pca.n_components_)

Step 4: Using PCA for change detection #

Now that you have the PCA components for your 1984 data, you'd repeat the same steps for your second time period's Landsat bands. Once you have pca_spatial_1984 and pca_spatial_20XX, you can:

• Calculate the difference between corresponding components (e.g., change_component1 = pca_spatial_20XX[...,0] - pca_spatial_1984[...,0])
• Threshold these difference arrays to identify statistically significant changes
• Visualize the change maps to spot land cover shifts

For example, a simple thresholding step might look like:

# Example: Threshold the first component difference to detect large changes
significant_change = np.abs(change_component1) > 2 * np.std(change_component1)

That's it! This workflow will give you the PCA components you need for change detection, and it's fully compatible with Jupyter Notebook for interactive exploration.

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