GDAL 实现数据属性查询

前言

❝

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

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

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

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

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

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

os.environ['PROJ_LIB'] = r'D:ProgramsPythonPython311Libsite-packagesosgeodataproj'