第一步是导入所需的库并加载数据。为了便于使访问和处理，波段图像被堆叠在大小为 850 x 1100 x 7（高 x 宽 x 波段数）的 3d numpy 阵列中。下面显示的彩色图像是红色、绿色和蓝色 (RGB) 波段图像的合成，再现了与我们所看到的相同的视图。
from IPython.display import Image, displayimport matplotlib.pyplot as pltimport pandas as pdimport seaborn as snsimport cv2import numpy as np # Read RGB image into an arrayimg = cv2.imread('band321.jpg')img_shape = img.shape[:2]print('image size = ',img_shape) # specify no of bands in the imagen_bands = 7 # 3 dimensional dummy array with zerosMB_img = np.zeros((img_shape[0],img_shape[1],n_bands)) # stacking up images into the arrayfor i in range(n_bands): MB_img[:,:,i] = cv2.imread('band'+str(i+1)+'.jpg', cv2.IMREAD_GRAYSCALE) # Let's take a look at sceneprint('\n\nDispalying colour image of the scene')plt.figure(figsize=(img_shape[0]/100,img_shape[1]/100))plt.imshow(img, vmin=0, vmax=255)plt.axis('off');
图像场景包括各种地表特征，例如水、建筑面积、森林和农田。
2. 数据探索
让我们看一下不同特征的单个波段图像的反射率，并尝试深入了解波段图像中的特征。
import matplotlib.pyplot as pltimport matplotlib.gridspec as grid fig,axes = plt.subplots(2,4,figsize=(50,23),sharex='all', sharey='all')fig.subplots_adjust(wspace=0.1, hspace=0.15)fig.suptitle('Intensities at Different Bandwidth in the visible and Infra-red spectrum', fontsize=30)axes = axes.ravel() for i in range(n_bands): axes[i].imshow(MB_img[:,:,i],cmap='gray', vmin=0, vmax=255) axes[i].set_title('band '+str(i+1),fontsize=25) axes[i].axis('off')fig.delaxes(axes[-1])
我们的变量是图像二维数组，需要转换为一维向量以便于矩阵计算。让我们创建一个大小为 935000 X 7（图像中的像素数 X 波段数）的变量矩阵，并将这些一维向量存储在其中。
# Convert 2d band array in 1-d to make them as feature vectors and StandardizationMB_matrix = np.zeros((MB_img[:,:,0].size,n_bands)) for i in range(n_bands): MB_array = MB_img[:,:,i].flatten() # covert 2d to 1d array MB_arrayStd = (MB_array - MB_array.mean())/MB_array.std() MB_matrix[:,i] = MB_arrayStdMB_matrix.shape;
在这一步中，实现数据压缩和降维。如果我们观察特征值，我们会发现这些值完全不同。这些值为我们提供了特征向量或方向的重要性顺序，即沿着特征向量的轴具有最大的特征值，是最重要的 PC 轴，依此类推。下一步是根据特征值从高到低对特征向量进行排序，按重要性顺序重新排列主成分。我们需要在有序特征向量的方向上转换数据，而有序特征向量又会产生主成分。
# Ordering Eigen values and vectorsorder = EigVal.argsort()[::-1]EigVal = EigVal[order]EigVec = EigVec[:,order] #Projecting data on Eigen vector directions resulting to Principal Components PC = np.matmul(MB_matrix,EigVec) #cross product
6. 主成分验证
依赖性检查
我们能够成功地生产主要成分。现在，让我们验证 PC 以检查它们是否能够减少冗余，并检查实现数据压缩的程度。我们将创建散点图来可视化原始波段中的成对关系，并将其与 PC 的成对关系进行比较以测试是否存在依赖性。
# Generate Paiplot for original data and transformed PCs Bandnames = ['Band 1','Band 2','Band 3','Band 4','Band 5','Band 6','Band 7']a = sns.pairplot(pd.DataFrame(MB_matrix, columns = Bandnames), diag_kind='kde',plot_kws={"s": 3})a.fig.suptitle("Pair plot of Band images") PCnames = ['PC 1','PC 2','PC 3','PC 4','PC 5','PC 6','PC 7']b = sns.pairplot(pd.DataFrame(PC, columns = PCnames), diag_kind='kde',plot_kws={"s": 3})b.fig.suptitle("Pair plot of PCs")
#Information Retained by Principal Components plt.figure(figsize=(8,6))plt.bar([1,2,3,4,5,6,7],EigVal/sum(EigVal)*100,align='center',width=0.4, tick_label = ['PC1','PC2','PC3','PC4','PC5','PC6','PC7'])plt.ylabel('Variance (%)')plt.title('Information retention');
上面绘制的以百分比表示的特征值条形图为我们提供了每个 PC 中保留的信息。请注意，最后一个 PC 的特征值很小且不那么重要，这就是降维的作用所在。如果我们选择保留 93% 信息的前三个相关组件，那么最终数据可以从 7 维减少到 3 维，而不会丢失太多信息。
7. 将 PC 转换回图像
是时候将一维 PC 重塑回原始图像形状并将PCs在 0 到 255 之间进行归一化，这与原始图像范围相同，以使图像可视化成为可能。
# Rearranging 1-d arrays to 2-d arrays of image size PC_2d = np.zeros((img_shape[0],img_shape[1],n_bands))for i in range(n_bands): PC_2d[:,:,i] = PC[:,i].reshape(-1,img_shape[1]) # normalizing between 0 to 255PC_2d_Norm = np.zeros((img_shape[0],img_shape[1],n_bands))for i in range(n_bands): PC_2d_Norm[:,:,i] = cv2.normalize(PC_2d[:,:,i], np.zeros(img_shape),0,255 ,cv2.NORM_MINMAX)
让我们直观地确定压缩量。
fig,axes = plt.subplots(2,4,figsize=(50,23),sharex='all', sharey='all')fig.subplots_adjust(wspace=0.1, hspace=0.15)fig.suptitle('Intensities of Principal Components ', fontsize=30) axes = axes.ravel()for i in range(n_bands): axes[i].imshow(PC_2d_Norm[:,:,i],cmap='gray', vmin=0, vmax=255) axes[i].set_title('PC '+str(i+1),fontsize=25) axes[i].axis('off')fig.delaxes(axes[-1])
# Comparsion of RGB and Image produced using first three bandsfig,axes = plt.subplots(1,2,figsize=(50,23), sharex='all', sharey='all') fig.subplots_adjust(wspace=0.1, hspace=0.15) axes[0].imshow(MB_img[:,:,0:3].astype(int))axes[0].axis('off'); axes[1].imshow(PC_2d_Norm[:,:,:3][:,:,[0,2,1]].astype(int))axes[1].axis('off');
