基于GDAL的Python实现遥感影像PCA的代码

# -*- coding: utf-8 -*- """ Created -0on Sun Feb 28 10:04:26 2016 PCA source code @author: MM PCA基本步骤： 对数据进行归一化处理（代码中并非这么做的，而是直接减去均值） 计算归一化后的数据集的协方差矩阵 计算协方差矩阵的特征值和特征向量 保留最重要的k个特征（通常k要小于n），也可以自己制定，也可以选择一个阈值，然后通过前k个特征值之和减去后面n-k个特征值之和大于这个阈值，则选择这个k 找出k个特征值对应的特征向量 将m * n的数据集乘以k个n维的特征向量的特征向量（n * k）,得到最后降维的数据。 """ import numpy as np import gdal import os from osgeo import gdal_array as ga import matplotlib.pyplot as plt # 计算均值,要求输入数据为numpy的矩阵格式，行表示样本数，列表示特征 def meanX(dataX): return np.mean(dataX, axis=0) # axis=0表示按照列来求均值，如果输入list,则axis=1 # 计算方差,传入的是一个numpy的矩阵格式，行表示样本数，列表示特征 def variance(X): m, n = np.shape(X) mu = meanX(X) muAll = np.tile(mu, (m, 1)) X1 = X - muAll variance = 1. / m * np.diag(X1.T * X1) return variance # 标准化,传入的是一个numpy的矩阵格式，行表示样本数，列表示特征 def normalize(X): m, n = np.shape(X) mu = meanX(X) muAll = np.tile(mu, (m, 1)) X1 = X - muAll X2 = np.tile(np.diag(X.T * X), (m, 1)) XNorm = X1 / X2 return XNorm """ 参数： - XMat：传入的是一个numpy的矩阵格式，行表示样本数，列表示特征 - k：表示取前k个特征值对应的特征向量 返回值： - finalData：参数一指的是返回的低维矩阵，对应于输入参数二 - reconData：参数二对应的是移动坐标轴后的矩阵 """ def pca(XMat, k): average = meanX(XMat) ##计算均值,为进行归一化做准备 m, n = np.shape(XMat) #行，列（7） print (m,n) data_adjust = [] avgs = np.tile(average, (m, 1)) #重复某个数组。比如tile(A,n)，功能是将数组A重复n次，构成一个新的数组 data_adjust = XMat - avgs # 对数据进行归一化处理（代码中并非这么做的，而是直接减去均值） covX = np.cov(data_adjust.T) # 计算协方差矩阵 featValue, featVec = np.linalg.eig(covX) # 求解协方差矩阵的特征值和特征向量 print (np.shape(featVec)) index = np.argsort(-featValue) # 按照featValue进行从大到小排序 #确定k的主成分的个数 sum=0 max=0 for i in range(len(index)): sum=sum+featValue[index[i]] if max/sum<0.99: max=max+featValue[index[i]] k=k+1 finalData = [] print ("k = ",k) if k > n: print ("k must lower than feature number") return else: # 注意特征向量是列向量，而numpy的二维矩阵(数组)a[m][n]中，a[1]表示第1行值 selectVec = np.matrix(featVec[index[:k]]) # 找出k个特征值对应的特征向量 print ("特征向量维数") print (np.shape(selectVec)) #finalData = data_adjust * selectVec.T reconData = (data_adjust * selectVec.T) + average[:,6] #将m * n的数据集乘以k个n维的特征向量的特征向量（m * k）,得到最后降维的数据。 print( np.shape(reconData)) array = [] array = np.array(array) # 列表转数组 for j in range(k): array = np.append(array, reconData[:,j]) #将reconData从二维（3列）数组转化从一维数组（1维） print (np.shape(array)) array1 = array.reshape(k, im_height, im_width) #将 一维数组 转成 3维矩阵（k,高（行），宽（列）） out = ga.SaveArray(array1, os.path.join(path, "after.img"), format="GTiff", prototype=img) #array2=array.reshape(im_height,im_width,k) #print (np.shape(array2)) #out1=ga.SaveArray(array2,os.path.join(path,"after22.img"),format="GTiff",prototype=img) #plt.imshow(array2, ), plt.show() print ("主成分分析结果：主成分数、图像大小",np.shape(array1)) #return finalData, reconData return reconData #- XMat：传入的是一个numpy的矩阵格式，行表示样本数，列表示特征 # 根据数据集 path = "G:/Experiment/kaixueshixi/任务数据/1 主成分分析" input = os.listdir(path) data = [] for i in input: if i == 'before.img': img = gdal.Open(os.path.join(path, i)) im_width = img.RasterXSize # 栅格矩阵的列数 im_height = img.RasterYSize # 栅格矩阵的行数 im_bands = img.RasterCount # 波段数 im_data = img.ReadAsArray(0, 0, im_width, im_height) # 获取数据 将数据写成数组，对应栅格矩阵 #im_data:(波段数，行，列) for j in range(im_bands): picture = im_data[j,:,:].flatten() # 变成一维数组 data.append(picture) #将各个波段合并在一个数组里，一行为一个波段（n个波段n行） #(7,rows*cols) print ("data维数", np.shape(data)) data = np.mat(data) #数组转换为矩阵 data1 = data.T #进行转置 #或者 data1=np.transpose(data) print ("data1维数", np.shape(data1)) XMat = data1 # - XMat：是一个numpy的矩阵格式，行表示样本数，列表示特征（每一列表示一个波段的信息，共n个波段） k = 0 # k：表示取前k个特征值对应的特征向量 pca(XMat,k)