GDAL 实现数据空间查询

前言

❝

在GIS开发中，空间查询和属性查询都是常见的基础操作，是每一个GISer都要掌握的必备技能。实现高效的数据查询功能可以提升用户体验，提升数据可视化效率。

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

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

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

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

ifhasDs is None:
print("数据源打开异常，请检查路径！")
returnFalse

# 获取图层
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插件中的投影库版本，需要更换投影库或者升级版本。

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

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