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

PCA基本步骤： 对数据进行归一化处理（代码中并非这么做的，而是直接减去均值） 计算归一化后的数据集的协方差矩阵 计算协方差矩阵的特征值和特征向量 保留最重要的k个特征（通常k要小于n），也可以自己制定，也可以选择一个阈值，然后通过前k个特征值之和减去后面n-k个特征值之和大于这个阈值，则选择这个k 找出k个特征值对应的特征向量 将m * n的数据集乘以k个n维的特征向量的特征向量（n * k）,得到最后降维的数据。

# -*- 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 （mn）
    data_adjust = []
    avgs = np.tile（average （m 1））   #重复某个数组。比如tile（An），功能是将数组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_heightim_widthk）
        #print （np.shape（array2））
        #out1=ga.SaveArray（array2os.path.join（path“after22.img“）format=“GTiff“prototyp

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       5746  2018-09-09 11:31  PCA.py

----------- ---------  ---------- -----  ----

                 5746                    1