登录GDAL 实现数据空间查询
^ 关注我，带你一起学GIS ^
前言
❝
在GIS开发中，空间查询和属性查询都是常见的基础操作，是每一个GISer都要掌握的必备技能。实现高效的数据查询功能可以提升用户体验，提升数据可视化效率。
在之前的文章中讲了如何使用GDAL或者ogr2ogr工具将txt以及csv文本数据转换为Shp格式，本篇教程在之前一系列文章的基础上讲解如何使用GDAL实现数据空间查询功能。
如果你还没有看过，建议从以上内容开始。
1. 开发环境
本文使用如下开发环境，以供参考。
时间：2025年
系统：Windows 11
Python：3.11.7
GDAL：3.11.1
2. 空间查询
在GDAL中，有两个图层方法可以用于实现空间查询，分别是SetSpatialFilter和SetSpatialFilterRect。此种查询方式为直接在源数据上操作，返回结果为查询图层。
"""
参数
-geom：用于几何查询的几何对象
"""
SetSpatialFilter(geom)
"""
参数：
-minx：最小x坐标
-miny：最小y坐标
-maxx：最大x坐标
-maxy：最大y坐标
"""
SetSpatialFilterRect(minx,miny,maxx,maxy)
传入None值时重置查询。
sourceLayer.SetSpatialFilter(None)
在进行正式查询前和以前的文章一样，首先获取数据驱动，添加数据源，获取图层数据。定义一个方法SpatialFilter用于实现空间查询，该方法接受两个参数，一个sourcePath传递源数据图层路径，另一个selectPath用于定义查询图层路径。
"""
说明：图层属性过滤
参数：
-sourcePath：待查询 Shp 文件路径
-selectPath：查询 Shp 文件路径
"""
def SpatialFilter(sourcePath,selectPath):
在以下代码中完成图层数据的读取操作。
# 注册所有驱动
ogr.RegisterAll()
# 添加数据驱动
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
checkFilePath(sourcePath,shpDriver)
checkFilePath(selectPath,shpDriver)
# 打开数据源
sourceDs = shpDriver.Open(sourcePath)
selectDs = shpDriver.Open(selectPath)
hasDs = sourceDs and selectDs
if hasDs is None:
print("数据源打开异常，请检查路径！")
return False
# 获取图层
sourceLayer = sourceDs.GetLayer(0)
selectLayer = selectDs.GetLayer(0)
下面对几何查询以及矩形查询两种实现方式进行介绍。
2.1. 几何查询
通过图层方法GetNextFeature获取第一个要素对象。
# 获取查询几何对象
queryFeat = selectLayer.GetNextFeature()
queryGeom = queryFeat.GetGeometryRef()
print(f"查询要素id：{queryFeat.GetField('Id')}")
print(f"查询几何对象：{queryGeom}")
将几何对象传入SetSpatialFilter方法进行空间过滤。
# 空间过滤
sourceLayer.SetSpatialFilter(queryGeom)
查询完成之后传入None值结束查询。
# 结束查询
sourceLayer.SetSpatialFilter(None)
以下为几何查询部分代码。
# 获取查询几何对象
queryFeat = selectLayer.GetNextFeature()
queryGeom = queryFeat.GetGeometryRef()
print(f"查询要素id：{queryFeat.GetField('Id')}")
print(f"查询几何对象：{queryGeom}")
# 获取要素数量
featureCount = sourceLayer.GetFeatureCount()
print(f"所有要素数量：{featureCount}")
# 空间过滤
sourceLayer.SetSpatialFilter(queryGeom)
queryFeatCount = sourceLayer.GetFeatureCount()
print(f"查询要素数量：{queryFeatCount}")
# 结束查询
sourceLayer.SetSpatialFilter(None)
finalFeatCount = sourceLayer.GetFeatureCount()
print(f"重置查询后要素数量：{finalFeatCount}")
print("n~~~~~~~方式一：结束几何查询~~~~~~~")
若有兴趣，还可以创建自定义几何对象进行空间查询。
# 自定义Geometry查询对象
customGeom = ogr.Geometry(ogr.wkbPolygon)
ringGeom = ogr.Geometry(ogr.wkbLinearRing)
ringGeom.AddPoint(102.884350,32.501570)
ringGeom.AddPoint(105.025865,74.974949)
ringGeom.AddPoint(50.417235,55.701314)
ringGeom.AddPoint(50.417235,32.501570)
ringGeom.CloseRings()
customGeom.AddGeometry(ringGeom)
# 获取要素数量
featureCount = sourceLayer.GetFeatureCount()
print(f"所有要素数量：{featureCount}")
# 空间过滤
sourceLayer.SetSpatialFilter(customGeom)
queryFeatCount = sourceLayer.GetFeatureCount()
print(f"查询要素数量：{queryFeatCount}")
# 结束查询
sourceLayer.SetSpatialFilter(None)
finalFeatCount = sourceLayer.GetFeatureCount()
print(f"重置查询后要素数量：{finalFeatCount}")
print("n~~~~~~~方式一：结束几何查询~~~~~~~")
如下为几何查询输出结果：
如下为该要素在ArcGIS中查询显示结果。
2.2. 矩形查询
矩形查询主要调用SetSpatialFilterRect方法，传入x、y坐标即可。
# 获取要素数量
featureCount2 = sourceLayer.GetFeatureCount()
print(f"所有要素数量：{featureCount2}")
# 空间过滤
sourceLayer.SetSpatialFilterRect(-16.9994,23.1235,42.9111,49.5852,)
queryFeatCount2 = sourceLayer.GetFeatureCount()
如下为矩形查询输出结果：
3. 注意事项
在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开发 相关内容，欢迎关注
