GDAL python教程（3）——过滤器，简单的空间分析，函数和模块
===============================================================


属性过滤器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

空间过滤器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

统计每种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的：

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

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

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

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

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