矢量資料空間查詢

文章目錄

• ​​開篇​​
• ​​案例一​​

• ​​案例說明​​
• ​​代碼示範​​
• ​​方法總結​​

• ​​案例二​​

• ​​案例說明​​
• ​​代碼示範​​
• ​​方法總結​​

作者：阿振

開篇

在前面四篇部落格中我們主要講了對于空間矢量資料的屬性資料的增删改查，在這篇博文中我們要講解空間查詢–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中。

案例二

案例說明

代碼示範

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，為最小距離。