GDAL 实现矢量合并

前言

❝

矢量数据作为数据处理的半壁江山，在日常工作中涉及到多种操作，矢量数据合并也是一项常用操作，该功能涉及到两个及以上的数据源在几何对象与属性对象之间的合并操作。

在之前的文章中讲了如何使用GDAL或者ogr2ogr工具将txt以及csv文本数据转换为Shp格式，本篇教程在之前一系列文章的基础上讲解如何使用GDAL 实现矢量合并。

• GDAL 简介
• GDAL 下载安装
• GDAL 开发起步
• GDAL 实现 GIS 数据读取转换（全）

如果你还没有看过，建议从以上内容开始。

1. 开发环境

本文使用如下开发环境，以供参考。

2.Union合并方法

矢量图层Layer对象上具有一个方法Union可用于数据裁剪。该方法接收合并图层和结果图层两个必填参数，剩下三个可选参数，最后返回一个整型值。

"""
说明：GDAL 矢量合并
参数：
-unionLayer：用于裁剪的矢量图层
-resultLayer：生成结果数据图层
"""
def Union(unionLayer,resultLayer):

3. 矢量图层合并

老规矩定义一个图层合并方法LayerUnion用于数据合并，该方法接收如下三个参数。

"""
说明：GDAL 图层合并操作
参数：
-sourcePath： 源文件Shp数据路径
-unionPath：用于合并的Shp数据路径
-resultPath：用于存放合并结果的数据路径
"""
def LayerUnion(sourcePath,unionPath,resultPath):

然后添加数据驱动，只不过在添加数据驱动前使用checkFilePath方法检查文件路径是否存在，之后判断数据驱动是否正常。

# 检查文件是否存在
checkFilePath(sourcePath)
checkFilePath(unionPath)
checkFilePath(resultPath)

# 添加数据驱动
shpDriver = ogr.GetDriverByName("ESRI Shapefile")

# 检查数据驱动是否正常
checkDriver(shpDriver)

在打开Shp数据源时传递第二个可选参数，该参数为布尔类型，将其值设置为True意为以可写模式打开数据源。

"""
以可写模式打开数据源
"""

# 获取数据源
# sourceDs = shpDriver.Open(sourcePath)
# unionDs = shpDriver.Open(unionPath)
sourceDs = shpDriver.Open(sourcePath,True)
clipDs = shpDriver.Open(unionPath,True)

以源数据图层sourceLayer属性结构创建结果图层，获取源图层空间参考以及几何数据类型信息，并创建结果图层属性字段。

# 获取图层
sourceLayer = sourceDs.GetLayer(0)
srs = sourceLayer.GetSpatialRef()
geomType = sourceLayer.GetGeomType()

unionLayer = unionDs.GetLayer(0)

# 创建输出数据源
resultDs = shpDriver.CreateDataSource(resultPath)

# 根据源数据结构创建输出图层
resultLayer = resultDs.CreateLayer("resultLayer",srs,geomType)

# 添加数据结构
featureDefn = sourceLayer.GetLayerDefn()
fieldCount = featureDefn.GetFieldCount()

# 创建要素
# resultFeature = ogr.Feature(featureDefn)

# 添加字段（复制源图层字段定义）
foriinrange(fieldCount):
fieldDefn = featureDefn.GetFieldDefn(i)
resultLayer.CreateField(fieldDefn)

在源数据图层上调用方法Union进行图层合并，返回值为一个整型数值。合并完成后关闭所有数据源并返回结果。

print(f"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~图层合并开始~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
"""
图层裁剪
"""
result = sourceLayer.Union(unionLayer,resultLayer,[],GetProgress,GetData)
print(f"裁剪结果：{result}")

# 关闭数据源
sourceDs = unionDs = resultDs = None
print(f"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~图层合并结束~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
returnresult

GetProgress为获取裁剪进度函数。

"""
说明：定义进度获取函数
参数：
-complete：完成度（0.0-1.0）
-message：进度消息
-userData：用户数据
"""
def GetProgress(complete,message,userData):
print(f"进度：{complete * 100:.1f}% - {message}")
# 返回1表示继续处理，返回0表示取消
return1

GetData为获取裁剪数据函数。

def GetData(data):
print("裁剪数据：")

数据合并信息输出如下。

裁剪数据在ArcGIS中显示如下。

4. 注意事项

在windows开发环境中同时安装GDAL与PostGIS，其中投影库PROJ的环境变量指向PostGIS的安装路径，在运行GDAL程序时，涉及到要素、几何与投影操作时会导致异常。具体意思为GDAL不支持PostGIS插件中的投影库版本，需要更换投影库或者升级版本。

解决办法为修改PROJ的环境变量到GDAL支持的版本或者在GDAL程序开头添加以下代码：

os.environ['PROJ_LIB'] = r'D:\Programs\Python\Python311\Libsite-packages\osgeo\data\proj'