登录
GDAL 实现空间分析
^ 关注我，带你一起学GIS ^
前言
❝
空间分析基于空间关系，常见的如叠加、求差、取反等。在GIS开发中，空间分析作为其灵魂支柱，具有重要意义，是每一个GISer都必须要掌握的技能。
在之前的文章中讲了如何使用GDAL或者ogr2ogr工具将txt以及csv文本数据转换为Shp格式，本篇教程在之前一系列文章的基础上讲解如何使用GDAL 实现空间分析。
如果你还没有看过，建议从以上内容开始。
1. 开发环境
本文使用如下开发环境，以供参考。
时间：2025年
系统：Windows 11
Python：3.11.7
GDAL：3.11.1
2. 空间分析
定义一个方法SpatialAnalysis用于判断要素间空间关系，该方法接收两个参数，其中sourcePath指向源数据路径，resultPath指向结果文件路径。
"""说明：GDAL 空间分析操作参数： -sourcePath：Shp 数据路径 -resultPath：生成结果数据路径"""def SpatialAnalysis(sourcePath,resultPath):
接着是基础的常规操作，添加数据驱动，获取数据源，获取目标图层。
# 检查文件是否存在checkFilePath(sourcePath)# 添加数据驱动shpDriver = ogr.GetDriverByName("ESRI Shapefile")# 获取数据源shpDs = shpDriver.Open(sourcePath)# 获取图层shpLayer = shpDs.GetLayer(0)srs = shpLayer.GetSpatialRef()geomType = shpLayer.GetGeomType()
方法checkFilePath用于检查文件路径是否正常，接受一个路径参数。
"""说明：检查文件路径是否正常参数： -filePath：文件数据路径"""def checkFilePath(filePath): if os.path.exists(filePath): print("文件数据路径存在") else: print("文件数据路径不存在，请检查！")
本文基于以下数据进行空间分析操作。
对涉及面的操作分别获取面0和面1两个要素。GDAL中的几何处理操作，如叠加、缓冲、求差等均返回一个新的几何对象。
"""获取第一个图层要素"""# feature = shpLayer.GetNextFeature()# 根据FId获取指定要素feature1 = shpLayer.GetFeature(0)print(f"n面状要素~~~~~Id：{feature1.GetField('Id')}n")# 获取要素IdfeatId1 = feature1.GetField("Id")# 获取几何要素geom1 = feature1.GetGeometryRef()# print(f"面状几何：{geom1}")""" 获取第二个要素"""feature2 = shpLayer.GetFeature(1)geom2 = feature2.GetGeometryRef()featId2 = feature2.GetField("Id")
3. 空间叠加
Geometry对象上具有一个方法Intersection用于空间叠加操作，该方法接收另一个几何对象，并返回一个新的Geometry对象。
""" 空间操作之【叠加】"""resultGeom = geom1.Intersection(geom2)# print(f"面要素【{featId1}】与 面要素【{featId2}】相交结果：{resultGeom}")
如下面0和面1相交结果为图中高亮部分。
4. 空间缓冲
Geometry对象上具有一个方法Buffer用于空间缓冲区生成，该方法接收另一个几何对象，并返回一个新的Geometry对象。
""" 空间操作之【缓冲】"""resultGeom = geom1.Buffer(1)# print(f"面要素【{featId1}】缓冲结果：{resultGeom}")
如下面0缓冲结果为图中高亮部分。
如下线段缓冲结果为图中高亮部分。
如下点缓冲结果为图中高亮部分。
5. 空间联合
Geometry对象上具有一个方法Union用于合并几何对象，该方法接收另一个几何对象，并返回一个新的Geometry对象。
""" 空间操作之【联合】"""resultGeom = geom1.Union(geom2)# print(f"面要素【{featId1}】联合结果：{resultGeom}")
如下面0和面1联合结果为图中高亮部分。
6. 空间求差
Geometry对象上具有一个方法Difference用于求取目标几何对象差集，该方法接收另一个几何对象，并返回一个新的Geometry对象。
""" 空间操作之【求差】"""resultGeom = geom1.Difference(geom2)# print(f"面要素【{featId1}】联合结果：{resultGeom}")
如下面0和面1求差结果为图中高亮部分。
7. 对称差集
Geometry对象上具有一个方法SymmetricDifference用于几何对象交集取反操作，该方法接收另一个几何对象，并返回一个新的Geometry对象。
""" 空间操作之【求对称差集】"""resultGeom = geom1.SymmetricDifference(geom2)print(f"面要素【{featId1}】联合结果：{resultGeom}")
如下面0和面1求对称差集结果为图中高亮部分。
还有一些几何操作，如求等、求距离等内容本教程并未涉及，留给感兴趣的读者自行实践。
8. 创建新图层
几何运算完成之后，需要将分析结果保存到新的Shp图层。通过数据驱动方法CreateDataSource创建Shp数据源，再使用数据源方法CreateLayer创建图层即可。
""" 创建新图层"""targetDs = shpDriver.CreateDataSource(resultPath)targetLayer = targetDs.CreateLayer("result_layer",srs,ogr.wkbPolygon)featureDefn = shpLayer.GetLayerDefn()tarFeat = ogr.Feature(featureDefn)fieldCount = featureDefn.GetFieldCount()tarFeat.SetGeometry(resultGeom)for i in range(fieldCount): fieldDefn = featureDefn.GetFieldDefn(i) name = fieldDefn.GetName() value = feature1.GetField(i) print(f"字段名称：{name}，字段值：{value}") # tarFeat.SetField(i,value) tarFeat.SetField(name,value)targetLayer.CreateFeature(tarFeat)
9. 注意事项
在windows开发环境中同时安装GDAL与PostGIS，其中投影库PROJ的环境变量指向PostGIS的安装路径，在运行GDAL程序时，涉及到要素、几何与投影操作时会导致异常。具体意思为GDAL不支持PostGIS插件中的投影库版本，需要更换投影库或者升级版本。
RuntimeError: PROJ: proj_identify: D:Program FilesPostgreSQL13sharecontribpostgis-3.5projproj.db contains DATABASE.LAYOUT.VERSION.MINOR = 2 whereas a number >= 5 is expected. It comes from another PROJ installation.
解决办法为修改PROJ的环境变量到GDAL支持的版本或者在GDAL程序开头添加以下代码：
os.environ['PROJ_LIB'] = r'D:\Programs\Python\Python311\Libsite-packages\osgeo\data\proj'
视频效果
❝
OpenLayers示例数据下载，请在公众号后台回复：ol数据
全国信息化工程师－GIS 应用水平考试资料，请在公众号后台回复：GIS考试
❝
GIS之路 公众号已经接入了智能 助手，可以在对话框进行提问，也可以直接搜索历史文章进行查看。
都看到这了，不要忘记点赞、收藏 + 关注 哦 ！
本号不定时更新有关 GIS开发 相关内容，欢迎关注
