Author:v0_persistent
地理信息科学专业
把成果分享化为持续学习的动力!
前言
计算机遥感图像分类是通过模式识别理论,利用计算机将遥感图像分为若干地物类别的方法。提取信息类型包括分类、变化监测、物理量的提取、指标提取、特定地物和状态提取等,在多领域有广泛的应用。
本文将使用python进行遥感图像分类。
图像分类原理:依据和方法
统计模式识别的关键是提取待识别模式的一组统计特征值,然后按一定准则作出决策,对图像予以识别。
遥感图像分类的主要依据是地物的光谱特征,即地物电磁波辐射的多波段测量值,测量值可用作遥感图像分类的特征变量。python可以作为一种决策工具,将分类的准则和步骤写成代码令其执行,实现遥感图像分类。
本文将采用监督分类的最小距离分类法实现图像分类。
python实现图像分类
数据准备
已经过大气校正的遥感数据(本文以landset8数据为例)、ROI文件
获取训练区样本数据
(若已有ROI文件跳过此步)
在envi中,对已经进行过大气校正的遥感图像进行感兴趣区的选择,将ROI导出为dat文件
不同类别地物对应不同的value值,此处用于后面的python分类。
进行最小距离分类
工具包的准备和数据引入
用到rasterio、matplotlib、numpy工具包,引入原遥感图像的7个波段数据和ROI图像7波段数据。
import rasterio as rio
import matplotlib.pyplot as plt
import numpy as np
import matplotlib
matplotlib.use('TkAgg')
# 原始图像7波段数据
data = rio.open('D:\\envi_practice\\遥感制图\\flaash_result\\image_rc_ac.dat')
band1 = data.read(1)
band2 = data.read(2)
band3 = data.read(3)
band4 = data.read(4)
band5 = data.read(5)
band6 = data.read(6)
band7 = data.read(7)
# 感兴趣区7波段数据
train_data = rio.open('D:\\envi_practice\\遥感制图\\tx_fenlei\\train_roi\\train_roi.dat')
# train_roi = train_data.read(1)
# plt.figure()
# plt.imshow(train_roi, cmap='summer', vmin=0, vmax=20)
# plt.show()
提取相应地物的栅格
ROI文件中不同的地物已经对应了不同的value值,此处借助ROI来提取原遥感图像的栅格。
# 水体
b1_water_mask = train_data == 3
# plt.imshow(b1_water_mask, cmap='grey')
b1_water = band1[train_roi == 3]
b2_water = band2[train_roi == 3]
b3_water = band3[train_roi == 3]
b4_water = band4[train_roi == 3]
b5_water = band5[train_roi == 3]
b6_water = band6[train_roi == 3]
b7_water = band7[train_roi == 3]
# 垂直堆叠 vertical stack
# vstack接受一个数组序列
train_water = np.vstack((b1_water, b2_water, b3_water, b4_water, b5_water, b6_water, b7_water))
train_water_center = train_water.mean(axis=1)
# 以此类推
# 建筑
b1_building_mask = train_data == 1
# plt.imshow(b1_building_mask, cmap='grey')
b1_building = band1[train_roi == 1]
b2_building = band2[train_roi == 1]
b3_building = band3[train_roi == 1]
b4_building = band4[train_roi == 1]
b5_building = band5[train_roi == 1]
b6_building = band6[train_roi == 1]
b7_building = band7[train_roi == 1]
train_building = np.vstack((b1_building, b2_building, b3_building, b4_building, b5_building, b6_building, b7_building))
train_building_center = train_building.mean(axis=1)
# 森林
b1_forest_mask = train_data == 2
# plt.imshow(b1_forest_mask, cmap='grey')
b1_forest = band2[train_roi == 2]
b2_forest = band2[train_roi == 2]
b3_forest = band3[train_roi == 2]
b4_forest = band4[train_roi == 2]
b5_forest = band5[train_roi == 2]
b6_forest = band6[train_roi == 2]
b7_forest = band7[train_roi == 2]
train_forest = np.vstack((b1_forest, b2_forest, b3_forest, b4_forest, b5_forest, b6_forest, b7_forest))
train_forest_center = train_forest.mean(axis=1)
# 渔场
b1_fishing_mask = train_data == 4
# plt.imshow(b1_fishing_mask, cmap='grey')
b1_fishing = band1[train_roi == 4]
b2_fishing = band2[train_roi == 4]
b3_fishing = band3[train_roi == 4]
b4_fishing = band4[train_roi == 4]
b5_fishing = band5[train_roi == 4]
b6_fishing = band6[train_roi == 4]
b7_fishing = band7[train_roi == 4]
train_fishing = np.vstack((b1_fishing, b2_fishing, b3_fishing, b4_fishing, b5_fishing, b6_fishing, b7_fishing))
train_fishing_center = train_fishing.mean(axis=1)
# 农地
b1_farmland_mask = train_data == 5
# plt.imshow(b1_farmland_mask, cmap='grey')
b1_farmland = band1[train_roi == 5]
b2_farmland = band2[train_roi == 5]
b3_farmland = band3[train_roi == 5]
b4_farmland = band4[train_roi == 5]
b5_farmland = band5[train_roi == 5]
b6_farmland = band6[train_roi == 5]
b7_farmland = band7[train_roi == 5]
train_farmland = np.vstack((b1_farmland, b2_farmland, b3_farmland, b4_farmland, b5_farmland, b6_farmland, b7_farmland))
train_farmland_center = train_farmland.mean(axis=1)
# 展示各类地物的光谱曲线
plt.plot(train_water_center, c='blue', label='water')
plt.plot(train_building_center, c='gray', label='building')
plt.plot(train_forest_center, c='green', label='forest')
plt.plot(train_fishing_center, label='fishing')
plt.plot(train_farmland_center, c='coral', label='farmland')
plt.title('the center of all class')
plt.legend()
plt.show()
各类地物光谱曲线如图所示
逐个像元分类生成分类图
row_sum = band1.shape[0]
col_sum = band1.shape[1]
# 创建一个与原遥感图像大小相同的空白数组,用于生成分类图像
image_class = np.zeros((row_sum, col_sum))
# 对逐个像元进行分类
for row in range(band1.shape[0]):
for col in range(band1.shape[1]):
# 取出每个像元的光谱
r = np.array([band1[row, col], band2[row, col], band3[row, col], band4[row, col],
band5[row, col], band6[row, col], band7[row, col]])
# 欧式距离求和
dist_water = np.sum((r - train_water_center) ** 2)
dist_building = np.sum((r - train_building_center) ** 2)
dist_forest = np.sum((r - train_forest_center) ** 2)
dist_fishing = np.sum((r - train_fishing_center) ** 2)
dist_farmland = np.sum((r - train_farmland_center) ** 2)
# 距离数组
distance = np.array([dist_water, dist_building, dist_forest, dist_fishing, dist_farmland])
# 找出数组中最小值出现的位置
min_distance = distance.argmin()
# 为像元分类
image_class[row, col] = min_distance + 1
plt.figure()
plt.imshow(image_class, cmap='jet')
plt.show()
最终成果如图
点赞关注支持 -v- ❤️
若有错误,虚心请教