4. GDAL python教程(3)——过滤器,简单的空间分析,函数和模块?

4.1. 属性过滤器Attribute filters?

Layer对象有一个方法叫SetAttributeFilter(<where_clause>)可以将Layer中符合某一条件的Feature过滤出来。设定了Filter之后就可以用GetNextFeature()方法依次取出符合条件的Feature了。SetAttributeFilter(None)可以清楚一个Filter。例如

>>> layer.GetFeatureCount()
42
>>> layer.SetAttributeFilter("cover = 'shrubs'")
>>> layer.GetFeatureCount()
6
>>> layer.SetAttributeFilter(None)
>>> layer.GetFeatureCount()
42

4.2. 空间过滤器Spatial filters?

有两种。一种是SetSpatialFilter(<geom>),过滤某一类型的Feature,例如参数中填Polygon,就是选出Layer中的所有Polygon。

另外还有SetSpatialFilterRect(<minx>, <miny>, <maxx>, <maxy>),参数输入四个坐标,可以选中方框内的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>),凭借SQL的强大功能,可以执行更复杂的任务, 例如下面这段代码,就是选择cover类型为grass的Feature,并且按id号降序排列。

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(<result_layer>)是将查询结果释放,在执行下一条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类型

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

4.3. 统计每种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

4.4. 空间操作?

4.4.1. Intersect判断两个要素是否相交?

poly2.Intersect(poly1)

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

4.4.2. Disjoint判断两个要素是否不相交?

poly2.Disjoint(poly1)

返回1表示不相交,返回0表示相交,跟Intersect正好相反

4.4.3. Touch表示相邻(擦边)?

poly2.Touches(poly1)

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

4.4.4. Crosses穿越,一般是一条线穿过一个多边形?

poly2.Crosses(line)

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

4.4.5. Within包含,一个要素完全被另一个要素圈起来了?

ptB.Within(poly1)

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

4.4.6. Contains包含,跟Within正好相反?

poly1.Contains(ptB)

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

4.4.7. Overlaps重叠,好像只有两个多边形之间才能overlap?

poly2.Overlaps(poly3)

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

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

多边形的:

交:poly3.Intersection(poly2)

并:poly3.Union(poly2)

差:poly3.Difference(poly2)

补:poly3.SymmetricDifference(poly2)

geometry的:

<geom>.Buffer(<distance>) 给geometry加buffer,就是把点线变成多边形,变粗了

<geom1>.Equal(<geom2>) 两个geometry相等吗?

<geom1>.Distance(<geom2>) 返回两个geometry之间的最短距离

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

python的函数function,异常exception和模块module可以从任何一本python教材上找到,在此不多赘述。