利用OGR处理几何要素

OGR常量表示的不同几何类型 虽然许多几何要素只存在于二维(2D)的直角坐标系平面(X和Y坐标)，也有三维(3D)带有Z值得几何对象。这些Z值通常用来表示高程，但也可以用于表示其它数据，如每年温度得最大值。从技术上讲，这些几何要素在OGR中被认为是2.5D，而不是3D，因为OGR进行空间操作时，不考虑Z值。注意，虽然可以将Z值添加到二维几何结构中，但将数据写入文件时，Z值将被忽略。

from osgeo import ogr ogr.Geometry(ogr.wkbPoint) 几何类型2D常量2.5常量点wkbPointwkbPoint25D多点wkbMultiPointwkbMultiPoint25D线wkbLineStringwkbLineString25D多线wkbMultiLineStringwkbMultiLineString25D多边形环wkbLinearRingn/a面wkbPolygonwkbPolygon25D多面wkbMultiPolygonwkbMultiPolygon25D几何集合wkbGeometryCollectionwkbGeometryCollection25D

处理点集

创建并编辑一个点

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter firepit = ogr.Geometry(ogr.wkbPoint) firepit.AddPoint(59.5,11.5) print(firepit) x,y = firepit.GetX(), firepit.GetY() print(str(x)+','+str(y)) vp = VectorPlotter(True) vp.plot(firepit,'bo') POINT (59.5 11.5 0) 59.5,11.5

更改一个已经存在的点

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter firepit = ogr.Geometry(ogr.wkbPoint) firepit.AddPoint(59.5,11.5) print(firepit) x,y = firepit.GetX(), firepit.GetY() print(str(x)+','+str(y)) # 更改点坐标 firepit.AddPoint(59.5,13) print(firepit) x,y = firepit.GetX(), firepit.GetY() print(str(x)+','+str(y)) vp = VectorPlotter(True) vp.plot(firepit,'rs') POINT (59.5 11.5 0) 59.5,11.5 POINT (59.5 13.0 0) 59.5,13.0

创建一个2.5D的点

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter firepit = ogr.Geometry(ogr.wkbPoint25D) firepit.AddPoint(59.5,11.5,2) print(firepit) x,y,z = firepit.GetX(), firepit.GetY(), firepit.GetZ() print(str(x)+','+str(y) + ','+str(z)) vp = VectorPlotter(True) vp.plot(firepit,'bo') POINT (59.5 11.5 2) 59.5,11.5,2.0

多点作为一个几何类型

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter faucets = ogr.Geometry(ogr.wkbMultiPoint) faucet = ogr.Geometry(ogr.wkbPoint) faucet.AddPoint(67.5, 16) faucets.AddGeometry(faucet) # 添加到多点对象 faucet.AddPoint(73,31) faucets.AddGeometry(faucet) # 添加到多点对象 faucet.AddPoint(91,24.5) faucets.AddGeometry(faucet) # 添加到多点对象 #vp.clear() vp.plot(faucets,'bo') vp.zoom(-5) print(faucets) MULTIPOINT (67.5 16.0 0,73 31 0,91.0 24.5 0)

for i in range(faucets.GetGeometryCount()): pt = faucets.GetGeometryRef(i) pt.AddPoint(pt.GetX() + 2, pt.GetY()) vp.plot(faucets,'rs') vp.zoom(-5)

处理线要素

创建并编辑单条线

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter sidewalk = ogr.Geometry(ogr.wkbLineString) sidewalk.AddPoint(54, 37) sidewalk.AddPoint(62, 35.5) sidewalk.AddPoint(70.5, 38) sidewalk.AddPoint(74.5, 41.5) vp.plot(sidewalk,'b--') print(sidewalk) LINESTRING (54 37 0,62.0 35.5 0,70.5 38.0 0,74.5 41.5 0)

for i in range(sidewalk.GetPointCount()): sidewalk.SetPoint(i,sidewalk.GetX(i), sidewalk.GetY(i)+1) vp.plot(sidewalk,'r--')

print(sidewalk.GetPoints()) [(54.0, 38.0, 0.0), (62.0, 36.5, 0.0), (70.5, 39.0, 0.0), (74.5, 42.5, 0.0)] vertices = sidewalk.GetPoints() vertices[2:2] = [(66.5,35)] print(vertices) new_sidewalk = ogr.Geometry(ogr.wkbLineString) for vertex in vertices: new_sidewalk.AddPoint(*vertex) vp.plot(new_sidewalk,'g:') [(54.0, 38.0, 0.0), (62.0, 36.5, 0.0), (66.5, 35), (70.5, 39.0, 0.0), (74.5, 42.5, 0.0)]

import os data_dir = "D:\中国交通通信信息中心\书籍配套代码\python地理数据处理\资源\osgeopy-data\osgeopy-data\osgeopy-data-misc\osgeopy-data" ds = ogr.Open(os.path.join(data_dir, 'misc', 'line-example.geojson')) lyr = ds.GetLayer() feature = lyr.GetFeature(0) line = feature.geometry().Clone() vp.clear() vp.plot(line, 'b-')

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter import os data_dir = "D:\中国交通通信信息中心\书籍配套代码\python地理数据处理\资源\osgeopy-data\osgeopy-data\osgeopy-data-misc\osgeopy-data" ds = ogr.Open(os.path.join(data_dir, 'misc', 'line-example.geojson')) lyr = ds.GetLayer() feature = lyr.GetFeature(0) line = feature.geometry().Clone() vertices = line.GetPoints() vertices = line.GetPoints() vertices[26:26] = [(87, 57)] vertices[19:19] = [(95, 38), (97, 43), (101, 42)] vertices[11:11] = [(121, 18)] vertices[5:5] = [(67, 32), (74, 30)] new_line = ogr.Geometry(ogr.wkbLineString) for vertex in vertices: new_line.AddPoint(*vertex) vp.plot(new_line, 'b--')

线创建点

def line_to_point_layer(ds, line_name, pt_name): """Creates a point layer from vertices in a line layer.""" # Delete the point layer if it exists. if ds.GetLayer(pt_name): ds.DeleteLayer(pt_name) # Get the line layer and its spatial reference. line_lyr = ds.GetLayer(line_name) sr = line_lyr.GetSpatialRef() # Create a point layer with the same SR as the lines # and copy the field definitions from the line to # the point layer. pt_lyr = ds.CreateLayer(pt_name, sr, ogr.wkbPoint) pt_lyr.CreateFields(line_lyr.schema) # Create a feature and geometry to use over and over. pt_feat = ogr.Feature(pt_lyr.GetLayerDefn()) pt_geom = ogr.Geometry(ogr.wkbPoint) # Loop through all of the lines. for line_feat in line_lyr: # Copy the attribute values from the line to the # new feature. atts = line_feat.items() for fld_name in atts.keys(): pt_feat.SetField(fld_name, atts[fld_name]) # Loop through the line's vertices and for each one, # set the coordinates on the point geometry, add # it to the feature, and use the feature to create # a new feature in the point layer. for coords in line_feat.geometry().GetPoints(): pt_geom.AddPoint(*coords) pt_feat.SetGeometry(pt_geom) pt_lyr.CreateFeature(pt_feat) from osgeo import ogr from ospybook.vectorplotter import VectorPlotter import os import ospybook as pb data_dir = "D:\中国交通通信信息中心\书籍配套代码\python地理数据处理\资源\osgeopy-data\osgeopy-data\osgeopy-data-misc\osgeopy-data" ds = ogr.Open(os.path.join(data_dir, 'misc', 'line-example.geojson')) lyr = ds.GetLayer() feature = lyr.GetFeature(0) line = feature.geometry().Clone() pb.print_layers(os.path.join(data_dir, 'misc', 'line-example.geojson')) 0: line-example (LineString)

多线作为一个几何类型

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter path1 = ogr.Geometry(ogr.wkbLineString) path1.AddPoint(61.5, 29) path1.AddPoint(63, 20) path1.AddPoint(62.5, 16) path1.AddPoint(60, 13) path2 = ogr.Geometry(ogr.wkbLineString) path2.AddPoint(60.5, 12) path2.AddPoint(68.5, 13.5) path3 = ogr.Geometry(ogr.wkbLineString) path3.AddPoint(69.5, 33) path3.AddPoint(80, 33) path3.AddPoint(86.5, 22.5) paths = ogr.Geometry(ogr.wkbMultiLineString) paths.AddGeometry(path1) paths.AddGeometry(path2) paths.AddGeometry(path3) # Take a look at the multiline. vp.clear() vp.plot(paths, 'b-') print(paths) MULTILINESTRING ((61.5 29.0 0,63 20 0,62.5 16.0 0,60 13 0),(60.5 12.0 0,68.5 13.5 0),(69.5 33.0 0,80 33 0,86.5 22.5 0))

处理多边形

创建并编辑单多边形

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter ring = ogr.Geometry(ogr.wkbLinearRing) ring.AddPoint(58, 38.5) ring.AddPoint(53, 6) ring.AddPoint(99.5, 19) ring.AddPoint(73, 42) yard = ogr.Geometry(ogr.wkbPolygon) yard.AddGeometry(ring) yard.CloseRings() # Take a look at the polygon. Setting fill=False makes the polygon hollow when # it is drawn. vp = VectorPlotter(True) vp.plot(yard, fill=False, edgecolor='blue') print(yard) POLYGON ((58.0 38.5 0,53 6 0,99.5 19.0 0,73 42 0,58.0 38.5 0))

ring = yard.GetGeometryRef(0) for i in range(ring.GetPointCount()): ring.SetPoint(i, ring.GetX(i) - 5, ring.GetY(i)) vp.plot(yard, fill=False, ec='red', linestyle='dashed')

ring = yard.GetGeometryRef(0) vertices = ring.GetPoints() vertices[2:3] = ((90, 16), (90, 27)) for i in range(len(vertices)): ring.SetPoint(i, *vertices[i]) vp.plot(yard, fill=False, ec='black', ls='dotted', linewidth=3)

从多边形创建线

def poly_to_line_layer(ds, poly_name, line_name): """Creates a line layer from a polygon layer.""" # Delete the line layer if it exists. if ds.GetLayer(line_name): ds.DeleteLayer(line_name) # Get the polygon layer and its spatial reference. poly_lyr = ds.GetLayer(poly_name) sr = poly_lyr.GetSpatialRef() # Create a line layer with the same SR as the polygons # and copy the field definitions from the polygons to # the line layer. line_lyr = ds.CreateLayer(line_name, sr, ogr.wkbLineString) line_lyr.CreateFields(poly_lyr.schema) # Create a feature to use over and over. line_feat = ogr.Feature(line_lyr.GetLayerDefn()) # Loop through all of the polygons. for poly_feat in poly_lyr: # Copy the attribute values from the polygon to the # new feature. atts = poly_feat.items() for fld_name in atts.keys(): line_feat.SetField(fld_name, atts[fld_name]) # Loop through the rings in the polygon. poly_geom = poly_feat.geometry() for i in range(poly_geom.GetGeometryCount()): ring = poly_geom.GetGeometryRef(i) # Create a new line using the ring's vertices. line_geom = ogr.Geometry(ogr.wkbLineString) for coords in ring.GetPoints(): line_geom.AddPoint(*coords) # Insert the new line feature. line_feat.SetGeometry(line_geom) line_lyr.CreateFeature(line_feat)

复合多边形形成一个几何类型

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter box1 = ogr.Geometry(ogr.wkbLinearRing) box1.AddPoint(87.5, 25.5) box1.AddPoint(89, 25.5) box1.AddPoint(89, 24) box1.AddPoint(87.5, 24) garden1 = ogr.Geometry(ogr.wkbPolygon) garden1.AddGeometry(box1) box2 = ogr.Geometry(ogr.wkbLinearRing) box2.AddPoint(89, 23) box2.AddPoint(92, 23) box2.AddPoint(92,22) box2.AddPoint(89,22) garden2 = ogr.Geometry(ogr.wkbPolygon) garden2.AddGeometry(box2) gardens = ogr.Geometry(ogr.wkbMultiPolygon) gardens.AddGeometry(garden1) gardens.AddGeometry(garden2) gardens.CloseRings() # Take a look at the multipolygon. vp.clear() vp.plot(gardens, fill=False, ec='blue') vp.zoom(-1) print(gardens) MULTIPOLYGON (((87.5 25.5 0,89.0 25.5 0,89 24 0,87.5 24.0 0,87.5 25.5 0)),((89 23 0,92 23 0,92 22 0,89 22 0,89 23 0)))

创建并编辑带洞多边形：甜甜圈

from osgeo import ogr from ospybook.vectorplotter import VectorPlotter lot = ogr.Geometry(ogr.wkbLinearRing) lot.AddPoint(58, 38.5) lot.AddPoint(53, 6) lot.AddPoint(99.5, 19) lot.AddPoint(73, 42) house = ogr.Geometry(ogr.wkbLinearRing) house.AddPoint(67.5, 29) house.AddPoint(69, 25.5) house.AddPoint(64, 23) house.AddPoint(69, 15) house.AddPoint(82.5, 22) house.AddPoint(76, 31.5) yard = ogr.Geometry(ogr.wkbPolygon) yard.AddGeometry(lot) yard.AddGeometry(house) yard.CloseRings() # Take a look at the polygon. vp.clear() vp.plot(yard, 'yellow') print(yard) POLYGON ((58.0 38.5 0,53 6 0,99.5 19.0 0,73 42 0,58.0 38.5 0),(67.5 29.0 0,69.0 25.5 0,64 23 0,69 15 0,82.5 22.0 0,76.0 31.5 0,67.5 29.0 0))

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)) vp.plot(yard, fill=False, hatch='x', color='blue')