GDAL笔记-chapter6
这一章是利用OGR处理几何要素。
首先介绍一下OGR常量表示的不同几何类型。其中注意线是wkbLineString\多边形环是wkbLinearRing，特别注意多边形环用于多边形中。
点 wkbPoint
多点 wkbMultiPoint
线 wkbLineString
多线 wkbMultiLineString
多边形环 wkbLinearRing
面 wkbPolygon
多面 wkbMultiPolygon
几何集合 wkbGeometryCollection
1.点
from osgeo import ogr firepit = ogr.Geometry(ogr.wkbPoint) #create a null geometry firepit.AddPoint(59.5, 11.5) #add a point firepit.SetPoint(59.5, 10) #rewrite a point
2.多点
###创建多点集合类型 faucets = ogr.Geometry(ogr.wkbMultiPoint) faucet = ogr.Geometry(ogr.wkbPoint) faucet.AddPoint(67.5, 16) faucets.AddGeometry(faucet) faucet.AddPonit(12, 34) #重用几何对象，单点会覆盖前点 faucets.AddGeometry(faucet) #重用几何对象 faucet.AddPoint(0, 12) faucets.AddGeometry(faucet)
####GetGeometryRef获取点的索引进行更新 faucets.GetGeometryRef(1).AddPoint(75, 32) #获取多点几何对象的索引，并且更新坐标 for i in range(faucets.GetGeometryCount()): #遍历所有点 pt = faucets.GetGeometryRef(i) pt.AddPoint(pt.GetX()+2, pt.GetY())
3.线
####创建线 sidewalk = ogr.Geometry(ogr.wkbLineString) sidewalk.AddPoint(54, 12) sidewalk.AddPoint(56, 23) sidewalk.AddPoint(12, 23) sidewalk.SetPoint(2, 12, 12) #改变索引为2的点的坐标
#GetGeometryCount确定几何形状集合中包含的几何图形数目 for i in range(sidewalk.GetPointCount()): sidewalk.SetPoint(i, sidewalk.GetX(i), sidewalk.GetY(i)+2) #GetPointCount方法返回单一几何类型的顶点个数
4.由后向前插入线的顶点可以不用更改索引号
#插入点时候，由后向前插入不用更改索引 vertices = line.GetPoints() vertices[26:26] = [(12, 23)] vertices[19:19] = [(11, 13),(14, 56), (12, 97)] vertices[11:11] = [(28, 23)] vertices[5:5] = [(12, 53)] new_line = ogr.Geometry(ogr.wkbLineString) for vertex in vertices: new_line.AddPoint(*vertex) #不创建副本直接新增点 vertices = sidewalk.GetPoints() vertices[2:2] = [(12, 45)] for i in range(len(vertices)): sidewalk.SetPoint(i, *vertices[i]) #原来的线缺少一个新点，SetPoint方法在请求的索引处新增了一个点，初始化为(0，0)，然后再更新
5.从线图层中创建点图层
检查新图层是否存在--获取线图层--获取线图层空间参考并且创建点图层--通过schema获取属性列表，创建点图层属性字段--获取图层定义并创建点图层要素--创建点图层几何体--遍历线图层要素，获取键值对元组--
遍历线要素的点--点要素添加几何体--点图层添加要素
#从线创建点 def line_to_point_layer(ds, line_name, pt_name): if ds.GetLayer(pt_name): ds.DeleteLayer(pt_name) line_lyr = ds.GetLayer(line_name) #获取线图层 sr = line_lyr.GetSpatialRef() #获取线图层的空间参考 pt_lyr = ds.CreateLayer(pt_name, sr, ogr.wkbPoint) #创建点图层 pt_lyr.CreateFields(line_lyr.schema) #创建点图层的属性字段，lne_lyr.schema获取属性列表 pt_feat = ogr.Feature(pt_lyr.GetLayerDefn()) #获取图层定义，创建点图层要素 pt_geom = ogr.Geometry(ogr.wkbPoint) #创建点图层几何体 for line_feat in line_lyr: #遍历线图层的每个要素 atts = line_feat.items() #获取线要素的键值对元组 for fld_name in atts.keys(): #遍历属性值 pt_feat.SetField(fld_name, atts[fld_name]) #复制属性值 for coords in line_feat.geometry().GetPoints(): #遍历线要素的点 pt_geom.AddPoint(*coords) pt_feat.SetGeometry(pt_geom) #点图层添加几何体 pt_lyr.CreateFeature(pt_feat) #点图层添加要素
6.多线
#创建多线 path1 = ogr.Geometry(ogr.wkbLineString) path1.AddPoint( 0, 0 ) path1.AddPoint( 0, 0 ) path1.AddPoint( 0, 0) path2 = ogr.Geometry(ogr.wkbLineString) path2.AddPoint( 0, 0) path2.AddPoint( 0, 0) path2.AddPoint( 0, 0) paths = ogr.Geometry(ogr.wkbMultiLineString) paths.AddGeometry(path1) paths.AddGeometry(path2) paths.GetGeometryRef(0).SetPoint(2, 63, 22) #GetGeometry获取多线中的第1条线，SetPoint更新第三个点 paths.GetGeometryRef() #获取多线中的线 #移动多线的每个点 for i in range(paths.GetGeometryCount()): #首先遍历多线的单线 path = paths.GetGeometryRef(i) for j in range(path.GetPointCount()): #然后遍历单线的每个点 path.SetPoint(j, path.GetX(j)+2, path.GetY(j)-3)
7.多边形（由环构成，最后需要CloseRings闭合）
#创建多边形 ring = ogr.Geometry(ogr.wkbLinearRing) #先创建环 ring.AddPoint(56, 57) ring.AddPoint(52, 51) ring.AddPoint(52, 53) ring.AddPoint(56, 53) yard = ogr.Geometry(ogr.wkbPolygon) #创建多边形几何体并且添加环 yard.AddGeometry(ring) yard.CloseRings() #闭合多边形 #多边形修改顶点 ring = yard.GetGeometryRef(0) for i in range(ring.GetPointCount()): ring.SetPoint(i, ring.GetX(i)-5, ring.GetY(i)) ring = yard.GetGeometryRef(0) #到环 vertices = ring.GetPoints() vertices[2:3] = ((90, 26),(90, 16)) for i in range(len(vertices)): ring.SetPoint(i, *vertices[i])
8.从多边形图层创建线图层的函数
#从多边形图层创建线图层的函数 def poly_to_line_layer(ds, poly_name, line_name): if ds.GetLayer(line_name): ds.DeleteLayer(line_name) poly_lyr = ds.GetLayer(poly_name) sr = poly_lyr.GetSpatialRef() line_lyr = ds.CreateLayer(line_name, sr, ogr.wkbLineString) #线图层，不是环 line_lyr.CreateFields(poly_lyr.schema) line_feat = ogr.Feature(line_lyr.GetLayerDefn()) for poly_feat in poly_lyr: atts = poly_feat.items() for fld_name in atts.key(): line_feat.SetField(fld_name, atts[fld_name]) poly_geom = poly_feat.GetGeometry() for i in range(poly_geom.GetGeometryCount()): ring = poly_geom.GetGeometryRef(i) line_geom = ogr.Geometry(ogr.wkbLineString) for coords in ring.GetPoints(): line_geom.AddPoint(*coords) line_feat.SetGeometry(line_geom) line_lyr.CreateFeature(line_feat)
9.复合多边形
#复合多边形 box1 = ogr.Geometry(ogr.wkbLinearRing) #环 box1.AddPoint(0, 0) box1.AddPoint(0, 0) box1.AddPoint(0, 0) garden1 = ogr.Geometry(ogr.wkbPolygon) garden1.AddGeometry(box1) box2 = ogr.Geometry(ogr.wkbLinearRing) box2.AddPoint(0, 0) box2.AddPoint(0, 0) box2.AddPoint(0, 0) garden2 = ogr.Geometry(ogr.wkbPolygon) garden2.AddGeometry(box2) gardens = ogr.Geometry(ogr.wkbMultiPolygon) #创建复合多边形几何体并添加多边形 gardens.AddGeometry(garden1) gardens.AddGeometry(garden2) gardens.CloseRings() #闭合多边形 #更新复合多边形顶点 for i in range(gardens.GetGeometryCount()): ring = gardens.GetGeometryRef(i).GetGeometryRef(0) #第一个GetGeometryRef从复合多边形中获取多边形，第二个GetGeometryRef获取多边形的环 for j in range(ring.GetPoints()): ring.SetPoint(j, ring.GetX(j)+1, ring.GetY(j)+2)
10.带洞的多边形
#创建带洞的多边形 lot = ogr.Geometry(ogr.wkbLinearRing) #环 lot.AddPoint(0, 0) lot.AddPoint(0, 0) lot.AddPoint(0, 0) house = ogr.Geometry(ogr.wkbLinearRing) house.AddPoint(0, 0) house.AddPoint(0, 0) house.AddPoint(0, 0) yard = ogr.Geometry(ogr.wkbPolygon) yard.AddGeometry(lot) yard.AddGeometry(house) yard.CloseRings() #更新带洞多边形的点 for i in range(yard.GetGeometryCount()): #先遍历环 ring = yard.GetGeometryRef(i) for j in range(ring.GetPointCount): #取点 ring.SetPoint(j, ring.GetX(j)-5, ring.GetY(j))
