【Python&RS】基于矢量点读取遥感影像波段值&制作训练样本
在进行遥感定量反演或数据分析时，往往我们都具有矢量的真值，可能是点文件也可能是面文件，最重要的还是通过这个矢量获取影像中该区域的值，这样方便做波段分析以及后续的反演等流程。今天给大家分享一下如何通过点文件获取影像的波段值。
原创作者：RS迷途小书童
注意：栅格影像和矢量点文件都应具有相同坐标系！！！
1 获取label值
我这里分析时点矢量是具有多个字段的，这些字段都是标签值，或者可以说是测量的真值，如果你只有一种类型的真值可以自己修改一下。将所有的真值和坐标写入列表中并返回。
# -*- coding: utf-8 -*-
"""
@Time ： 2024/3/5 14:48
@Auth ： RS迷途小书童
@File ：Get_Raster_datas_from_Points.py
@IDE ：PyCharm
@Purpose：通过矢量点获取栅格数据的值（多波段）
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""
import os
import sys
import numpy as np
from datetime import datetime
from osgeo import gdal, osr, ogr
# 从点矢量中获取标签数据
def Get_label_point(train_image, train_shp, class_id):
"""
:param train_image: 输入训练的影像
:param train_shp: 输入训练的点矢量
:param class_id: 输入训练的点矢量对应类别的字段名
:return: x,y,类别值
"""
ogr.RegisterAll() # 注册所有的驱动
ds_image = gdal.Open(train_image)
ds_proj_image = osr.SpatialReference()
ds_proj_image.ImportFromWkt(ds_image.GetProjectionRef())
ds_shp = ogr.Open(train_shp)
if ds_shp is None:
print("打开文件【%s】失败！", train_shp)
return
layer = ds_shp.GetLayerByIndex(0) # 通过索引获取shapefile的第一个图层
if layer is None:
print("获取第%d个图层失败！\n")
return
ds_proj_shp = layer.GetSpatialRef()
transform = osr.CreateCoordinateTransformation(ds_proj_shp, ds_proj_image)
# # 通过索引获取shapefile的第一个图层
layer.ResetReading()
# 重置图层的读取位置到开始位置
feature = layer.GetNextFeature()
values = list()
while feature is not None:
value = feature.GetField(class_id) # 通过字段名获取该要素的类ID属性值
value1 = feature.GetField("1")
value2 = feature.GetField("2")
value3 = feature.GetField("3")
value4 = feature.GetField("4")
value5 = feature.GetField("5")
geometry = feature.GetGeometryRef()
x = geometry.GetX(0)
y = geometry.GetY(0)
point_transform = transform.TransformPoint(x, y, 0)
values.append([point_transform[0], point_transform[1], value, value1, value2, value3, value4, value5, x, y])
# 将坐标、字段加入到标签数组中
feature = layer.GetNextFeature() # 继续下一个要素的读取
ds = None
return values2 获取栅格值
通过上一步存储的坐标可以定位到栅格影像中的实际位置，通过波段的循环将所有波段值写入列表/数组或者表格中，我这里是写入表格中。
# -*- coding: utf-8 -*-
"""
@Time ： 2024/3/5 14:48
@Auth ： RS迷途小书童
@File ：Get_Raster_datas_from_Points.py
@IDE ：PyCharm
@Purpose：通过矢量点获取栅格数据的值（多波段）
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""
import os
import sys
import numpy as np
from datetime import datetime
from osgeo import gdal, osr, ogr
# 通过点矢量获取栅格各波段的值
def Get_train_data_point(train_image, train_shp, class_id, work_path):
"""
:param train_image: 输入训练的影像
:param train_shp: 输入训练的点矢量
:param class_id: 输入训练的点矢量对应类别的字段名
:param work_path: 工作空间
:return: 训练的特征值，训练的样本类别
"""
print("【数据准备阶段】")
print("[%s]获取样本点数据......" % datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
import openpyxl
wb = openpyxl.load_workbook(r'G:\RS迷途小书童/2.xlsx') # read_only=True 只读打开可以提高效率
ws_sheets = wb.sheetnames # 获取所有的sheet名称
ws = wb[ws_sheets[0]] # 激活指定的工作空间
os.chdir(work_path)
ds_image = gdal.Open(train_image)
ds_geo = ds_image.GetGeoTransform()
left_x, left_y, resolution_x, resolution_y = ds_geo[0], ds_geo[3], ds_geo[1], ds_geo[5]
point_values = Get_label_point(train_image, train_shp, class_id) # 读取点文件中的位置以及类别属性
ds_bands = ds_image.RasterCount # 获取影像的波段数
train_x, train_y = list(), list()
z = 2
for point_value in point_values:
# 遍历所有的样本点
row, col = int((point_value[0] - left_x) / resolution_x), int((point_value[1] - left_y) / resolution_y)
# 通过仿射变换参数，获取矢量点对应的栅格行列数
try:
arr_point = ds_image.ReadAsArray(row, col, 1, 1)
point_data = list()
for i in range(ds_bands): # 遍历该点所有的波段
point_data.append(int(arr_point[i])) # 添加每个波段的值
train_x.append(point_data)
ws.cell(z, 1).value = point_data[0]
ws.cell(z, 2).value = point_data[1]
ws.cell(z, 3).value = point_data[2]
"""ws.cell(z, 4).value = point_data[3]
ws.cell(z, 5).value = point_data[4]
ws.cell(z, 6).value = point_data[5]
ws.cell(z, 7).value = point_data[6]
ws.cell(z, 8).value = point_data[7]
ws.cell(z, 9).value = point_data[8]
ws.cell(z, 10).value = point_data[9]
ws.cell(z, 11).value = point_data[10]
ws.cell(z, 12).value = point_data[11]
ws.cell(z, 13).value = point_value[2]
ws.cell(z, 14).value = point_value[3]
ws.cell(z, 15).value = point_value[4]
ws.cell(z, 16).value = point_value[5]
ws.cell(z, 17).value = point_value[6]
ws.cell(z, 18).value = point_value[7]
ws.cell(z, 19).value = point_value[8]
ws.cell(z, 20).value = point_value[9]"""
# 将每个点的特征添加至数组，本质为二维数组
train_y.append(float(point_value[2]))
# print((point_value))
# 将每个点对应的类别加入数组
except Exception as e:
print(e)
z += 1
train_y_num = len(train_y)
if train_y_num != 0:
print("[%s]共获取%s个样本点数据......" % (datetime.now().strftime('%Y-%m-%d %H:%M:%S'), train_y_num))
else:
print("[%s]未获取到训练样本,请检查样本数据！" % datetime.now().strftime('%Y-%m-%d %H:%M:%S'))
sys.exit(0)
# print(train_x)
# print(train_y)
np.save("train.npy", train_x)
np.save("label.npy", train_y)
print("--------------------------------------------------------------------------------------")
wb.save(r"G:\RS迷途小书童/2.xlsx")
return train_x, train_y
后续在分享使用Python机器学习分类时，再和大家分享如何使用栅格、矢量面读取栅格值并制作成训练集和样本集。本质其实就是通过仿射变换矩阵和循环波段读取到该点的值，如果大家有什么问题也可以留言交流。
