线矢量分割面矢量

ytkz2024-01-152024-09-09

最近遇到一个问题：怎么用线矢量去切割面矢量。这个功能看似很简单，其实还是要一定的功夫才能完成它。所以也才有了这篇技术文章的记录。

目的：使用线段分割多边形，并将分割后的多边形保存为Shapefile文件

代码

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import geopandas as gpd
from shapely.geometry import Polygon, LineString
from shapely.ops import split

def split_polygon_with_lines(polygon, lines):
    # 使用每条线段分割多边形
    for line in lines:
        parts = [split(poly, line) for poly in polygon]
        # 将列表扁平化并从GeometryCollection中获取几何体
        polygon = [item for sublist in parts for item in sublist.geoms]
    return polygon

def save_polygons_to_shp(polygon, lines, filename):
    # 使用线段分割多边形
    splitted = split_polygon_with_lines([polygon], lines)
    # 将Shapely几何体列表转换为GeoDataFrame
    gdf = gpd.GeoDataFrame(geometry=splitted)
    # 设置GeoDataFrame的坐标参考系统为WGS 84
    gdf.crs = "EPSG:4326"
    # 将GeoDataFrame保存为Shapefile
    gdf.to_file(filename)

if __name__ == '__main__':
    # 创建一个多边形
    polygon = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)])
    # 创建一些线段对象
    line1 = LineString([(5, 0), (5, 10)])
    line2 = LineString([(0, 5), (9, 5), (9, 6), (0, 6)])
    # 将线段对象放入一个列表中
    lines = [line1, line2]
    # 使用线段分割多边形
    save_polygons_to_shp(polygon, lines, 'output.shp')

面矢量可视化如下：

线矢量可视化如下：

<<<<<<< HEAD

e8e02ff5542145568539ca02109bd75f955339c9

分割后结果，在arcgis打开如下：

代码解析

逐行解析以上代码：

import geopandas as gpd: 导入geopandas库，并将其简写为gpd。geopandas是一个用于处理地理数据的Python库。
from shapely.geometry import Polygon, LineString: 从shapely.geometry库中导入Polygon和LineString。shapely库用于操作和分析平面几何对象。
from shapely.ops import split: 从shapely.ops库中导入split函数，该函数用于分割几何对象。
def split_polygon_with_lines(polygon, lines):: 定义一个函数，该函数接收一个多边形和一组线段，然后使用每条线段分割多边形。
for line in lines:: 对输入的每条线段进行遍历。
parts = [split(poly, line) for poly in polygon]: 使用列表推导式，对每个多边形使用split函数进行分割，并将结果存储在parts列表中。
polygon = [item for sublist in parts for item in sublist.geoms]: 使用列表推导式从GeometryCollection中获取几何体，并将结果存储在polygon列表中。
return polygon: 返回分割后的多边形列表。
def save_polygons_to_shp(polygon, lines, filename):: 定义一个函数，该函数接收一个多边形、一组线段和一个文件名，然后将分割后的多边形保存为Shapefile。
splitted = split_polygon_with_lines([polygon], lines): 调用前面定义的函数，使用线段分割多边形，并将结果存储在splitted变量中。
gdf = gpd.GeoDataFrame(geometry=splitted): 将分割后的多边形列表转换为GeoDataFrame。
gdf.crs = "EPSG:4326": 设置GeoDataFrame的坐标参考系统为WGS 84。
gdf.to_file(filename): 将GeoDataFrame保存为Shapefile。
if __name__ == '__main__':: 如果当前脚本被直接运行，而不是作为模块导入，则执行以下代码。
polygon = Polygon([(0, 0), (10, 0), (10, 10), (0, 10)]): 创建一个多边形对象。
line1 = LineString([(5, 0), (5, 10)]) 和 line2 = LineString([(0, 5), (9, 5), (9, 6), (0, 6)]): 创建两条线段对象。
lines = [line1, line2]: 将两条线段对象放入一个列表中。
save_polygons_to_shp(polygon, lines, 'output.shp'): 调用前面定义的函数，使用线段分割多边形，并将结果保存为Shapefile。