前言
❝在GIS开发中，属性查询是非常普遍的操作，这是每一个GISer都要掌握的必备技能。实现高效的数据查询功能可以提升用户体验，提升数据可视化效率。
在之前的文章中讲了如何使用GDAL或者ogr2ogr工具将txt以及csv文本数据转换为Shp格式，本篇教程在之前一系列文章的基础上讲解如何使用GDAL实现空间数据属性查询功能。
GDAL 简介 GDAL 下载安装 GDAL 开发起步 GDAL 实现 GIS 数据读取转换（全）
如果你还没有看过，建议从以上内容开始。
1. 开发环境
本文使用如下开发环境，以供参考。
时间：2025年
系统：Windows 11
Python：3.11.7
GDAL：3.11.1
2. 属性查询
定义一个方法AttributeFilter用于实现空间数据属性查询，该方法接受一个shpPath路径参数。
"""
说明：图层属性过滤
参数：
-shpPath：Shp 文件路径
"""
def AttributeFilter(shpPath):
还是老规矩，第一步添加Shp数据驱动。
# 注册所有驱动
ogr.RegisterAll()
# 添加数据驱动
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
打开数据源，进行属性过滤。图层对象有个方法SetAttributeFilter可以进行简单的属性查询，只需传入字段名以及字段值，查询结束后将查询条件设置为None退出查询。使用此种方式是直接在源数据图层上进行操作。
with shpDriver.Open(shpPath) as ds:
ifds is None:
print("数据源打开异常，请检查路径！")
returnFalse
# 获取图层
layer = ds.GetLayer(0)# Shp图层默认值为0
# 获取要素数量
featureCount = layer.GetFeatureCount()
print(f"所有要素数量：{featureCount}")
# "学校"要素数量
layer.SetAttributeFilter("type='学校'")
featureCount1 = layer.GetFeatureCount()
print(f"'学校'要素数量：{featureCount1}")
layer.SetAttributeFilter(None)
featureCount2 = layer.GetFeatureCount()
print(f"退出查询后要素数量：{featureCount2}")
属性表中总共有45条记录，查询记录与其保持一致。
除了上面介绍的SetAttributeFilter方法外，还有另外一种方式，就是SQL语句查询。数据源对象具有一个方法ExecuteSQL用于查询指定数据，大多情况下只需要传入一条SQL语句即可。查询结束之后调用ReleaseResultSet方法退出查询。使用此种方式将会返回一个满足查询条件的新图层。
# 复杂查询
# queryLayer = ds.ExecuteSQL("SELECT * FROM geoPoint WHERE TYPE='饭店'")
# 升序
# queryLayer = ds.ExecuteSQL("SELECT * FROM geoPoint WHERE TYPE='饭店' ORDER BY Id ASC")
# 降序
queryLayer = ds.ExecuteSQL("SELECT * FROM geoPoint WHERE TYPE='饭店' ORDER BY Id DESC")
queryFeatCount = queryLayer.GetFeatureCount()
print(f"nSQL查询要素数量：{queryFeatCount}n")
forfeatureinqueryLayer:
# 获取字段
field = feature.GetField("id")
# print(f"要素Id：{field}")
fieldValue = feature.GetField("type")
foriinrange(fieldCount):
fieldDefn = featureDefn.GetFieldDefn(i)
fieldName = fieldDefn.GetName()
fieldType = fieldDefn.GetType()
fieldTypeName = fieldDefn.GetTypeName()
fieldValue = feature.GetField("type")
print(f"id：{field}，value：{fieldValue}n")
ds.ReleaseResultSet(queryLayer)
ds = None
选择类型为”饭店“的POI点，并以Id字段降序排列。
SELECT * FROM geoPoint WHERE TYPE='饭店'ORDER BY Id DESC
查询结果显示如下：与实际数据显示结果一致。
还可以统计以下每种类型要素的数量。
# 统计每种类型的要素数量
finalFeatCount = layer.GetFeatureCount()
print(f"n最终查询要素数量：{finalFeatCount}n")
countLayer = ds.ExecuteSQL("SELECT DISTINCT type FROM geoPoint")
newFeature = countLayer.GetNextFeature()
whilenewFeature:
# countTypeLayer = ds.ExecuteSQL("SELECT COUNT(*) FROM geoPoint WHERE type='"+newFeature.GetField('type')+"'")
countTypeLayer = ds.ExecuteSQL("SELECT COUNT(*) FROM geoPoint WHERE type='"+newFeature.GetField(0)+"'")
print(f"{newFeature.GetField(0)} {countTypeLayer.GetFeature(0).GetFieldAsString(0)}")
ds.ReleaseResultSet(countTypeLayer)
newFeature = countLayer.GetNextFeature()
ds.ReleaseResultSet(countLayer)
查询结果显示如下：最后不要忘记关闭数据源。
ds = None
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:ProgramsPythonPython311Libsite-packagesosgeodataproj'
