如果说按照属性筛选要素是带有数据库特征的话，那么，根据空间位置的筛选就是典型的GIS操作了。
在OGR中，使用了Spatial filters（空间过滤）这一术语表征这一功能。
OGR提供的空间过滤功能有两种，一种是 SetSpatialFilter(geom) ---过滤某一类型的Feature，如参数中的Polygon，效用就是选出Layer中的所有Polygon覆盖的要素（注意，只要相交即可，不必完全包含）。
空间择舍器
下面这段代码用了两套数据。 world_borders 是全球国界数据， spatial_filter是覆盖了非洲南部地区的一个多边形数据。
首先定义一个根据图层直接生成shape文件的函数，方便后面调用，再使用cover.shp中的多边形来筛选全球国界数据：
In [2]:
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中的多边形来选择全球国界数据：
In [3]:
driver = ogr.GetDriverByName("ESRI Shapefile")
world_shp = '/data/gdata/GSHHS_h.shp'
cover_shp = '/data/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)
/opt/conda/lib/python3.12/site-packages/osgeo/ogr.py:601: FutureWarning: Neither ogr.UseExceptions() nor ogr.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default. warnings.warn(
In [4]:
world_layer.GetFeatureCount()
Out[4]:
153113
In [5]:
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 [6]:
world_layer.SetSpatialFilter(None)
world_layer.GetFeatureCount()
Out[6]:
153113
In [7]:
world_layer.SetSpatialFilterRect(50, 60, 25, 35)
print(world_layer.GetFeatureCount())
1806
In [8]:
out_shp = '/tmp/world_spatial_filter.shp'
create_shp_by_layer(out_shp, world_layer)
Warning 1: Value 50654050.694499999 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width
SetSpatialFilter(None)则是清空空间属性的过滤器，可查看输出图层要素的数目。
In [9]:
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 [10]:
driver = ogr.GetDriverByName("ESRI Shapefile")
world_shp = '/data/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.694499999 of field area of feature 0 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 50654050.694499999 of field area of feature 1 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 29220969.727000002 of field area of feature 2 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 20154740.09 of field area of feature 3 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 17534164.066799998 of field area of feature 4 not successfully written. Possibly due to too larger number with respect to field width Warning 1: Value 13919140.565400001 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语句之前一定要先释放。为简便而言，可同样将查询的结果当成是数据来查看。
In [11]:
from osgeo import ogr
import os
shpfile = '/data/gdata/GSHHS_h.shp'
ds = ogr.Open(shpfile)
layer = ds.GetLayer(0)
lyr_count = layer.GetFeatureCount()
print(lyr_count)
153113
In [12]:
layer.SetAttributeFilter("AREA > 800000")
lyr_count = layer.GetFeatureCount()
print(lyr_count)
8
In [13]:
driver = ogr.GetDriverByName("ESRI Shapefile")
extfile = '/tmp/xx_filter_demo.shp'
if os.access( extfile, os.F_OK ):
driver.DeleteDataSource( extfile )
In [14]:
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()
以下是生成的结果：
这是一幅不完整的世界地图，因为它只包含了面积在80万平方公里以上的国家。对于这个结果，首先要注意的是有一些面积比较小的图斑也被选了出来，这些图斑是属于某些国家的，其面积定义是与某个大图斑一致的，也就是说，此处的面积是国家的面积，而并非是图斑的面积。另外要注意的是台湾并没有被选择出来：这是个政治问题，而非技术问题。