一文详解高光谱数据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) # 32
img_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))
参考:
https://www.spectralpython.net/index.html
https://zhuanlan.zhihu.com/p/625036862