在进行遥感定量反演或数据分析时，往往我们都具有矢量的真值，可能是点文件也可能是面文件，最重要的还是通过这个矢量获取影像中该区域的值，这样方便做波段分析以及后续的反演等流程。今天给大家分享一下如何通过点文件获取影像的波段值。
原创作者：RS迷途小书童
博客地址：https://blog.csdn.net/m0_56729804?type=blog
注意：栅格影像和矢量点文件都应具有相同坐标系！！！
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)returnlayer = ds_shp.GetLayerByIndex(0) # 通过索引获取shapefile的第一个图层if layer is None:print("获取第%d个图层失败！\n")returnds_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 = Nonereturn values
2 获取栅格值
通过上一步存储的坐标可以定位到栅格影像中的实际位置，通过波段的循环将所有波段值写入列表/数组或者表格中，我这里是写入表格中。
# -*- 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 openpyxlwb = 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 = 2for 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 += 1train_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机器学习分类时，再和大家分享如何使用栅格、矢量面读取栅格值并制作成训练集和样本集。本质其实就是通过仿射变换矩阵和循环波段读取到该点的值，如果大家有什么问题也可以留言交流。
