GDAL python教程基础篇（6）OGR空间滤波器

皮卡兵快跑

已于 2023-03-16 22:04:56 修改

阅读量627

点赞数

分类专栏： Python与开源GIS（GDAL、OSR、OGR）文章标签： python arcgis 开发语言

于 2023-03-15 23:05:16 首次发布

本文链接：https://blog.csdn.net/qq_45373844/article/details/129569525

版权

Python与开源GIS（GDAL、OSR、OGR）专栏收录该内容

14 篇文章 35 订阅

订阅专栏

1.空间过滤器

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

OGR提供的空间过滤功能有两种，一种是SetSpatialFilter(geom)—过滤某一类型的Feature，如参数中的Polygon，效用就是选出Layer中的所有Polygon覆盖的要素（注意，只要相交即可，不必完全包含）。另外还有 SetSpatialFilterRect(<minx>, <miny>, <maxx>, <maxy>) 参数，可输入四个坐标选中矩形内的要素。

1.1SetSpatialFilter

下面我们来做一个示例，首先这里有两份数据分别为全国的市级行政区划矢量文件和江苏省的行政区划。然后我们通过空间滤波器选择江苏省的下辖市，并将其写入单独的文件当中。
在这里插入图片描述
示例代码

from osgeo import ogr

provice_path = r'D:/work/全国行政区划/江苏省.shp'
city_path = r'D:/work/全国行政区划/84坐标行政区划/市84坐标.shp'

driver = ogr.GetDriverByName('ESRI Shapefile')    # 驱动

provice_ds = driver.Open(provice_path)    # 获取数据集
city_ds = driver.Open(city_path)

provice_layer = provice_ds.GetLayer(0)    # 获取图层
city_layer = city_ds.GetLayer(0)

print(city_layer.GetFeatureCount())   # 打印出全国一共有371个地级市

provice_feat = provice_layer.GetNextFeature()    # 获取要素

poly = provice_feat.GetGeometryRef()       # 获取几何
city_layer.SetSpatialFilter(poly)           # 空间滤波器获取包含在几何内的city

import os

out_shp = 'D:/work/全国行政区划/江苏省_市级行政区划.shp'       # 定义输出shp
if os.access(out_shp, os.F_OK):           # 判断是否已经存在该文件
    driver.DeleteDataSource(out_shp)          # 若存在该文件则删除
newds = driver. CreateDataSource (out_shp)      # 创建数据集
pt_layer = newds.CopyLayer (city_layer, '')        # 通过图层拷贝方法获取图层
newds.Destroy ()                                              # 销毁数据源，若不销毁无法写入磁盘

将生成的文件加载入ArcGIS Pro，其中黄色的图层即我们通过空间滤波器得到的江苏省的市级行政区划。这里由于
在这里插入图片描述

1.2 SetSpatialFilterRect

SetSpatialFilterRect(, , , )` 参数，可输入四个坐标选中矩形内的要素。

print(city_layer.GetFeatureCount())
>>>27
city_layer.SetSpatialFilter(poly)    #SetSpatialFilter(None)则是清空空间属性的过滤器，可查看输出图层要素的数目。
print(city_layer.GetFeatureCount())
>>>371
city_layer.SetSpatialFilterRect(118.055315, 29.877947, 121.431257, 31.4785669)
print(city_layer.GetFeatureCount())
>>>16

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

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

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

>>> from osgeo import ogr
>>> import os
>>> shpfile =  r'D:/work/全国行政区划/84坐标行政区划/市84坐标.shp'
>>> ds = ogr.Open(shpfile)
>>> layer = ds.GetLayer(0)
>>> lyr_count = layer.GetFeatureCount()
>>> print(lyr_count)
371
>>> layer.SetAttributeFilter("FID = 0")
>>> lyr_count = layer.GetFeatureCount()
>>> print(lyr_count)
1

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

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

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

>>> from osgeo import ogr
>>> import os
>>> 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语句之前一定要先释放。为简便而言，可同样将查询的结果当成是数据来查看。