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

  • g2_586989
    了解作者
  • 2.3KB
    文件大小
  • 文件格式
  • 0
    收藏次数
  • VIP专享
    资源类型
  • 0
    下载次数
  • 2022-04-28 05:16
    上传日期
PCA基本步骤: 对数据进行归一化处理(代码中并非这么做的,而是直接减去均值) 计算归一化后的数据集的协方差矩阵 计算协方差矩阵的特征值和特征向量 保留最重要的k个特征(通常k要小于n),也可以自己制定,也可以选择一个阈值,然后通过前k个特征值之和减去后面n-k个特征值之和大于这个阈值,则选择这个k 找出k个特征值对应的特征向量 将m * n的数据集乘以k个n维的特征向量的特征向量(n * k),得到最后降维的数据。
PCA.rar
  • PCA.py
    5.6KB
内容介绍
# -*- 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)
评论
    相关推荐
    • 高光谱遥感图像数据集
      常用的高光谱遥感图像数据集,包含Indians Pines、Botswana、KSC、PaviaU、Salinas及它们的ground truth矩阵。文件后缀为.mat
    • KITTI数据集点云地图构建资料.rar
      本文描述了如何通过KITTI数据集,读取激光雷达点云数据,并通过ground truth进行点云建图的过程。其代码的主要功能包括:1)点云文件的格式转换2)点云转换矩阵计算3)点云地图构建
    • pmf:MovieLens 100K上的概率矩阵分解
      在此项目中,我们使用MovieLens 100K数据集。 该数据集包含来自943位用户的1,682部电影的100,000个评分。 在此项目中,RMSE(均方根误差)用作度量。 我测试了2种不同的数据分割:密集和稀疏。 数据是随机拆分的,...
    • circuitdata:数据集
      连接组学数据集 此存档的目的是提供各种“计算系统”数据集的规范、清理、正确规范化的版本。 所有文件均以 SQLite 格式提供。 这是因为在处理这些数据集时很难同时获得干净的数据和元数据——亲爱的科学家们,请...
    • honst:根据您的规则修复您的数据集
      在编辑连接的矩阵数据时,需要保持数据完整性。 保持本地状态正确。 演示版 您可以在CodeSandbox上尽情玩乐: 概述 const data = [ { username : "johndoe" , name : "John" , surname : "Doe" , age : "22" } , ...
    • LFW-a 数据集mat文件.zip
      LFW-a 数据集子mat文件,已经做cropped处理切割人脸部分,并且resize成了30x15的大小。只包含样本个数大于等于20的类别。fea为样本矩阵,每行一个样本,gnd为类别标签。
    • 数据结构代码合集
      内含二叉树的遍历、哈夫曼编码、后缀表达式的实现、图的创建与遍历、图的应用、稀疏矩阵、线性表合并、一元多项式的乘法、加法、减法。
    • matrixTests:R包,用于在矩阵数据框的行列上计算多个假设检验
      使用“物种”作为组,对虹膜数据集的每一列进行Bartlett检验: col_bartlett( iris [, - 5 ], iris $ Species ) obs.tot obs.groups var.pooled df statistic pvalue Sepal.Length 150 3 0.26500816 2 16.005702 0...
    • ruse.js:为绘制大型数据集而构建的轻量级 WebGL 库
      ruse.js是一个轻量级的 WebGL 库,用于绘制大型数据集。 它能够一次绘制数百万个点,并在图之间平滑过渡。 依赖 ruse.js仅依赖于进行矩阵转换。 用法 使用常用的脚本标签将gl-matrix和ruse.js包含到项目中。 ...
    • 矩阵分解数据集
      矩阵分解所需数据集,里面数据集有6个,ratings.csv,