>>> from env_helper import info; info()
页面更新时间： 2024-07-23 22:54:16
运行环境：
Linux发行版本: Debian GNU/Linux 12 (bookworm)
操作系统内核: Linux-6.1.0-23-amd64-x86_64-with-glibc2.36
Python版本: 3.11.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)
Warning 1: Value 50654050.6944999993 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width
SetSpatialFilter(None)则是清空空间属性的过滤器，可查看输出图层要素的数目。
4.4.1. 在OGR中使用SQL语句进行查询¶
属性与空间的筛选可以看作是OGR的内置功能，这两种功能可以解决大部分实际应用问题。但是也有查询条件复杂的情况。针对这种情况，OGR也提供了灵活的解决方案：支持SQL查询语句。
例如执行SQL查询语句ExecuteSQL(<SQL>)，凭借SQL的强大功能，可执行更复杂的任务。例如下面这段代码，是从东北地区的分县数据中选择出吉林省的县级行政单位（对应的Prov_ID为22），并且按行政代码（）降序输出。
>>> 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)
Warning 1: Value 50654050.6944999993 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 50654050.6944999993 of field area of feature 1 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 29220969.7270000018 of field area of feature 2 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 20154740.0899999999 of field area of feature 3 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 17534164.0667999983 of field area of feature 4 not successfully written. Possibly due to too larger number with respect to field width
Warning 1: Value 13919140.5654000007 of field area of feature 152784 not successfully written. Possibly due to too larger number with respect to field width
上面使用的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()
以下是生成的结果：
图 4.2 image¶
这是一幅不完整的世界地图，因为它只包含了面积在80万平方公里以上的国家。对于这个结果，首先要注意的是有一些面积比较小的图斑也被选了出来，这些图斑是属于某些国家的，其面积定义是与某个大图斑一致的，也就是说，此处的面积是国家的面积，而并非是图斑的面积。另外要注意的是台湾并没有被选择出来：这是个政治问题，而非技术问题。