一文详解高光谱数据python处理包spectral(SPy)

一、基本操作

读取高光谱数据文件

import spectral
# 读取ENVI格式的高光谱图像
# image的后缀可以是.raw、.spe、.lan等
# 代码里img对象，类似于rasterio库的dataset对象，可以用它来读取高光谱数据
img = spectral.envi.read_envi(file='my_data.hdr', image='my_data')

加载数据

# 通过load函数获取数据
data = img.load() # 多维数组

读取元数据

# 读取元数据metadata，可以获取坐标系统和变换系数
metadata=img.metadata
projection = metadata['pseudo projection info'] 
rpc =metadata['rpc info']

获取高光谱图像的子集

# 获取高光谱图像的前10个波段
subset = img.extract(0, 10)# 读取单个波段
data0 = img.read_band(0)

获取波长信息

wavelength_units= img.bands.band_unit  # 获取波长单位
wl=img.bands.centers # 获取中心波长
bandwidth=img.bands.bandwidths # 获取光谱分辨率，半高全宽

显示高光谱图像并保存

 from spectral import *img = open_image('92AV3C.lan')view = imshow(img, (29, 19, 9)) # 显示save_rgb('rgb.jpg', img, [29, 19, 9]) # 保存

在这里插入图片描述

查看view的具体信息：

# 可以看到展示的三个RGB波段的索引以及值域
print(view)
ImageView object:Display bands       :  (29, 19, 9)Interpolation       :  <default>RGB data limits     :R: [2054.0, 6317.0]G: [2775.0, 7307.0]B: [3560.0, 7928.0]

显示光谱信息

可以通过以下方式开启交互功能：

import matplotlib  
matplotlib.use('WXAgg')  # 注意：这行代码必须在导入pyplot之前执行

交互功能开启后，加载图像，点击图像上的像素点，即可出现相应的光谱曲线：

# imshow函数可以显示三波段的图像，并且可以显示指定位置的光谱信息view = imshow(img, (29, 19, 9)) # 显示

在这里插入图片描述

二、光谱算法

k-means 聚类

# k-means 聚类
(m, c) = kmeans(image=img, nclusters=20, max_iterations=30)

在这里插入图片描述

查看聚类中心：

import matplotlib.pyplot as plt
plt.figure()
for i in range(c.shape[0]):plt.plot(c[i])
plt.grid()

在这里插入图片描述

监督分类

执行监督分类需要使用训练数据训练分类器，该分类器将样本与特定训练类相关联。要将类标签分配给具有 M 行和 N 列的图像中的像素，必须提供一个 MxN 整数值数组，其元素是相应训练类的索引。真值数组中的值为 0 表示未标记的像素（该像素未与训练类关联）。

# 加载并显示示例图像的地面实况地图
gt = open_image('92AV3GT.GIS').read_band(0)
v = imshow(classes=gt)# 通过调用 create_training_classes 来创建 TrainingClassSet 对象
classes = create_training_classes(img, gt)# 定义训练类后，创建一个监督式分类器，使用训练类对其进行训练，然后对图像进行分类
gmlc = GaussianClassifier(classes) # 高斯最大似然法# 训练分类器后，使用它来对与训练集具有相同光谱波段的图像进行分类
# 对训练图像进行分类并显示生成的分类图。
clmap = gmlc.classify_image(img)
v = imshow(classes=clmap)

在这里插入图片描述
上面的分类图显示了整个图像的分类结果。要仅查看真值像素的结果，必须屏蔽掉与训练类无关的所有像素。

gtresults = clmap * (gt != 0)
v = imshow(classes=gtresults)

在这里插入图片描述

如果分类结果良好，预计上面的分类图看起来与原始真实地图非常相似。要仅查看errors，必须屏蔽 gtresults 中与真实图像不匹配的所有元素。

 gterrors = gtresults * (gtresults != gt)v = imshow(classes=gterrors)

在这里插入图片描述
上图中误差图中的五个连续区域对应于 GaussianClassifier 忽略的真值类，因为它们的样本太少。

除了 GaussianClassifier，还提供了其他分类方法。

降维

处理具有数百个波段的高光谱图像可能会造成计算负担，并且由于所谓的“维数诅咒”，分类精度可能会受到影响。为了缓解这些问题，通常需要降低数据的维度。

# principal_components计算图像数据的主成分
# 并返回 PrincipalComponents 中的平均值、协方差、特征值和特征向量
pc = principal_components(img)
v = imshow(pc.cov)

在这里插入图片描述
在协方差矩阵显示中，较白的值表示较强的正协方差，较暗的值表示较强的负协方差，灰色值表示接近零的协方差。

为了使用主成分降低维度，可以按降序对特征值进行排序，然后保留足够的特征值(相应的特征向量)来捕获总图像方差的所需部分。然后，将图像像素投影到剩余的特征向量上，以降低图像像素的维数。选择保留至少 99.9% 的总图像差异。

 pc_0999 = pc.reduce(fraction=0.999)# How many eigenvalues are left?len(pc_0999.eigenvalues) # 32img_pc = pc_0999.transform(img)v = imshow(img_pc[:,:,:3], stretch_all=True)

在这里插入图片描述
现在，对缩减的主成分使用高斯最大似然分类器(GMLC)来针对训练数据进行训练和分类：

classes = create_training_classes(img_pc, gt)
gmlc = GaussianClassifier(classes)
clmap = gmlc.classify_image(img_pc)
clmap_training = clmap * (gt != 0)
v = imshow(classes=clmap_training)

在这里插入图片描述
相关的errors：

training_errors = clmap_training * (clmap_training != gt)
v = imshow(classes=training_errors)

在这里插入图片描述

三、应用

案例：植被与土壤的分类和识别

要求1：提取地面、无人机高光谱数据的反射率

解决方案：ASD（Analytical Spectral Devices）是一种用于地面光谱测量的设备，可以测量多个波长范围内的反射率。提取ASD地面光谱数据通常可以采用以下步骤：

"""
1 读取ASD数据文件：ASD数据通常以文本文件的形式存储，可以使用Python中的Pandas库中的read_csv()函数来读取。
"""
import pandas as pd"""
2 数据清洗和预处理：数据通常需要进行清洗和预处理，例如去除无效数据点，对数据进行插值、平滑和校准等。
"""
data = pd.read_csv('ASD_data.csv')# 删除无效数据点
data = data.dropna()# 对数据进行插值
data = data.interpolate()# 对数据进行平滑处理
from scipy.signal import savgol_filter
data['refl_smooth'] = savgol_filter(data['refl'], 11, 2)"""
3 可以使用Python中的Spectral库来进行光谱分析，例如计算反射率
"""# 计算反射率
data['reflectance'] = data['refl'] / data['white']

要求2：采用机器学习如SVM、RF、CNN的方法进行识别

支持向量机（SVM）是一种常用的分类算法，可以用于高光谱数据的分类、识别和回归。

"""
以无人机获取的高光谱图像为例
""""""
1 读取高光谱数据
"""import spectral
img = spectral.open_image('data.hdr')"""
2 数据预处理
"""
from sklearn.decomposition import PCA
X = img.load()
X = X.reshape((X.shape[0]*X.shape[1], X.shape[2]))
pca = PCA(n_components=30)
X = pca.fit_transform(X)"""
3 分割数据集
"""
from sklearn.model_selection import train_test_split
y = img.read_band(145)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)from sklearn.model_selection import train_test_split
y = img.read_band(145)
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)"""
4 SVM分类
"""
使用Python中的sklearn库中的SVC()函数来进行SVM分类：
from sklearn.svm import SVC
clf = SVC(kernel='linear')
clf.fit(X_train, y_train)"""
5 评估分类器
"""
from sklearn.metrics import classification_report
y_pred = clf.predict(X_test)
print(classification_report(y_test, y_pred))