【Python&RS】基于矢量点读取遥感影像波段值&制作训练样本

RS迷途小书童的博客 / 2024-08-05 / 原文

在进行遥感定量反演或数据分析时，往往我们都具有矢量的真值，可能是点文件也可能是面文件，最重要的还是通过这个矢量获取影像中该区域的值，这样方便做波段分析以及后续的反演等流程。今天给大家分享一下如何通过点文件获取影像的波段值。

原创作者：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)
        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 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 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