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
import os

#1、读取多光谱分辨率数据
def read_mul(input_path):
    dataset = gdal.Open(input_path)
    if dataset == None:
        print("输入文件无法打开")
    low_width = dataset.RasterXSize                                             #栅格矩阵的列数
    low_height = dataset.RasterYSize                                            #栅格矩阵的行数
    low_bands = dataset.RasterCount                                             #波段数
    lowim_data = dataset.ReadAsArray(0,0,low_width,low_height)                      #获取数据
    return low_width,low_height,low_bands,lowim_data,dataset

#2、读取全色分辨率数据
def read_pan(input_path):
    dataset_pan = gdal.Open(input_path)
    if dataset == None:
        print("输入文件无法打开")
    pan_width = dataset_pan.RasterXSize                                             #栅格矩阵的列数
    pan_height = dataset_pan.RasterYSize                                            #栅格矩阵的行数
    pan_bands = dataset_pan.RasterCount                                             #波段数
    pan_data = dataset_pan.ReadAsArray(0,0,pan_width ,pan_height)                      #获取数据
    return pan_width,pan_height,pan_bands,pan_data,dataset_pan

#3、计算均值,要求输入数据为numpy的矩阵格式，行表示样本数，列表示特征
def meanX(dataX):
    return np.mean(dataX, axis=0)  # axis=0表示按照列来求均值，如果输入list,则axis=1

#4计算方差,传入的是一个numpy的矩阵格式，行表示样本数，列表示特征
def variance(X):
    m, n = np.shape(X)
    mu = meanX(X)
    muAll = np.tile(mu, (m, 1))
    X1 = X - muAll
    variance = X1.T.dot(X1)/(m-1) 
#    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)) #m行n列的一个均值
    X1 = X - muAll#去中心化后的数据
    X2 = np.tile(np.diag(X.T * X), (m, 1))#np.diag创建一个对角矩阵
    XNorm = X1 / X2
    return XNorm
#参数：
#	- XMat：传入的是一个numpy的矩阵格式，行表示样本数，列表示特征    
#	- k：表示取前k个特征值对应的特征向量
#返回值：
#	- finalData：参数一指的是返回的低维矩阵，对应于输入参数二
#	- reconData：参数二对应的是移动坐标轴后的矩阵
#"""
def pca(XMat,k,height,width,prototype):
    m, n = np.shape(XMat)     #行，列表示波段数
    mean_vec = np.mean(XMat, axis=0)#对每一列求均值，输出结果为一行
    data_adjust=XMat - mean_vec
    cov_mat = np.cov(data_adjust.T)# 相当于cov_mat = (XMat - mean_vec).T.dot((XMat - mean_vec)) / (m-1)#计算协方差矩阵
    eig_vals, eig_vecs = np.linalg.eig(cov_mat)#计算协方差矩阵得特征值和特征向量
    index = np.argsort(-eig_vals)  # 按照featValue进行从大到小排序
    eig_vals=eig_vals[index]
    #确定k的主成分的个数
    sum=0
    max=0
    for i in range(len(index)):
        sum=sum+eig_vals[index[i]]
        if max/sum<0.99:
            max=max+eig_vals[index[i]]
            k=k+1
    finalData = []
    print ("k = ",k)
    if k > n:
        print ("k must lower than feature number")
    else:
        # 注意特征向量是列向量，而numpy的二维矩阵(数组)a[m][n]中，a[1]表示第1行值
        selectVec = np.matrix(eig_vecs[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, finalData[:,j])           #将reconData从二维（3列）数组转化从一维数组（1维）
        print (np.shape(array))
        array1 = array.reshape(k, height, width)             #将 一维数组  转成  3维矩阵（k,高（行），宽（列））
        out = ga.SaveArray(array1, os.path.join(path, "after.img"), format="GTiff", prototype=prototype)
        print ("主成分分析结果：主成分数、图像大小",np.shape(array1))
    #return finalData, reconData
    return  finalData

if __name__=="__main__":
    path="C:\\Users\\User\\Desktop\\image"
    mul_path="C:\\Users\\User\\Desktop\\image\\1.tif"
    high_path="C:\\Users\\User\\Desktop\\image\\2.tif"
    output_path="C:\\Users\\User\\Desktop\\image\\res.tif"
    low_width,low_height,low_bands,lowim_data,dataset=read_mul(mul_path)
    pan_width,pan_height,pan_bands,pan_data,dataset_pan=read_pan(high_path)
    if low_bands!=4:
        print("图像波段数有误，目前只能处理4波段图像！")
    if (pan_width != low_width and pan_height != low_height):
        print("两幅图像大小不一致")
    
#    mul_data = lowim_data.ReadAsArray(0, 0, low_width, low_height)  # 获取数据 将数据写成数组，对应栅格矩阵
    data=[]
    for j in range(low_bands):
           picture = lowim_data[j,:,:].flatten()  # 变成一维数组
           data.append(picture)            #将各个波段合并在一个数组里，一行为一个波段（n个波段n行）
           #(4,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个波段）
    pca(XMat, 1,low_height,low_width,dataset)