【GDAL学习】过滤器，简单的空间分析，函数和模块
目录
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
转载自：https://blog.csdn.net/qq_37935516/article/details/84192289
