空间过滤器(Spatial filters)

如果说按照属性筛选要素是带有数据库特征的话,那么,根据空间位置的筛选就是纯GIS了。在OGR中,使用了Spatial filters(空间过滤)这一术语表征这一功能。

OGR提供的空间过滤功能有两种,一种是SetSpatialFilter(<geom>)---过滤某一类型的Feature,如参数中的Polygon,效用就是选出Layer中的所有Polygon覆盖的要素(注意,只要相交即可,不必完全包含)。

下面这段代码用了两套数据。 world_borders 是全球国界数据, cover.shp是覆盖了非洲南部地区的一个多边形数据。

首先定义一个根据图层直接生成shape文件的函数,方便后面调用,再使用cover.shp中的多边形来筛选全球国界数据:

image

In [1]:
from osgeo import ogr
def create_shp_by_layer(shp, layer):
    outputfile = shp
    if os.access(outputfile, os.F_OK):
        driver.DeleteDataSource(outputfile)
    newds = driver. CreateDataSource ( outputfile )
    pt_layer = newds.CopyLayer ( layer, '')
    newds.Destroy ()

下面代码是使用cover.shp中的多边形来选择全球国界数据:

In [10]:
driver = ogr.GetDriverByName("ESRI Shapefile")
world_shp = '/gdata/world_borders.shp'
cover_shp = '/gdata/10740.shp'
world_ds = ogr.Open(world_shp)
cover_ds = ogr.Open(cover_shp)
world_layer = world_ds.GetLayer(0)
cover_layer = cover_ds.GetLayer(0)
In [11]:
world_layer.GetFeatureCount()
Out[11]:
3784
In [12]:
cover_feats = cover_layer.GetNextFeature()
poly = cover_feats.GetGeometryRef()
world_layer.SetSpatialFilter(poly)
out_shp = '/tmp/world_cover.shp'
create_shp_by_layer(out_shp, world_layer)

另外还有 SetSpatialFilterRect(<minx>, <miny>, <maxx>, <maxy>) 参数,可输入四个坐标选中矩形内的要素。

In [8]:
from osgeo import ogr
import os
shpfile = '/gdata/world_borders.shp'
ds = ogr.Open(shpfile)
world_layer = ds.GetLayer(0)
world_layer.SetSpatialFilterRect(50, 60, 25, 35)
print(world_layer.GetFeatureCount())
68
In [9]:
out_shp = '/gdata/world_spatial_filter.shp'
create_shp_by_layer(out_shp, world_layer)

SetSpatialFilter(None)则是清空空间属性的过滤器,可查看输出图层要素的数目。

在OGR中使用SQL语句进行查询

属性与空间的筛选可以看作是OGR的内置功能,这两种功能可以解决大部分实际应用问题。但是也有查询条件复杂的情况。针对这种情况,OGR也提供了灵活的解决方案:支持SQL查询语句。

例如执行SQL查询语句ExecuteSQL(<SQL>),凭借SQL的强大功能,可执行更复杂的任务。例如下面这段代码,是从东北地区的分县数据中选择出吉林省的县级行政单位(对应的Prov_ID22),并且按行政代码()降序输出。

In [13]:
from osgeo import ogr
import os

def create_shp_by_layer(shp, layer):
    outputfile = shp
    if os.access(outputfile, os.F_OK):
        driver.DeleteDataSource(outputfile)
    newds = driver. CreateDataSource ( outputfile )
    pt_layer = newds.CopyLayer ( layer, '')
    newds.Destroy ()
    
In [14]:
driver = ogr.GetDriverByName("ESRI Shapefile")
world_shp = '/gdata/world_borders.shp'
world_ds = ogr.Open(world_shp)
world_layer = world_ds.GetLayer()
world_layer_name = world_layer.GetName()
result = world_ds.ExecuteSQL("select * from %s " % (world_layer_name)) # )
resultFeat = result.GetNextFeature ()
out_shp = '/tmp/sql_res.shp'
create_shp_by_layer(out_shp, result)
world_ds.ReleaseResultSet(result)

上面使用的SQL语句与平常的SQL语句无太大区别,相同点是使用了SELECT语句与WHERE条件,不同点是在OGR中此语句会生成空间数据。

最后一句ReleaseResultSet(<result_layer>)是将查询结果释放,在执行下一条SQL语句之前一定要先释放。为简便而言,可同样将查询的结果当成是数据来查看。

In [15]:
# while resultFeat :
#     print (resultFeat.GetField('name'))
#     resultFeat = result.GetNextFeature ()

注意:ExecuteSQL是基于数据集进行的,而不是图层。ExecuteSQL是简单的SQL字符串,不能用作任何解译。

应用

统计草地(grass)覆盖(cover)的所有要素(Feature)数目:

In [2]:
result = dsSites.ExecuteSQL("select count(*) from sites where cover = '/gdata/grass'")
result.GetFeatureCount()
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-2-4b8f6321019c> in <module>()
      2 import os
      3 shpfile = '/gdata/world_borders.shp'
----> 4 result = dsSites.ExecuteSQL("select count(*) from sites where cover = '/gdata/grass'")
      5 result.GetFeatureCount()

NameError: name 'dsSites' is not defined
In [3]:
result.GetFeature(0).GetField(0)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-3-3bbc56bfb97c> in <module>()
----> 1 result.GetFeature(0).GetField(0)

NameError: name 'result' is not defined
In [4]:
dsSites.ReleaseResultSet(result)
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-174466a212ed> in <module>()
----> 1 dsSites.ReleaseResultSet(result)

NameError: name 'dsSites' is not defined

列出所有cover类型:

In [ ]:
result = ds.ExecuteSQL("select distinct cover from sites")
resultFeat = result.GetNextFeature()
while resultFeat:
    print resultFeat.GetField(0)
    resultFeat = result.GetNextFeature()
ds.ReleaseResultSet(result)

统计各种cover类型分别有多少个要素(Feature):

In [5]:
coverLayer = ds.ExecuteSQL('/gdata/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)                  
  File "<ipython-input-5-fb2faab008e9>", line 4
    cntLayer = ds.ExecuteSQL("select count( * from sites where cover = ‘ “ + coverFeat.GetField(0) + “ ‘ “)
           ^
IndentationError: expected an indented block

根据属性条件选择要素

根据属性选择数据库记录是数据库应用中必不可少的一项功能。地理数据也是如此,譬如筛选人口在百万以上的城市,面积在百万平方千米以上的国家等等。

Layer对象中有 SetAttributeFilter(<where_clause>)方法,可根据属性将Layer中符合某一条件的Feature过滤出来。设定了Filter之后用GetNextFeature()方法依次筛选出符合条件的Feature:

In [19]:
from osgeo import ogr
import os
shpfile = '/gdata/world_borders.shp'
ds = ogr.Open(shpfile)
layer = ds.GetLayer(0)
lyr_count = layer.GetFeatureCount()
print(lyr_count)
3784
In [20]:
layer.SetAttributeFilter("AREA > 800000")
lyr_count = layer.GetFeatureCount()
print(lyr_count)
1790

先输出图层中的要素数目,再使用过滤器选择面积在80万平方公里以上的国家(这里使用了“AREA”字段,这个字段是早在Shapefile中生成好的),因此输出图层要素数目时,根据筛选后的结果,数据量随之变少。

根据属性条件生成要素

然而,筛选要素的目的并不是为了简单地查看一下数目,而是要将此数据保存,以便后期使用。请仔细查看以下代码。注意,可在不同的层次上生成矢量数据。

In [21]:
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = '/tmp/xx_filter_demo.shp'
if os.access( extfile, os.F_OK ):
    driver.DeleteDataSource( extfile )
In [22]:
newds = driver.CreateDataSource(extfile)
layernew = newds.CreateLayer('rect',None,ogr.wkbPolygon)
feat = layer.GetNextFeature()
while feat is not None:
    layernew.CreateFeature(feat)
    feat = layer.GetNextFeature()
newds.Destroy()

以下是生成的结果:

image

这是一幅不完整的世界地图,因为它只包含了面积在80万平方公里以上的国家。对于这个结果,首先要注意的是有一些面积比较小的图斑也被选了出来,这些图斑是属于某些国家的,其面积定义是与某个大图斑一致的,也就是说,此处的面积是国家的面积,而并非是图斑的面积。另外要注意的是台湾并没有被选择出来:这是个政治问题,而非技术问题。