GDAL 实现矢量裁剪

前言

❝

矢量数据作为数据处理的半壁江山，在日常工作中涉及到多种操作，矢量数据裁剪尤其具有代表性和重要性，是常用操作，核心原理为从指定数据中提取出目标范围。

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

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

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

1. 开发环境

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

2.Clip裁剪方法

矢量图层Layer对象上具有一个方法Clip可用于数据裁剪。该方法接收裁剪图层和结果图层两个必填参数，剩下三个可选参数。

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

3. 矢量裁剪

定义一个图层裁剪方法LayerClip用于数据提取，该方法接收如下三个参数。

"""
说明：GDAL 空间分析操作
参数：
-sourcePath： 待裁剪Shp数据路径
-clipPath：用于裁剪的Shp数据路径
-resultPath：用于存放裁剪结果的数据路径
"""
def LayerClip(sourcePath,clipPath,resultPath):

在添加数据驱动前使用checkFilePath方法检查文件路径是否存在，之后判断数据驱动是否正常。在打开Shp数据源时传递第二个可选参数，该参数为布尔类型，将其值设置为True意为以可写模式打开数据源。

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

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

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

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

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

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

clipLayer = clipDs.GetLayer(0)

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

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

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

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

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

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

print(f"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~图层裁剪开始~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~")
"""
图层裁剪
"""
# result = sourceLayer.Clip(clipLayer,resultLayer)
result = sourceLayer.Clip(clipLayer,resultLayer,[],GetProgress,GetData)
print(f"裁剪结果：{result}")

# 关闭数据源
sourceDs = clipDs = 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'