>>> from env_helper import info; info()
页面更新时间: 2022-12-28 08:12:05
运行环境:
    Linux发行版本: Debian GNU/Linux 11 (bullseye)
    操作系统内核: Linux-5.10.0-20-amd64-x86_64-with-glibc2.31
    Python版本: 3.9.2

4.4. 空间过滤器(Spatial filters)

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

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

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

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

>>> import os
>>> 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中的多边形来选择全球国界数据:

>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> world_shp = '/gdata/GSHHS_h.shp'
>>> cover_shp = '/gdata/spatial_filter.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)
>>> world_layer.GetFeatureCount()
153113
>>> 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>) 参数,可输入四个坐标选中矩形内的要素。

>>> world_layer.SetSpatialFilter(None)
>>> world_layer.GetFeatureCount()
153113
>>> world_layer.SetSpatialFilterRect(50, 60, 25, 35)
>>> print(world_layer.GetFeatureCount())
1806
>>> out_shp = '/tmp/world_spatial_filter.shp'
>>> create_shp_by_layer(out_shp, world_layer)

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

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

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

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

>>> 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 ()
>>>
>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> world_shp = '/gdata/GSHHS_h.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语句之前一定要先释放。为简便而言,可同样将查询的结果当成是数据来查看。

4.4.2. 根据属性条件选择要素

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

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

>>> from osgeo import ogr
>>> import os
>>> shpfile = '/gdata/GSHHS_h.shp'
>>> ds = ogr.Open(shpfile)
>>> layer = ds.GetLayer(0)
>>> lyr_count = layer.GetFeatureCount()
>>> print(lyr_count)
153113
>>> layer.SetAttributeFilter("AREA > 800000")
>>> lyr_count = layer.GetFeatureCount()
>>> print(lyr_count)
8

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

根据属性条件生成要素

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

>>> driver = ogr.GetDriverByName("ESRI Shapefile")
>>> extfile = '/tmp/xx_filter_demo.shp'
>>> if os.access( extfile, os.F_OK ):
>>>     driver.DeleteDataSource( extfile )
>>> 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

图 4.2 image

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