目录

上一个主题

3.3. 使用OGR 创建Shapefile

下一个主题

3.5. 空间计算

关注公众号


常见问题

  1. Windows下的安装说明
  2. Jupyter免费在线实验环境
  3. 勘误与补充


>>> from helper import info; info()
待更新

3.4. 空间过滤器(Spatial filters)

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

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

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

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

image

图 3.2 image

>>> 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_c.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)
---------------------------------------------------------------------------

AttributeError                            Traceback (most recent call last)

<ipython-input-3-bb3c50025982> in <module>()
      5 cover_ds = ogr.Open(cover_shp)
      6 world_layer = world_ds.GetLayer(0)
----> 7 cover_layer = cover_ds.GetLayer(0)


AttributeError: 'NoneType' object has no attribute 'GetLayer'
>>> world_layer.GetFeatureCount()
>>> 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>) 参数,可输入四个坐标选中矩形内的要素。

>>> 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())
>>> out_shp = '/gdata/world_spatial_filter.shp'
>>> create_shp_by_layer(out_shp, world_layer)

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

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

image0

>>>  while resultFeat :
>>>      print (resultFeat.GetField('name'))
>>>      resultFeat = result.GetNextFeature ()

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

3.4.2. 应用

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

>>> result = dsSites.ExecuteSQL("select count(*) from sites where cover = '/gdata/grass'")
>>> result.GetFeatureCount()
>>> result.GetFeature(0).GetField(0)
>>> dsSites.ReleaseResultSet(result)

列出所有cover类型:

>>> 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):

>>> 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)

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

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

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

>>> 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)
>>> layer.SetAttributeFilter("AREA > 800000")
>>> lyr_count = layer.GetFeatureCount()
>>> print(lyr_count)

先输出图层中的要素数目,再使用过滤器选择面积在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

图 3.3 image

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