【GDAL学习】过滤器,简单的空间分析,函数和模块

1.属性过滤器

>>>import ogr,os
>>>os.chdir('E:/data/GDAL/ospy_data3')
>>>driver=ogr.GetDriverByName('ESRI Shapefile')
>>>ds=driver.Open('sites.shp',0)
>>>ly=ds.GetLayer()
>>>ly.GetFeatureCount()
Out[7]: 42
>>>ly.SetAttributeFilter("cover='shrubs'")
Out[8]: 0
>>>ly.GetFeatureCount()
Out[9]: 6
>>>ly.SetAttributeFilter(None)
Out[10]: 0
>>>ly.GetFeatureCount()
Out[11]: 42

2.空间过滤器

>>>import ogr
>>>import os
>>>os.chdir('E:/data/GDAL/ospy_data3')
>>>driver = ogr.GetDriverByName('ESRI Shapefile')
>>>ds1 = driver.Open('cache_towns.shp', 0)
>>>ds2 = driver.Open('sites.shp', 0)
>>>ly1 = ds1.GetLayer()
>>>ly2 = ds2.GetLayer()
>>>fAreas=ly1.GetNextFeature()
>>>poly=fAreas.GetGeometryRef()
>>>ly2.GetFeatureCount()
Out[12]: 42
>>>ly2.SetSpatialFilter(poly)
>>>ly2.GetFeatureCount()
Out[14]: 0
>>>ly2.SetSpatialFilterRect(460000,4590000,490000,4600000)
>>>ly2.GetFeatureCount()
Out[16]: 4
>>>ly2.SetSpatialFilter(None)
>>>ly2.GetFeatureCount()
Out[18]: 42

更强大的查询功能,执行SQL代码,查找出cover类型为grass的Feature并以ID降序排列。

最后一句ReleaseResultSet(<result_layer>)是将查询结果释放,在执行下一条SQL语句之前一定要先释放。

>>>import ogr,os
...os.chdir('E:/data/GDAL/ospy_data3')
...driver=ogr.GetDriverByName('ESRI Shapefile')
...ds=driver.Open('sites.shp',0)
...ly=ds.GetLayer()
...ly.GetFeatureCount()
Out[2]: 42
>>>result=ds.ExecuteSQL("select * from sites where cover = 'grass' order by id desc")
>>>resultFeat=result.GetNextFeature()
>>>while resultFeat:
    ...print(resultFeat.GetField('id'))
    ...resultFeat=result.GetNextFeature()
    
42
40
37
31
28
24
17
13
11
10
4
>>>ds.ReleaseResultSet(result)

 接上面代码。统计cover为grass的所有Feature的数目

>>>result=ds.ExecuteSQL("select count(*) from sites where cover = 'grass'")
>>>result.GetFeatureCount()
Out[8]: 1
>>>result.GetFeature(0).GetField(0)
Out[9]: 11
>>>ds.ReleaseResultSet(result)

 接上面代码,列出所有的cover类型

>>>result=ds.ExecuteSQL("select distinct cover from sites")
>>>resultFeat=result.GetNextFeature()
>>>while resultFeat:
    ...print(resultFeat.GetField(0))
    ...resultFeat = result.GetNextFeature()
    
shrubs
trees
rocks
grass
bare
water
>>>ds.ReleaseResultSet(result)

接上面代码,统计每种cover类型各有多少个Feature

>>>coverLayer = ds.ExecuteSQL('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)+' '+cntLayer.GetFeature(0).GetFieldAsString(0))
    ...ds.ReleaseResultSet(cntLayer)
    ...coverFeat = coverLayer.GetNextFeature()
    ...
shrubs 6
trees 11
rocks 6
grass 11
bare 6
water 2
>>>ds.ReleaseResultSet(coverLayer)

 3.空间操作

1)Intersect:判断两个要素是否相交,Polygon与Polygon,Polygon与Point,Polygon与Line,Line与Point

>>> poly2.Intersect(poly1)
0
>>> poly2.Intersect(poly3)
1
>>> poly2.Intersect(poly2)
1
>>> poly1.Intersect(ptA)
0
>>> poly1.Intersect(ptB)
1
>>> poly1.Intersect(line)
1
>>> poly3.Intersect(line)
1
>>> line.Intersect(ptB)
1

2)Disjoint:判断两个要素是否不相交,Polygon与Polygon,Polygon与Point,Polygon与Line,Line与Point

>>> poly2.Disjoint(poly1)
1
>>> poly2.Disjoint(poly3)
0
>>> poly1.Disjoint(ptA)
1
>>> poly1.Disjoint(ptB)
0
>>> poly1.Disjoint(line)
0
>>> poly3.Disjoint(line)
0
>>> line.Disjoint(ptB)

3)Touches:表示相邻(擦边),Polygon与Polygon,Polygon与Point,Polygon与Line

>>> poly2.Touches(poly1)
0
>>> poly2.Touches(poly3)
0
>>> poly1.Touches(line)
0
>>> poly1.Touches(ptB)
0
>>> poly3.Touches(line)
1

4)Crosses:穿越,一般是一条线穿过一个多边形,一般是Polygon与Line

>>> poly2.Crosses(poly1)
0
>>> poly2.Crosses(poly3)
0
>>> poly2.Crosses(line)
1
>>> poly3.Crosses(line)
0
>>> poly1.Crosses(line)
1
>>> line.Crosses(ptB)
0

5)Within:包含于,一个要素完全被另一个要素圈起来了,Polygon与Polygon,Polygon与Point,Polygon与Line

>>> poly3.Within(poly2)
0
>>> line.Within(poly2)
0
>>> ptA.Within(poly1)
0
>>> ptB.Within(poly1)
1
>>> poly1.Within(ptB)
0

6)Contains:包含,跟Within正好相反,一个要素完全包含另一个要素,Polygon与Polygon,Polygon与Point,Polygon与Line

>>> poly3.Contains(poly2)
0
>>> line.Contains(poly2)
0
>>> poly2.Contains(line)
0
>>> poly1.Contains(ptA)
0
>>> poly1.Contains(ptB)
1
>>> ptB.Contains(poly1)
0

7)Overlaps:重叠,只有两个多边形之间才能overlap,Polygon与Polygon

>>> poly2.Overlaps(poly1)
0
>>> poly2.Overlaps(poly3)
1
>>> poly2.Overlaps(line)
0
>>> poly3.Overlaps(line)
0
>>> poly1.Overlaps(ptB)
0

用week3中的数据测试一下,下图中多边形的FID是18,在cache_towns.shp文件中;点的FID是40,cover类型是tree,在sites.shp文件中。

 

>>>import ogr
...import os
...os.chdir('E:/data/GDAL/ospy_data3')
...driver = ogr.GetDriverByName('ESRI Shapefile')
...ds1 = driver.Open('sites.shp', 0)
...if ds1 is None:
...    print('Cannot open this file:', ds1.name)
...pyLy=ds1.GetLayer()
...treeFe=pyLy.GetFeature(40)
...treePt=treeFe.GetGeometryRef()
...ds2=driver.Open('cache_towns.shp',0)
...if ds2 is None:
...    print('Cannot open this file:',ds2.name)
...tnLy=ds2.GetLayer()
...townFe=tnLy.GetFeature(18)
...townPg=townFe.GetGeometryRef()
>>>townPg.Intersect(treePt)
Out[3]: True
>>>treePt.Intersect(townPg)
Out[4]: True
>>>townPg.Disjoint(treePt)
Out[5]: False
>>>townPg.Touches(treePt)
Out[6]: False
>>>townPg.Crosses(treePt)
Out[7]: False
>>>townPg.Within(treePt)
Out[8]: False
>>>townPg.Contains(treePt)
Out[9]: True

Assignment 3a:Use filters and buffers 

代码1:自己写的,没有注释

  • Use an attribute filter to restrict cache_towns.shp to Nibley (“name” field) 
  • Buffer the Nibley geometry by 1500 
  • Use the new geometry with a spatial filter on sites.shp to find all sites within 1500 meters of Nibley
  • Print out the “id” value for those sites
# Assignment 3a
import ogr
import os

os.chdir('E:/data/GDAL/ospy_data3')
driver = ogr.GetDriverByName('ESRI Shapefile')
ds1 = driver.Open('cache_towns.shp', 0)
if ds1 is None:
    print('Cannot open this file:', ds1.name)
ly1 = ds1.GetLayer()
ly1.SetAttributeFilter("NAME='Nibley'")
NibleyFeat = ly1.GetFeature(0)
NibleyGeom = NibleyFeat.GetGeometryRef()
bufferGeom = NibleyGeom.Buffer(1500)

ds2 = driver.Open('sites.shp', 0)
if ds2 is None:
    print('Cannot open this file:', ds2.name)
ly2 = ds2.GetLayer()
ly2.SetSpatialFilter(bufferGeom)
# ly2.GetFeatureCount()
siteFeat = ly2.GetNextFeature()
while siteFeat:
    print(siteFeat.GetField('ID'))
    siteFeat.Destroy()
    siteFeat = ly2.GetNextFeature()

ds2.Destroy()
ds1.Destroy()

 代码2:参考答案,有详细注释

# Assignment 3a
# import modules
import ogr, os, sys

# set the working directory
os.chdir('E:/data/GDAL/ospy_data3')

# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# open sites.shp and get the layer
sitesDS = driver.Open('sites.shp', 0)
if sitesDS is None:
    print('Could not open sites.shp')
    sys.exit(1)
sitesLayer = sitesDS.GetLayer()

# open cache_towns.shp and get the layer
townsDS = driver.Open('cache_towns.shp', 0)
if townsDS is None:
    print('Could not open cache_towns.shp')
    sys.exit(1)
townsLayer = townsDS.GetLayer()

# use an attribute filter to restrict cache_towns.shp to "Nibley"
townsLayer.SetAttributeFilter("NAME = 'Nibley'")

# get the Nibley geometry and buffer it by 1500
nibleyFeature = townsLayer.GetFeature(0)
nibleyGeom = nibleyFeature.GetGeometryRef()
bufferGeom = nibleyGeom.Buffer(1500)

# use bufferGeom as a spatial filter on sites.shp to get all sites
# within 1500 meters of Nibley
sitesLayer.SetSpatialFilter(bufferGeom)

# loop through the remaining features in sites.shp and print their
# id values
siteFeature = sitesLayer.GetNextFeature()
while siteFeature:
    print(siteFeature.GetField('ID'))
    siteFeature.Destroy()
    siteFeature = sitesLayer.GetNextFeature()

# close the shapefiles
sitesDS.Destroy()
townsDS.Destroy()

数据下载:https://download.csdn.net/download/qq_37935516/10791607

 

转载自:https://blog.csdn.net/qq_37935516/article/details/84192289

You may also like...