老样子最近有项目需要将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文件，这个就是合并之后的结果。