矢量数据空间查询

---

版权声明：本文为博主原创文章，转载请注明原文出处！

作者：阿振

写作时间：2020-06-14 周天

---

开篇

在前面四篇博客中我们主要讲了对于空间矢量数据的属性数据的增删改查，在这篇博文中我们要讲解空间查询–GIS系统很重要的一项功能。空间查询就是根据地物的空间位置进行查询的一种数据检索方式。比如，我们要查询一条河流经的城市；一个公园内的所有路灯；离当前位置最近的公共卫生间等等都属于常用的空间查询。

OGC简单要素规范定义了空间几何体之间的空间关系，包括Equals，Disjoint，Intersects，Touches，Crosses，Within，Contains，Overlaps，Relate，LocateAlong，LocateBetween。感兴趣的同学可以从OGC官网下载下来看看。

现有的空间数据库例如Oracle Spatial，PostGIS，SQL Server都根据OGC简单要素规范提供了对空间查询的支持，他们有差异地在标准SQL语句中添加了空间关系查询的功能。

本文主要介绍如何使用GDAL库对空间数据进行空间查询，常用的方法可以概括为三大类：

1. 第一类就是使用支持空间查询的SQL语句进行查询，但是这种方式只对某些特定种类的数据源可以使用，有些数据源不一定支持。
2. 第二类是使用GDAL提供的SetSpatialFilter()方法增加过滤条件。但是这种方式只能是选择给定范围的空间地位，类似于Within或者Contains的功能，不能实现其他类型的空间关系查询。
3. 第三类就是读取每个Feature要素包含的Geometry对象，根据Geometry的空间关系手动进行筛选。因为GDAL中的Geometry对象基本上实现了OGC简单要素规范定义的空间关系，所以这种方式最灵活，本文主要介绍如何使用这种方式进行空间查询。

案例一

案例说明

我们现在有省的面状数据以及每个城市的点数据，我们需要找到湖北省内的所有城市。

实现思路是先从省的面状数据中找出湖北省，然后遍历城市的点数据看是否落在湖北省境内。

代码演示

import ogr
ogr.UseExceptions()

ds_province: ogr.DataSource = ogr.Open('../data/Provinces.shp')
l_province: ogr.Layer = ds_province.GetLayer()
# 使用filter()方法找出湖北省
f_hubei: ogr.Feature = next(filter(lambda f: '湖北' in f.GetField('NAME'), l_province))

ds_city: ogr.DataSource = ogr.Open('../data/Cities.shp')
l_city: ogr.Layer = ds_city.GetLayer()
# 使用filter()方法过滤出落在湖北省境内的所有市
selected = filter(lambda f: f.GetGeometryRef().Within(f_hubei.GetGeometryRef()), l_city)
for city in selected:
    print(city.GetField('name'))
del ds_province
del ds_city

方法总结

1. 使用ogr.Open()函数读取Shapefile数据，使用GetLayer()获取当前图层，图层中包含了所有的Feature要素。
2. 使用Python的内置filter()函数对省进行过滤，通过NAME字段找出湖北省。filter()函数的第一个参数是一个自定义函数，第二个参数是一个可迭代对象iterable。该函数会遍历可迭代对象将满足第一个自定义函数的值过滤出来。通过next()方法拿到迭代器的当前值，即湖北省的Feature对象。
3. 继续使用filter()函数对城市的点数据进行筛选，这里通过Feature的GetGeometryRef()方法获得要素代表的几何体，然后调用Geometry的Within()方法判断该城市是否落在湖北省对应的Geometry中。

案例二

案例说明

我们将使用城市的点数据获取离武汉市最近的三座城市。

实现的思路是首先从数据中找到武汉市，然后计算每个城市到武汉市的距离并排序，对排好序的Feature选择前三即可。

代码演示

import ogr
ogr.UseExceptions()

ds: ogr.DataSource = ogr.Open('../data/Cities.shp')
cities: ogr.Layer = ds.GetLayer()
# 使用filter()方法找出武汉市
city: ogr.Feature = next(filter(lambda f: '武汉' in f.GetField('name'), cities))
# 调用ResetReading()方法特别重要，如果不ResetReading的话后面的对Feature的遍历会出错
cities.ResetReading()
# 根据每个市到武汉市的距离进行排序
selected = sorted(cities, key=lambda f: f.GetGeometryRef().Distance(city.GetGeometryRef()))
for i in range(1, 4):
    print(selected[i].GetField('name'))
del ds

方法总结

1. 跟案例一一样，我们使用Python的内置filter()函数对市进行过滤，通过NAME字段找出武汉市。
2. 需要特别注意了，当我们遍历完一遍Layer的Feature以后需要调用ResetReading()对迭代器重新归位，否则后面要继续进行要素遍历的话会出错。
3. 接着我们使用Python内置函数sorted()根据每个城市到武汉市的距离进行排序。sorted()函数包含三个参数（后两个可选），第一个参数是一个可迭代对象iterable，第二个参数是用于自定义排序的函数，第三个参数指定是否逆序。sorted()函数的返回值是一个list对象。
4. 对于距离的计算，我们首先使用GetGeometryRef()函数获得要素对应的空间几何体，然后再使用Geometry对象的Distance()函数进行。
5. 计算完以后我们从第二个元素进行输出，因为第一个元素肯定是武汉市，武汉市到武汉市的距离为0，为最小距离。