1.屬性過濾器

>>>import ogr,os
>>>os.chdir('E:/data/GDAL/ospy_data3')
>>>driver=ogr.GetDriverByName('ESRI Shapefile')
>>>ds=driver.Open('sites.shp',0)
>>>ly=ds.GetLayer()
>>>ly.GetFeatureCount()
Out[7]: 42
>>>ly.SetAttributeFilter("cover='shrubs'")
Out[8]: 0
>>>ly.GetFeatureCount()
Out[9]: 6
>>>ly.SetAttributeFilter(None)
Out[10]: 0
>>>ly.GetFeatureCount()
Out[11]: 42
 

2.空間過濾器

>>>import ogr
>>>import os
>>>os.chdir('E:/data/GDAL/ospy_data3')
>>>driver = ogr.GetDriverByName('ESRI Shapefile')
>>>ds1 = driver.Open('cache_towns.shp', 0)
>>>ds2 = driver.Open('sites.shp', 0)
>>>ly1 = ds1.GetLayer()
>>>ly2 = ds2.GetLayer()
>>>fAreas=ly1.GetNextFeature()
>>>poly=fAreas.GetGeometryRef()
>>>ly2.GetFeatureCount()
Out[12]: 42
>>>ly2.SetSpatialFilter(poly)
>>>ly2.GetFeatureCount()
Out[14]: 0
>>>ly2.SetSpatialFilterRect(460000,4590000,490000,4600000)
>>>ly2.GetFeatureCount()
Out[16]: 4
>>>ly2.SetSpatialFilter(None)
>>>ly2.GetFeatureCount()
Out[18]: 42
 

更強大的查詢功能，執行SQL程式碼，查找出cover型別為grass的Feature並以ID降序排列。

最後一句ReleaseResultSet(<result_layer>)是將查詢結果釋放，在執行下一條SQL語句之前一定要先釋放。

>>>import ogr,os
...os.chdir('E:/data/GDAL/ospy_data3')
...driver=ogr.GetDriverByName('ESRI Shapefile')
...ds=driver.Open('sites.shp',0)
...ly=ds.GetLayer()
...ly.GetFeatureCount()
Out[2]: 42
>>>result=ds.ExecuteSQL("select * from sites where cover = 'grass' order by id desc")
>>>resultFeat=result.GetNextFeature()
>>>while resultFeat:
    ...print(resultFeat.GetField('id'))
    ...resultFeat=result.GetNextFeature()
    
42
40
37
31
28
24
17
13
11
10
4
>>>ds.ReleaseResultSet(result)
 

 接上面程式碼。統計cover為grass的所有Feature的數目

>>>result=ds.ExecuteSQL("select count(*) from sites where cover = 'grass'")
>>>result.GetFeatureCount()
Out[8]: 1
>>>result.GetFeature(0).GetField(0)
Out[9]: 11
>>>ds.ReleaseResultSet(result)

 接上面程式碼，列出所有的cover型別

>>>result=ds.ExecuteSQL("select distinct cover from sites")
>>>resultFeat=result.GetNextFeature()
>>>while resultFeat:
    ...print(resultFeat.GetField(0))
    ...resultFeat = result.GetNextFeature()
    
shrubs
trees
rocks
grass
bare
water
>>>ds.ReleaseResultSet(result)

接上面程式碼，統計每種cover型別各有多少個Feature

>>>coverLayer = ds.ExecuteSQL('select distinct cover from sites')
>>>coverFeat = coverLayer.GetNextFeature()
>>>while coverFeat:
    ...cntLayer = ds.ExecuteSQL("select count(*) from sites where cover = '" + coverFeat.GetField(0) +"'")
    ...print(coverFeat.GetField(0)+' '+cntLayer.GetFeature(0).GetFieldAsString(0))
    ...ds.ReleaseResultSet(cntLayer)
    ...coverFeat = coverLayer.GetNextFeature()
    ...
shrubs 6
trees 11
rocks 6
grass 11
bare 6
water 2
>>>ds.ReleaseResultSet(coverLayer)

 3.空間操作

1）Intersect：判斷兩個要素是否相交，Polygon與Polygon，Polygon與Point，Polygon與Line，Line與Point

>>> poly2.Intersect(poly1)
0
>>> poly2.Intersect(poly3)
1
>>> poly2.Intersect(poly2)
1
>>> poly1.Intersect(ptA)
0
>>> poly1.Intersect(ptB)
1
>>> poly1.Intersect(line)
1
>>> poly3.Intersect(line)
1
>>> line.Intersect(ptB)
1

2）Disjoint：判斷兩個要素是否不相交，Polygon與Polygon，Polygon與Point，Polygon與Line，Line與Point

>>> poly2.Disjoint(poly1)
1
>>> poly2.Disjoint(poly3)
0
>>> poly1.Disjoint(ptA)
1
>>> poly1.Disjoint(ptB)
0
>>> poly1.Disjoint(line)
0
>>> poly3.Disjoint(line)
0
>>> line.Disjoint(ptB)

3）Touches：表示相鄰（擦邊），Polygon與Polygon，Polygon與Point，Polygon與Line

>>> poly2.Touches(poly1)
0
>>> poly2.Touches(poly3)
0
>>> poly1.Touches(line)
0
>>> poly1.Touches(ptB)
0
>>> poly3.Touches(line)
1

4）Crosses：穿越，一般是一條線穿過一個多邊形，一般是Polygon與Line

>>> poly2.Crosses(poly1)
0
>>> poly2.Crosses(poly3)
0
>>> poly2.Crosses(line)
1
>>> poly3.Crosses(line)
0
>>> poly1.Crosses(line)
1
>>> line.Crosses(ptB)
0

5）Within：包含於，一個要素完全被另一個要素圈起來了，Polygon與Polygon，Polygon與Point，Polygon與Line

>>> poly3.Within(poly2)
0
>>> line.Within(poly2)
0
>>> ptA.Within(poly1)
0
>>> ptB.Within(poly1)
1
>>> poly1.Within(ptB)
0

6）Contains：包含，跟Within正好相反，一個要素完全包含另一個要素，Polygon與Polygon，Polygon與Point，Polygon與Line

>>> poly3.Contains(poly2)
0
>>> line.Contains(poly2)
0
>>> poly2.Contains(line)
0
>>> poly1.Contains(ptA)
0
>>> poly1.Contains(ptB)
1
>>> ptB.Contains(poly1)
0

7）Overlaps：重疊，只有兩個多邊形之間才能overlap，Polygon與Polygon

>>> poly2.Overlaps(poly1)
0
>>> poly2.Overlaps(poly3)
1
>>> poly2.Overlaps(line)
0
>>> poly3.Overlaps(line)
0
>>> poly1.Overlaps(ptB)
0

用week3中的資料測試一下，下圖中多邊形的FID是18，在cache_towns.shp檔案中；點的FID是40，cover型別是tree，在sites.shp檔案中。

>>>import ogr
...import os
...os.chdir('E:/data/GDAL/ospy_data3')
...driver = ogr.GetDriverByName('ESRI Shapefile')
...ds1 = driver.Open('sites.shp', 0)
...if ds1 is None:
...    print('Cannot open this file:', ds1.name)
...pyLy=ds1.GetLayer()
...treeFe=pyLy.GetFeature(40)
...treePt=treeFe.GetGeometryRef()
...ds2=driver.Open('cache_towns.shp',0)
...if ds2 is None:
...    print('Cannot open this file:',ds2.name)
...tnLy=ds2.GetLayer()
...townFe=tnLy.GetFeature(18)
...townPg=townFe.GetGeometryRef()
>>>townPg.Intersect(treePt)
Out[3]: True
>>>treePt.Intersect(townPg)
Out[4]: True
>>>townPg.Disjoint(treePt)
Out[5]: False
>>>townPg.Touches(treePt)
Out[6]: False
>>>townPg.Crosses(treePt)
Out[7]: False
>>>townPg.Within(treePt)
Out[8]: False
>>>townPg.Contains(treePt)
Out[9]: True

Assignment 3a:Use filters and buffers 

程式碼1：自己寫的，沒有註釋

• Use an attribute filter to restrict cache_towns.shp to Nibley (“name” field) 
• Buffer the Nibley geometry by 1500 
• Use the new geometry with a spatial filter on sites.shp to find all sites within 1500 meters of Nibley
• Print out the "id" value for those sites

# Assignment 3a
import ogr
import os

os.chdir('E:/data/GDAL/ospy_data3')
driver = ogr.GetDriverByName('ESRI Shapefile')
ds1 = driver.Open('cache_towns.shp', 0)
if ds1 is None:
    print('Cannot open this file:', ds1.name)
ly1 = ds1.GetLayer()
ly1.SetAttributeFilter("NAME='Nibley'")
NibleyFeat = ly1.GetFeature(0)
NibleyGeom = NibleyFeat.GetGeometryRef()
bufferGeom = NibleyGeom.Buffer(1500)

ds2 = driver.Open('sites.shp', 0)
if ds2 is None:
    print('Cannot open this file:', ds2.name)
ly2 = ds2.GetLayer()
ly2.SetSpatialFilter(bufferGeom)
# ly2.GetFeatureCount()
siteFeat = ly2.GetNextFeature()
while siteFeat:
    print(siteFeat.GetField('ID'))
    siteFeat.Destroy()
    siteFeat = ly2.GetNextFeature()

ds2.Destroy()
ds1.Destroy()

 程式碼2：參考答案，有詳細註釋

# Assignment 3a
# import modules
import ogr, os, sys

# set the working directory
os.chdir('E:/data/GDAL/ospy_data3')

# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# open sites.shp and get the layer
sitesDS = driver.Open('sites.shp', 0)
if sitesDS is None:
    print('Could not open sites.shp')
    sys.exit(1)
sitesLayer = sitesDS.GetLayer()

# open cache_towns.shp and get the layer
townsDS = driver.Open('cache_towns.shp', 0)
if townsDS is None:
    print('Could not open cache_towns.shp')
    sys.exit(1)
townsLayer = townsDS.GetLayer()

# use an attribute filter to restrict cache_towns.shp to "Nibley"
townsLayer.SetAttributeFilter("NAME = 'Nibley'")

# get the Nibley geometry and buffer it by 1500
nibleyFeature = townsLayer.GetFeature(0)
nibleyGeom = nibleyFeature.GetGeometryRef()
bufferGeom = nibleyGeom.Buffer(1500)

# use bufferGeom as a spatial filter on sites.shp to get all sites
# within 1500 meters of Nibley
sitesLayer.SetSpatialFilter(bufferGeom)

# loop through the remaining features in sites.shp and print their
# id values
siteFeature = sitesLayer.GetNextFeature()
while siteFeature:
    print(siteFeature.GetField('ID'))
    siteFeature.Destroy()
    siteFeature = sitesLayer.GetNextFeature()

# close the shapefiles
sitesDS.Destroy()
townsDS.Destroy()

資料下載：https://download.csdn.net/download/qq_37935516/10791607