中商外贸网站牛推网

        老样子最近有项目需要将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 = Nonefor 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 valueoutput_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 = Noneoutput_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 ogrdef Merge_shp(path):print("-----------------合并shp-----------------")path_lists = os.listdir(path)src_proj = Nonefor 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 valueoutput_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 = Noneoutput_ds = Nonedef 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 = Noneif __name__ == "__main__":path1 = r"G:/2/"Merge_shp(path1)Get_polygon_area(path1+"Merge.shp")

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