【Python&GIS】基于Python批量合并矢量数据

​

　　老样子最近有项目需要将N个矢量文件合并成一个，总不能用ArcGIS一个个导入吧。所以我就想着用Python编个程序实现批量合并矢量。我之前也发了一些关于Python操作矢量数据的文章：【Python&GIS】Python处理矢量数据的基本操作（查询、修改、删除、新建）,如果大家感兴趣可以去我的主页看看，给我点个关注！

一、导入库

import os
from osgeo import ogr

二、合并shp

        整体的思路就是创建一个空的shp资源，然后遍历文件夹中所有的shp，然后针对每一个shp再遍历它的要素，将每个要素写入创建的新shp中。需要注意的是最后需要释放内存，不然数据不会写入shp。

def Merge_shp(path):
    print("-----------------合并shp-----------------")
    path_lists = os.listdir(path)
    src_proj = None
    for c in path_lists:
        if c.endswith(".shp"):
            print(c)
            ds = ogr.Open(path+c)
            layer = ds.GetLayer()
            # 打开需要转换的矢量数据，获取图层
            src_proj = layer.GetSpatialRef()
            break
    # 获取其源坐标信息
    output_file = path+"Merge.shp"
    driver = ogr.GetDriverByName('ESRI Shapefile')
    output_ds = driver.CreateDataSource(output_file)
    output_layer = output_ds.CreateLayer("Shp", srs=src_proj, geom_type=ogr.wkbMultiPolygon)
    new_field = ogr.FieldDefn('value', ogr.OFTReal)  # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value
    output_layer.CreateField(new_field)
    for i in range(0, len(path_lists)):
        if path_lists[i].endswith(".shp"):
            print("正在合并%s......" % path_lists[i])
            input_ds = ogr.Open(path + path_lists[i])
            input_layer = input_ds.GetLayer()
            for feature in input_layer:
                output_layer.CreateFeature(feature.Clone())
            input_ds = None
    output_ds = None

三、获取要素面积

        我这里给shp添加了一个新的字段，用来计算面积。大家视情况而定，这个可以没有。

def Get_polygon_area(path_shp):
    """
    :param path_shp: 输入矢量文件
    :return:
    """
    print("---------------获取矢量面积---------------")
    driver = ogr.GetDriverByName("ESRI Shapefile")  # 创建数据驱动
    ds = driver.Open(path_shp, 1)  # 创建数据资源
    layer = ds.GetLayer()
    new_field = ogr.FieldDefn("Area", ogr.OFTReal)  # 创建新的字段
    # new_field.SetWidth(32)
    # new_field.SetPrecision(16)
    layer.CreateField(new_field)
    for feature in layer:
        geom = feature.GetGeometryRef()
        geom2 = geom.Clone()
        area = geom2.GetArea()  # 默认为平方米
        # area = area / 1000000 # 转化为平方公里
        feature.SetField("Area", area)
        # feature.GetField('Area')
        layer.SetFeature(feature)
    ds = None

四、完整代码

# -*- coding: utf-8 -*-
"""
@Time ： 2023/10/20 11:56
@Auth ： RS迷途小书童
@File ：Vector Data Batch Merge.py
@IDE ：PyCharm
@Purpose：矢量数据批量合并成一个文件并计算面积
"""
import os
from osgeo import ogr


def Merge_shp(path):
    print("-----------------合并shp-----------------")
    path_lists = os.listdir(path)
    src_proj = None
    for c in path_lists:
        if c.endswith(".shp"):
            print(c)
            ds = ogr.Open(path+c)
            layer = ds.GetLayer()
            # 打开需要转换的矢量数据，获取图层
            src_proj = layer.GetSpatialRef()
            break
    # 获取其源坐标信息
    output_file = path+"Merge.shp"
    driver = ogr.GetDriverByName('ESRI Shapefile')
    output_ds = driver.CreateDataSource(output_file)
    output_layer = output_ds.CreateLayer("Shp", srs=src_proj, geom_type=ogr.wkbMultiPolygon)
    new_field = ogr.FieldDefn('value', ogr.OFTReal)  # 给目标shp文件添加一个字段，用来存储原始栅格的pixel value
    output_layer.CreateField(new_field)
    for i in range(0, len(path_lists)):
        if path_lists[i].endswith(".shp"):
            print("正在合并%s......" % path_lists[i])
            input_ds = ogr.Open(path + path_lists[i])
            input_layer = input_ds.GetLayer()
            for feature in input_layer:
                output_layer.CreateFeature(feature.Clone())
            input_ds = None
    output_ds = None


def Get_polygon_area(path_shp):
    """
    :param path_shp: 输入矢量文件
    :return:
    """
    print("---------------获取矢量面积---------------")
    driver = ogr.GetDriverByName("ESRI Shapefile")  # 创建数据驱动
    ds = driver.Open(path_shp, 1)  # 创建数据资源
    layer = ds.GetLayer()
    new_field = ogr.FieldDefn("Area", ogr.OFTReal)  # 创建新的字段
    # new_field.SetWidth(32)
    # new_field.SetPrecision(16)
    layer.CreateField(new_field)
    for feature in layer:
        geom = feature.GetGeometryRef()
        geom2 = geom.Clone()
        area = geom2.GetArea()  # 默认为平方米
        # area = area / 1000000 # 转化为平方公里
        feature.SetField("Area", area)
        # feature.GetField('Area')
        layer.SetFeature(feature)
    ds = None


if __name__ == "__main__":
    path1 = r"G:/2/"
    Merge_shp(path1)
    Get_polygon_area(path1+"Merge.shp")

        大家在使用时只需要将你要合并的shp全部放入一个文件夹，然后再将main中文件夹路径改成自己的就行了，程序最后会在该目录下生成Merge.shp文件，这个就是合并之后的结果。