目錄

• • 寫在前面
  • 1. 屬性過濾條件
  • 2. 空間過濾條件
  • 3. 使用SQL創(chuàng)建臨時圖層
  • 4. 利用過濾條件

寫在前面

??過濾條件可以將不想要的要素拋棄，通過過濾條件可以選出符合特定條件的要素，也可以通過空間范圍限定要素，這樣就可以簡單地處理感興趣的數據。

1. 屬性過濾條件

??過濾條件需要一個條件語句，類似于SQL語句中的Where子句。如：

‘Population < 50000’ ‘Population = 50000’ ‘Name = “Tokyo”’

??注：比較字符串時，需要用引號將字符串包住，并且保證它們與包住整個查詢字符串的引號不同。否則，將會導致語法錯誤。

??測試某個對象不等于另一個值，可以使用！=或<>，也可以使用AND或者OR語句：

'(Population > 25000) AND (Population < 50000)' '(Population > 50000) OR (Place_type = "County Seat")'

語句1：選出人口大于25000，但少于50000的要素；
語句2：選出人口超過50000，或者Place_type 為 County Seat 的要素（或同時滿足）。

??使用NOT來否定條件，NULL用于指示一個空值或屬性表中沒有數據值：

'(Population < 50000) OR NOT (Place_type = "County Seat")' 'County NOT NULL'

語句1：選出人口小于50000，或Place_type 不是 County Seat 的要素（或同時滿足）；
語句2：選出County 不為空的要素。

??檢查某個值是否在另外兩個值之間，用 BETWEEN ：

'Population BETWEEN 25000 AND 50000' '(Population > 25000) AND (Population < 50000)'

語句1：選出人口在25000和50000之間的要素；
語句2：也是選出人口在25000和50000之間的要素。

??檢驗某個值是否與其他多個不同值相同：

'Type_code IN (4, 3, 7)' '(Type_code = 4) OR (Type_code = 3) OR (Type_code = 7)'

??也適用于字符串：

'Place_type IN ("Populated Place", "County Seat")'

??下劃線匹配任何單個字符，百分號匹配任意數量的字符，使用LIKE進行字符串匹配（不區(qū)分大小寫）：

'Name LIKE "%Seattle%"'

??使用LIKE操作符的匹配實例：

方式匹配不匹配

_eattle	Seattle	Seattle WA
Seattle%	Seattle, Seattle WA	North Seattle
%Seattle%	Seattle, Seattle WA, North Seattle	Tacoma
Sea%le	Seattle	Seattle WA
Sea_le	Seatle (note misspelling)	Seattle

??使用Python交互式窗口測試，需要加載VectorPlotter類來交互式地繪制選擇的結果：

>>> import os >>> import sys >>> from osgeo import ogr >>> import ospybook as pb >>> import matplotlib.pyplot as plt >>> from ospybook.vectorplotter import VectorPlotter >>> data_dir = r'E:\Google chrome\Download' >>> vp = VectorPlotter(True) >>> ds = ogr.Open(os.path.join(data_dir, 'global')) >>> lyr = ds.GetLayer('ne_50m_admin_0_countries')# Get the countries shapefile layer >>> vp.plot(lyr, fill=False) # fill=False表示只繪制出國家邊界 >>> vp.plot(lyr, fill=False) >>> pb.print_attributes(lyr, 4, ['name'], geom=False) # 查看圖層屬性FID name 0 Aruba 1 Afghanistan 2 Angola 3 Anguilla 4 of 241 features

??查看亞洲包含的要素：
??SetAttributeFilter（） 函數

import os import sys from osgeo import ogr import ospybook as pb from ospybook.vectorplotter import VectorPlotterdata_dir = r'E:\Google chrome\Download'ds = ogr.Open(os.path.join(data_dir, 'global')) lyr = ds.GetLayer('ne_50m_admin_0_countries')asia=lyr.SetAttributeFilter("continent = 'Asia'") print(asia) # asia 為0時，表示查詢語句被執(zhí)行成功count_asia=lyr.GetFeatureCount() print(count_asia)# 亞洲總共有多少個記錄vp = VectorPlotter(True) vp.plot(lyr, 'y') vp.draw() 0 53

??繪制結果：

??查看屬性過濾器效果：

>>> pb.print_attributes(lyr, 4, ['name'], geom=False)FID name 1 Afghanistan 7 United Arab Emirates 9 Armenia 17 Azerbaijan 4 of 53 features >>> lyr.GetFeature(2).GetField('name') # FID并沒有改變 'Angola'

??當設置另一個屬性過濾條件，它不是創(chuàng)建當前過濾圖層要素的子集，而是將過濾條件應用于整個圖層：

lyr.SetAttributeFilter('continent = "South America"') vp.plot(lyr, 'b') 5 241

??可以同時使用空間和屬性過濾條件優(yōu)化結果，后面會提到。
清除屬性過濾條件來獲取所有要素：

>>> lyr.SetAttributeFilter(None) 0 >>> lyr.GetFeatureCount() 241

2. 空間過濾條件

??空間過濾條件可以使用空間范圍，而不是屬性值來篩選要素，可以用來選擇位于另一個要素內或邊界框內部的要素。
??使用天然地球shapefile文件來選擇出中國城市。首先，打開數據源文件，使用屬性過濾條件篩選中國，并獲得對應的要素記錄和幾何對象：

>>> from osgeo import ogr >>> from ospybook.vectorplotter import VectorPlotter >>> import ospybook as pb >>> ds=ogr.Open(r"E:\Google chrome\Download\global") >>> country_lyr = ds.GetLayer('ne_50m_admin_0_countries') >>> vp = VectorPlotter(True) >>> vp.plot(country_lyr, fill=False) >>> country_lyr.SetAttributeFilter('name = "China"') 5 >>> feat = country_lyr.GetNextFeature() >>> China = feat.geometry().Clone()

??打開居住區(qū)圖層，將所有城市數據表示為紅點：

>>> city_lyr = ds.GetLayer('ne_50m_populated_places') >>> city_lyr.GetFeatureCount() 1249 >>> vp.plot(city_lyr, 'r.')

# 選取中國城市 # 打開數據源文件夾獲取國家邊界圖層，然后使用屬性過濾條件篩選出中國，并獲得對應的要素記錄和幾何對象 # 使用中國邊界，根據空間過濾條件，篩選出居住區(qū)圖層中的中國城市。from osgeo import ogr from ospybook.vectorplotter import VectorPlotter import ospybook as pb# 打開圖層 ds=ogr.Open(r"E:\Google chrome\Download\global") lyr=ds.GetLayer('ne_50m_admin_0_countries')vp=VectorPlotter(True) vp.plot(lyr,fill=False)# 根據屬性過濾條件篩選出中國 """ 屬性過濾條件篩選后會返回一個要素，lyr.GetNextFeature()獲得這個要素 Clone() 克隆這個要素，這樣即使要素在內存中被刪除后，仍然可以使用""" country = lyr.SetAttributeFilter("name = 'China'") feat = lyr.GetNextFeature() China = feat.geometry().Clone()# 加載城市圖層，并繪制 city_lyr=ds.GetLayer('ne_50m_populated_places') vp.plot(city_lyr,'r.')# 根據空間過濾條件，篩選出中國的城市 country_filter=city_lyr.SetSpatialFilter(China) country_count=city_lyr.GetFeatureCount() print(country_count) vp.plot(city_lyr,'bo') vp.draw()

??使用空間過濾條件和屬性查詢函數組合來優(yōu)化選擇，找出人口在100萬以上的城市，并用正方形表示：

# 找出中國人口超過100萬的城市 country_100 = city_lyr.SetAttributeFilter('pop_min > 1000000') country_100_count = city_lyr.GetFeatureCount() vp.plot(city_lyr, 'rs') vp.draw()

??一共有95個城市人口數超過100萬人口：

??查看世界上有多少個城市人口數大于100萬：

city_lyr.SetSpatialFilter(None) city_lyr.GetFeatureCount() vp.plot(city_lyr, 'm^', markersize=8) vp.draw()# 246

??繪制結果：

SetSpatialFilterRect(minx, miny, maxx, maxy): 創(chuàng)建矩形進行空間過濾。

country_lyr.SetSpatialFilterRect(110, -50, 160, 0) vp.plot(country_lyr, 'y') vp.draw()

??繪制結果：

??繪制澳大利亞的最大矩形可以作為全球國界圖層的一個空間過濾條件。

3. 使用SQL創(chuàng)建臨時圖層

??主要使用EcecuteSQL函數對數據源進行更復雜的查詢，作用與數據源，而不是圖層，可以操作多個圖層?？臻g過濾條件可以作為可選項。
EcecuteSQL（statement，[spatialFilter]，[dialect]）

statement：要使用的SQL語句

spatialFilter：可選項，針對結果做空間過濾，默認沒有過濾

dialect用于說明SQL語句用什么標準的字符串，可用選項為：OGRSQL和SQLite。默認使用SQL標準。

??1. 按照人口數降序排列返回結果：

ds = ogr.Open(os.path.join(data_dir, 'global')) sql = '''SELECT ogr_geom_area as area, name, pop_estFROM 'ne_50m_admin_0_countries' ORDER BY POP_EST DESC''' lyr = ds.ExecuteSQL(sql) pb.print_attributes(lyr, 3) # 此處輸出報錯

??未能解決的問題：

??init文件顯示：

def print_attributes(lyr_or_fn, n=None, fields=None, geom=True, reset=True):"""Print attribute values in a layer.lyr_or_fn - OGR layer object or filename to datasource (will use 1st layer)n - optional number of features to print; default is allfields - optional list of case-sensitive field names to print; defaultis allgeom - optional boolean flag denoting whether geometry type is printed;default is Truereset - optional boolean flag denoting whether the layer should be reset- to the first record before printing; default is True"""lyr, ds = _get_layer(lyr_or_fn)if reset:lyr.ResetReading()n = n or lyr.GetFeatureCount()geom = geom and lyr.GetGeomType() != ogr.wkbNonefields = fields or [field.name for field in lyr.schema]data = [['FID'] + fields]if geom:data[0].insert(1, 'Geometry')feat = lyr.GetNextFeature()while feat and len(data) <= n:data.append(_get_atts(feat, fields, geom))feat = lyr.GetNextFeature()lens = map(lambda i: max(map(lambda j: len(str(j)), i)), zip(*data))format_str = ''.join(map(lambda x: '{{:<{}}}'.format(x + 4), lens))for row in data:try:print(format_str.format(*row))except UnicodeEncodeError:e = sys.stdout.encodingprint(codecs.decode(format_str.format(*row).encode(e, 'replace'), e))print('{0} of {1} features'.format(min(n, lyr.GetFeatureCount()), lyr.GetFeatureCount()))if reset:lyr.ResetReading()

??shapefile文件沒有任何內置的SQL引擎，若要查詢的數據源自身具有SQL引擎支持，則會使用原生的SQL版本。

??2. 使用SQLite版本的SQL標準從數據庫中獲取信息：

ds = ogr.Open(os.path.join(data_dir, 'global','natural_earth_50m.sqlite')) sql = '''SELECT geometry, area(geometry) AS area, name, pop_estFROM countries ORDER BY pop_est DESC LIMIT 3''' lyr = ds.ExecuteSQL(sql) pb.print_attributes(lyr)

??3. 使用函數連接多個圖層的屬性：

# ne_50m_populated_places pp：圖層重命名 # ON pp.adm0_a3 = c.adm0_a3 基于公共字段連接兩表 ds = ogr.Open(os.path.join(data_dir, 'global')) sql = '''SELECT pp.name AS city, pp.pop_min AS city_pop,c.name AS country, c.pop_est AS country_popFROM ne_50m_populated_places ppLEFT JOIN ne_50m_admin_0_countries cON pp.adm0_a3 = c.adm0_a3WHERE pp.adm0cap = 1''' lyr = ds.ExecuteSQL(sql) pb.print_attributes(lyr, 3, geom=False)

??4. 查看多個圖層相關數據（城市人口和全國人口）：

ds = ogr.Open(os.path.join(data_dir, 'global')) sql = '''SELECT pp.name AS city, pp.pop_min AS city_pop,c.name AS country, c.pop_est AS country_popFROM ne_50m_populated_places ppLEFT JOIN ne_50m_admin_0_countries cON pp.adm0_a3 = c.adm0_a3WHERE pp.adm0cap = 1 AND c.continent = "Asia"''' lyr = ds.ExecuteSQL(sql) pb.print_attributes(lyr, 3)

??5. 使用SQLite標準查看shapefile數據源（SQLite標準也適用于其他數據源）：

lyr = ds.ExecuteSQL(sql, dialect='SQLite')

??6. 合并幾何對象：

# Plot the counties in California. ds = ogr.Open(os.path.join(data_dir, 'US')) sql = 'SELECT * FROM countyp010 WHERE state = "CA"' lyr = ds.ExecuteSQL(sql) vp = VectorPlotter(True) vp.plot(lyr, fill=False)sql = 'SELECT st_union(geometry) FROM countyp010 WHERE state = "CA"' lyr = ds.ExecuteSQL(sql, dialect='SQLite') vp.plot(lyr, 'w')

4. 利用過濾條件

??1. 創(chuàng)建所需要的數據：

ds = ogr.Open(os.path.join(data_dir, 'global'), 1)in_lyr = ds.GetLayer('ne_50m_populated_places') in_lyr.SetAttributeFilter("FEATURECLA = 'Admin-0 capital'")out_lyr = ds.CopyLayer(in_lyr, 'capital_cities2') out_lyr.SyncToDisk()

??2. 只選擇需要的要素：

sql = """SELECT NAME, ADM0NAME FROM ne_50m_populated_placesWHERE FEATURECLA = 'Admin-0 capital'""" in_lyr2 = ds.ExecuteSQL(sql) out_lyr2 = ds.CopyLayer(in_lyr2, 'capital_cities3') out_lyr2.SyncToDisk()

總結

以上是生活随笔為你收集整理的Python地理数据处理 六：使用OGR过滤数据的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。