【Python&GIS】Python线矢量等距离取点/线等分取点&点创建矢量面

不多说，这是之前项目需求的代码，已经是去年的了一直没来的及发，今天抽出来一丢丢的空挡发一下。主要就是利用线矢量等距离生成点矢量，或者直接将线矢量等分生成点矢量，这个需求其实极限一下就是线转点了（将距离设置小一点）。顺便将点生成矩形面的代码也给出来，这里的矩形就直接中心点往外扩的固定距离，可以按自己的需求修改。

1.线等分取点代码

这里注意一下，我是提前知道的了线段的长度，所以可以直接用num_points = int(length/150) + 1 取到我希望的间距，这个代码是用来等分线段的，但是你计算好后也可以实现等距离！

# -*- coding: utf-8 -*-
"""
@Time ： 2024/11/22 17:10
@Auth ： RS迷途小书童
@File ：Vector Line Select Points.py
@IDE ：PyCharm
@Purpose：线矢量数据等距离取点/线等分取点，另加入点创建矢量面代码
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""
import numpy as np
import geopandas as gpd
from shapely import geometry
from shapely.geometry import Point, LineString, Polygon


def create_point(folder_path):  # 线矢量等距离取点
    """
    :param folder_path: 输入需要等分的矢量线文件
    :return: None
    """
    # shp_files = [f for f in os.listdir(folder_path) if f.endswith('.shp')]
    # os.mkdir(r"Dom_clip/%s" % files[:2])
    gdf = gpd.read_file(folder_path)  # 读取shp文件
    lines = gdf[gdf.geometry.type == 'LineString']  # 选择线要素
    line = LineString(lines.geometry.values[0].coords)  # 提取线的坐标并转换为LineString对象
    length = line.length  # 计算线的长度
    print("当前线矢量长度为：", length)
    num_points = int(length/150) + 1  # 设定你想要等分的点数
    distances = np.linspace(0, length, num_points)
    points = line.interpolate(distances)  # 计算等距点的坐标
    for point in points:  # 打印点的坐标
        print(point.x, point.y)
        number = 1  # 创建面的编号
        create_shp(point.x, point.y, number)  # 调用创建面的函数
        number += 1  # 创建面的编号

2.点创建面矩形代码

我这里直接是从上一步获得的中心点上下左右同时外扩获取的，你们可以按照需求来。

# -*- coding: utf-8 -*-
"""
@Time ： 2024/11/22 17:10
@Auth ： RS迷途小书童
@File ：Vector Line Select Points.py
@IDE ：PyCharm
@Purpose：线矢量数据等距离取点/线等分取点，另加入点创建矢量面代码
@Web：博客地址:https://blog.csdn.net/m0_56729804
"""
import numpy as np
import geopandas as gpd
from shapely import geometry
from shapely.geometry import Point, LineString, Polygon


def create_shp(x, y, number1):  # 创建面shp
    """
    :param x: 输入面矢量中心点x
    :param y: 输入面矢量中心点y
    :param number1: 保存新矢量的编号
    :return: None
    """
    distance = 78  # 矩形的大小
    shp = gpd.GeoSeries([geometry.Polygon([(x - distance, y + distance), (x + distance, y + distance),
                                           (x + distance, y - distance), (x - distance, y - distance)])],
                        crs='EPSG:32651')  # 指定坐标系为WGS84/UTM 51N
    # 左上、右上、左下、右下   32651右&上为正
    shp.to_file(r'彭俊喜/%s.shp' % number1, driver='ESRI Shapefile', encoding='utf-8')
    # 导出数据为shapefile文件