Python GDAL | 教程之：过滤器，简单的空间分析，函数和模块

layer对象有一个方法叫setattributefilter()可以将layer中符合某一条件的feature过滤出来。设定了filter之后就可以用getnextfeature()方法依次取出符合条件的feature了。setattributefilter(none)可以清楚一个filter。例如

>>> layer.getfeaturecount()

42

>>> layer.setattributefilter(“cover = ‘shrubs'”)

>>> layer.getfeaturecount()

6

>>> layer.setattributefilter(None)

>>> layer.getfeaturecount()

42

 

空间过滤器spatial filters

有两种。一种是setspatialfilter()，过滤某一类型的feature，例如参数中填polygon，就是选出layer中的所有polygon

另外还有setspatialfilterrect(, , , )，参数输入四个坐标，可以选中方框内的feature

setspatialfilter(none)一样是清空空间属性过滤器。

例如下面这段代码，layerareas 是polygon，layersites是point

>>> featareas = layerareas.getnextfeature()

>>> poly = featareas.getgeometryref()

>>> layersites.getfeaturecount()

42

>>> layersites.setspatialfilter(poly)

>>> layersites.getfeaturecount()>>> layersites.getfeaturecount()

33

>>> layersites.setspatialfilterrect(460000, 4590000, 490000, 4600000)

>>> layersites.getfeaturecount()

4

>>> layersites.setspatialfilter(none)

>>> layersites.getfeaturecount()

42

还有更复杂的filter，例如执行sql查询语句executesql()，凭借sql的强大功能，可以执行更复杂的任务，例如下面这段代码，就是选择cover类型为grass的feature，并且按id号降序排列。

import ogr,os
os.chdir(r'F:\Python+gdal\7weeks数据\7weeks数据\ospy_data3\ospy3_data')
driver = ogr.GetDriverByName('ESRI Shapefile')
inds = ogr.Open('sites.shp',0)
resultslayer = inds.ExecuteSQL("select * from sites where cover ='grass' order by id desc")
feature = resultslayer.GetNextFeature()
while feature:
    print(feature.GetField('id'))
    feature = resultslayer.GetNextFeature()
inds.ReleaseResultSet(resultslayer)



result = dssites.executesql(“select * from sites where cover = ‘grass’ order by id desc”)

resultfeat = result.getnextfeature()

while resultfeat :

print resultfeat.getfield(‘id’)print resultfeat.getfield(‘id’)

resultfeat = result.getnextfeature()

dssites.releaseresultset(result)

42

40

4

最后一句releaseresultset()是将查询结果释放，在执行下一条sql语句之前一定要先释放。

下面的例子，统计了cover为grass的所有feature的数目

>>> result = dssites.executesql(“select count(*) from sites where cover = ‘grass'”)

>>> result.getfeaturecount()

11

>>> result.getfeature(0).getfield(0)

11

>>> dssites.releaseresultset(result)

列出所有不同的cover类型

resuts3 = inds.ExecuteSQL("select distinct cover from sites")
feature3 = resuts3.GetNextFeature()
while feature3:
    cntlayer = inds.ExecuteSQL("select * from sites where cover = '"+ feature3.GetField('cover')+"'")
    cnt = cntlayer.GetFeatureCount()
    print("%s: %s"%(feature3.GetField('cover'),str(cnt)))
    inds.ReleaseResultSet(cntlayer)
    feature3 = resuts3.GetNextFeature()
inds.ReleaseResultSet(resuts3)


result = ds.executesql(“select distinct cover from sites”)

resultfeat = result.getnextfeature()

while resultfeat:

print resultfeat.getfield(0)

resultfeat = result.getnextfeature()

ds.releaseresultset(result)

shrubs

trees

rocks

grass

bare

water

统计每种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) + ‘ ‘ +print coverfeat.getfield(0) + ‘ ‘ + cntlayer.getfeature(0).getfieldasstring(0)

ds.releaseresultset(cntlayer)

coverfeat = coverlayer.getnextfeature()

ds.releaseresultset(coverlayer)

shrubs 6

trees 11

rocks 6

grass 11

bare 6

water 2

intersect判断两个要素是否相交

poly2.intersect(poly1)

返回0表示不相交，返回1表示相交

disjoint判断两个要素是否不相交

poly2.disjoint(poly1)

返回1表示不相交，返回0表示相交，跟intersect正好相反

touch表示相邻（擦边）

poly2.touches(poly1)

返回0表示不擦边，返回1表示擦边

crosses穿越，一般是一条线穿过一个多边形

poly2.crosses(line)

返回0表示不穿过，返回1表示穿过

within包含，一个要素完全被另一个要素圈起来了

ptb.within(poly1)

返回0表示点在多边形外，返回1表示点在多边形内

contains包含，跟within正好相反

poly1.contains(ptb)

就是把主调对象和参数换一下啦

overlaps重叠，好像只有两个多边形之间才能overlap

poly2.overlaps(poly3)

返回0表示不重叠，返回1表示重叠

下面看看简单的地理数据处理geoprocessing

多边形的：

交：poly3.intersection(poly2)

并：poly3.union(poly2)

差：poly3.difference(poly2)

补：poly3.symmetricdifference(poly2)

geometry的：

.buffer() 给geometry加buffer，就是把点线变成多边形，变粗了

.equal() 两个geometry相等吗？

.distance() 返回两个geometry之间的最短距离

.getenvelope() 信封，有意思，其实就是用一个方框框住这个几何形状，返回四个角的坐标(minx, maxx, miny, maxy)