python 根据矢量点在栅格中采样

3D_DLW

已于 2023-08-12 16:37:30 修改

阅读量177

点赞数

分类专栏：遥感和地理信息文章标签： python 计算机视觉 opencv

于 2022-05-21 22:15:14 首次发布

本文链接：https://blog.csdn.net/DLW__/article/details/124903922

版权

工具同时被 2 个专栏收录

17 篇文章 0 订阅

订阅专栏

遥感和地理信息

14 篇文章 4 订阅

订阅专栏

该博客介绍了如何使用Python将含有地理坐标的.shp文件中的点与.tif栅格图像的值进行对应。首先获取.tif文件的坐标范围和像元大小，然后读取.shp文件中的点坐标和标签，接着将地理坐标转换为像素坐标，最后从栅格图像中提取对应位置的值。整个过程涉及地理信息系统和图像处理技术。

摘要由CSDN通过智能技术生成

前言

我有一个.shp文件，里面存了很多矢量点，每个矢量点有一个标签值。
我还有一个.tif文件，里面有一张栅格图片。
如何使用 python 将 .tif 中栅格的值与 .shp 中点的值一一对应？

方法

1 获取图像坐标

由于我的.shp文件和.tif文件都是以地理坐标定义的，在对其进行处理前需要获取一些参数。

使用ArcMap或其他方法查看tiff的坐标范围和像元大小
在这里插入图片描述

得：

windows = (529000, 2500000, 530000, 2501000)  # 左 下 右 上
pixel_size = 0.25

2 读取点坐标和标签

def read_point_from_shp(path, label_name='label'):
    file = shapefile.Reader(path)
    # 确认shp类型
    assert file.shapeTypeName == 'POINT'  # shp类型
    # 读取点坐标和标签
    points = []
    labels = []
    shape_records = file.shapeRecords()
    for shape_record in shape_records:
        points.append(shape_record.shape.points[0])
        labels.append(shape_record.record[label_name])
    points = np.array(points)
    labels = np.array(labels)
    file.close()
    return points, labels

# 读取点数据
points, labels = read_point_from_shp(shp_path)
# 读取栅格数据
img = imageio.imread(img_path)

3 地理坐标转像素坐标

def coo_geo2pix(points, boundary, pixel_size=0.25):
    """
    地理坐标转像素坐标
    :param points: 点集(地理坐标) np.array(n, 2)
    :param boundary: 栅格边界(地理坐标) [左 下 右 上]
    :param pixel_size: 像元大小
    :return: 点集(栅格坐标)
    """
    pos = np.zeros_like(points)
    pos[:, 0] = boundary[3] - points[:, 1]
    pos[:, 1] = points[:, 0] - boundary[0]
    pos = (pos / pixel_size).astype(np.int32)
    return pos

# 定义边界信息
windows = (529910, 2538625, 530517.5, 2539218.25)  # 左 下 右 上
# 地理坐标转像素坐标，提供像元大小
pos = coo_geo2pix(points, windows, pixel_size=0.25)
# 越界检查
assert np.all(pos.min(0) > 0) and np.all(pos.max(0) < img.shape)

4 提取栅格中的值

samples = []
for p in pos:
    samples.append(img[p[0], p[1]])
samples = np.array(samples)

全部代码

import imageio
import numpy as np
import shapefile


def coo_geo2pix(points, boundary, pixel_size=0.25):
    """
    地理坐标转像素坐标
    :param points: 点集(地理坐标) np.array(n, 2)
    :param boundary: 栅格边界(地理坐标) [左 下 右 上]
    :param pixel_size: 像元大小
    :return: 点集(栅格坐标)
    """
    pos = np.zeros_like(points)
    pos[:, 0] = boundary[3] - points[:, 1]
    pos[:, 1] = points[:, 0] - boundary[0]
    pos = (pos / pixel_size).astype(np.int32)
    return pos


def read_point_from_shp(path, label_name='label'):
    file = shapefile.Reader(path)
    # 读取属性表头
    fields = file.fields  # 字段[name, type, size, decimal]
    assert file.shapeTypeName == 'POINT'  # shp类型
    # 读取点坐标和标签
    points = []
    labels = []
    shape_records = file.shapeRecords()
    for shape_record in shape_records:
        points.append(shape_record.shape.points[0])
        labels.append(shape_record.record[label_name])
    points = np.array(points)
    labels = np.array(labels)
    file.close()
    return points, labels


def sample_by_shp(img_path, shp_path, windows, pixel_size):
    # 读取点数据
    points, labels = read_point_from_shp(shp_path)

    # 读取栅格数据
    img = imageio.imread(img_path)

    # 地理坐标转像素坐标
    pos = coo_geo2pix(points, windows, pixel_size=pixel_size)
    assert np.all(pos.min(0) > 0) and np.all(pos.max(0) < img.shape)
    
    # 提取栅格值
    samples = []
    for p in pos:
        samples.append(img[p[0], p[1]])
    samples = np.array(samples)

    return labels, samples


if __name__ == '__main__':
    img_path = r'iamge.png'
    shp_path = r'label.shp'
    windows = (529000, 2500000, 530000, 2501000)  # 左 下 右 上
    pixel_size = 0.25
    labels, samples = sample_by_shp(img_path, shp_path, windows, pixel_size)