# GDAL的几何操作

Date: 2015-08-13 12:43
Summary: 一段代码一个小功能，简单清晰又实用。翻译自英文的cookbook。Thanks the author for sharing us such a wonderful tutorial！！！

# GDAL的几何操作

## 创建一个点

``````from osgeo import ogr
point = ogr.Geometry(ogr.wkbPoint)
print point.ExportToWkt()
``````

## 创建一条线

``````from osgeo import ogr
line = ogr.Geometry(ogr.wkbLineString)
print line.ExportToWkt()
``````

## 创建多边形

``````from osgeo import ogr

# Create ring
ring = ogr.Geometry(ogr.wkbLinearRing)

# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)

print poly.ExportToWkt()
``````

## 创建带洞的多边形

``````from osgeo import ogr

# Create outer ring 外环线
outRing = ogr.Geometry(ogr.wkbLinearRing)

# Create inner ring 内环线
innerRing = ogr.Geometry(ogr.wkbLinearRing)

# Create polygon
poly = ogr.Geometry(ogr.wkbPolygon)

print poly.ExportToWkt()
``````

## 创建MultiPoint

``````from osgeo import ogr

multipoint = ogr.Geometry(ogr.wkbMultiPoint)

point1 = ogr.Geometry(ogr.wkbPoint)

point2 = ogr.Geometry(ogr.wkbPoint)

point3 = ogr.Geometry(ogr.wkbPoint)

print multipoint.ExportToWkt()
``````

## 创建MultiLineString

``````from osgeo import ogr

multiline = ogr.Geometry(ogr.wkbMultiLineString)

line1 = ogr.Geometry(ogr.wkbLineString)

line1 = ogr.Geometry(ogr.wkbLineString)

print multiline.ExportToWkt()
``````

## 创建MultiPolygon

``````from osgeo import ogr

multipolygon = ogr.Geometry(ogr.wkbMultiPolygon)

# Create ring #1
ring1 = ogr.Geometry(ogr.wkbLinearRing)

# Create polygon #1
poly1 = ogr.Geometry(ogr.wkbPolygon)

# Create ring #2
ring2 = ogr.Geometry(ogr.wkbLinearRing)

# Create polygon #2
poly2 = ogr.Geometry(ogr.wkbPolygon)

print multipolygon.ExportToWkt()
``````

## 创建几何体集合

``````from osgeo import ogr

# Create a geometry collection
geomcol =  ogr.Geometry(ogr.wkbGeometryCollection)

point = ogr.Geometry(ogr.wkbPoint)

line = ogr.Geometry(ogr.wkbLineString)

print geomcol.ExportToWkt()
``````

## 从WKT创建Geometry

``````from osgeo import ogr

wkt = "POINT (1120351.5712494177 741921.4223245403)"
point = ogr.CreateGeometryFromWkt(wkt)
print "%d,%d" % (point.GetX(), point.GetY())
``````

## 从GeoJSON创建Geometry

``````from osgeo import ogr

geojson = """{"type":"Point","coordinates":[108420.33,753808.59]}"""
point = ogr.CreateGeometryFromJson(geojson)
print "%d,%d" % (point.GetX(), point.GetY())
``````

## 从GML创建Geometry

``````from osgeo import ogr

gml = """<gml:Point xmlns:gml="http://www.opengis.net/gml"><gml:coordinates>108420.33,753808.59</gml:coordinates></gml:Point>"""
point = ogr.CreateGeometryFromGML(gml)
print "%d,%d" % (point.GetX(), point.GetY())
``````

## 从WKB创建Geometry

``````from osgeo import ogr
from base64 import b64decode

wkb = b64decode("AIAAAAFBMkfmVwo9cUEjylouFHrhAAAAAAAAAAA=")
point = ogr.CreateGeometryFromWkb(wkb)
print "%d,%d" % (point.GetX(), point.GetY())
``````

## 统计Geometry中的点数

``````from osgeo import ogr

wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
geom = ogr.CreateGeometryFromWkt(wkt)
print "Geometry has %i points" % (geom.GetPointCount())
``````

## 统计Geometry中的几何体数

``````from osgeo import ogr

wkt = "MULTIPOINT (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
geom = ogr.CreateGeometryFromWkt(wkt)
print "Geometry has %i geometries" % (geom.GetGeometryCount())
``````

## 遍历Geometry中的几何体

``````from osgeo import ogr

wkt = "MULTIPOINT (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
geom = ogr.CreateGeometryFromWkt(wkt)
for i in range(0, geom.GetGeometryCount()):
g = geom.GetGeometryRef(i)
print "%i). %s" %(i, g.ExportToWkt())
``````

## 遍历Geometry中的点

``````from osgeo import ogr

wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
geom = ogr.CreateGeometryFromWkt(wkt)
for i in range(0, geom.GetPointCount()):
# GetPoint returns a tuple not a Geometry
pt = geom.GetPoint(i)
print "%i). POINT (%d %d)" %(i, pt[0], pt[1])
``````

## 求Geometry的缓冲区

``````from osgeo import ogr

wkt = "POINT (1198054.34 648493.09)"
pt = ogr.CreateGeometryFromWkt(wkt)
bufferDistance = 500
poly = pt.Buffer(bufferDistance)
print "%s buffered by %d is %s" % (pt.ExportToWkt(), bufferDistance, poly.ExportToWkt())
``````

## 计算Geometry包围盒

``````from osgeo import ogr

wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
geom = ogr.CreateGeometryFromWkt(wkt)
# Get Envelope returns a tuple (minX, maxX, minY, maxY)
env = geom.GetEnvelope()
print "minX: %d, minY: %d, maxX: %d, maxY: %d" %(env[0],env[2],env[1],env[3])
``````

## 计算Geometry的面积

``````from osgeo import ogr

wkt = "POLYGON ((1162440.5712740074 672081.4332727483, 1162440.5712740074 647105.5431482664, 1195279.2416228633 647105.5431482664, 1195279.2416228633 672081.4332727483, 1162440.5712740074 672081.4332727483))"
poly = ogr.CreateGeometryFromWkt(wkt)
print "Area = %d" % poly.GetArea()
``````

## 计算Geometry的长度

``````from osgeo import ogr

wkt = "LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)"
geom = ogr.CreateGeometryFromWkt(wkt)
print "Length = %d" % geom.Length()
``````

## 获取Geometry的类型字符串

``````from osgeo import ogr

wkts = [
"POINT (1198054.34 648493.09)",
"LINESTRING (1181866.263593049 615654.4222507705, 1205917.1207499576 623979.7189589312, 1227192.8790041457 643405.4112779726, 1224880.2965852122 665143.6860159477)",
"POLYGON ((1162440.5712740074 672081.4332727483, 1162440.5712740074 647105.5431482664, 1195279.2416228633 647105.5431482664, 1195279.2416228633 672081.4332727483, 1162440.5712740074 672081.4332727483))"
]

for wkt in wkts:
geom = ogr.CreateGeometryFromWkt(wkt)
print geom.GetGeometryName()
``````

## 计算两个Geometry的交集

``````from osgeo import ogr

wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)

intersection = poly1.Intersection(poly2)

print intersection.ExportToWkt()
``````

## 计算两个Geometry的并集

``````from osgeo import ogr

wkt1 = "POLYGON ((1208064.271243039 624154.6783778917, 1208064.271243039 601260.9785661874, 1231345.9998651114 601260.9785661874, 1231345.9998651114 624154.6783778917, 1208064.271243039 624154.6783778917))"
wkt2 = "POLYGON ((1199915.6662253144 633079.3410163528, 1199915.6662253144 614453.958118695, 1219317.1067437078 614453.958118695, 1219317.1067437078 633079.3410163528, 1199915.6662253144 633079.3410163528)))"

poly1 = ogr.CreateGeometryFromWkt(wkt1)
poly2 = ogr.CreateGeometryFromWkt(wkt2)

union = poly1.Union(poly2)

print poly1
print poly2
print union.ExportToWkt()
``````

## Geometry保存为GeoJSON

``````from osgeo import ogr

# Create test polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
poly = ogr.Geometry(ogr.wkbPolygon)

# Create the output Driver
outDriver = ogr.GetDriverByName('GeoJSON')

# Create the output GeoJSON
outDataSource = outDriver.CreateDataSource('test.geojson')
outLayer = outDataSource.CreateLayer('test.geojson', geom_type=ogr.wkbPolygon )

# Get the output Layer's Feature Definition
featureDefn = outLayer.GetLayerDefn()

# create a new feature
outFeature = ogr.Feature(featureDefn)

# Set new geometry
outFeature.SetGeometry(poly)

# Add new feature to output Layer
outLayer.CreateFeature(outFeature)

# destroy the feature
outFeature.Destroy

# Close DataSources
outDataSource.Destroy()

from osgeo import ogr

# Create test polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
poly = ogr.Geometry(ogr.wkbPolygon)

geojson = poly.ExportToJson()
print geojson
``````

## Geometry保存为WKT

``````from osgeo import ogr

# Create test polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
geom_poly = ogr.Geometry(ogr.wkbPolygon)

# Export geometry to WKT
wkt = geom_poly.ExportToWkt()
print wkt
``````

## Geometry保存为KML

``````from osgeo import ogr

# Create test polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
geom_poly = ogr.Geometry(ogr.wkbPolygon)

kml = geom_poly.ExportToKML()
print kml
``````

## Geometry保存为WKB

``````from osgeo import ogr

# Create test polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
geom_poly = ogr.Geometry(ogr.wkbPolygon)

# Export geometry to WKT
wkb = geom_poly.ExportToWkb()
print wkb
``````

## polygon转成multipolygon

``````from osgeo import ogr

# Given a test polygon
poly_wkt= "POLYGON ((1179091.164690328761935 712782.883845978067257,1161053.021822647424415 667456.268434881232679,1214704.933941904921085 641092.828859039116651,1228580.428455505985767 682719.312399842427112,1218405.065812198445201 721108.180554138729349,1179091.164690328761935 712782.883845978067257))"
geom_poly = ogr.CreateGeometryFromWkt(poly_wkt)

# Force a polygon to multipolygon
if geom_poly.GetGeometryType() == ogr.wkbPolygon:
geom_poly = ogr.ForceToMultiPolygon(geom_poly)
# if are iterating over features just to update the geometry
# to multipolygon you can update the geometry using "feature.SetGeometryDirectly(geom_poly)"

# Then export geometry to WKT
wkt = geom_poly.ExportToWkt()
print wkt
``````

## 四分多边形并且创建中心点

``````import ogr

# Given a test polygon
poly_Wkt= "POLYGON((-107.42631019589980212 40.11971708125970082,-107.42455436683293613 40.12061219666851741,-107.42020981542387403 40.12004414402532859,-107.41789122063043749 40.12149008687303819,-107.41419947746419439 40.11811617239460048,-107.41915181585792993 40.11761695654455906,-107.41998470913324581 40.11894245264452508,-107.42203317637793702 40.1184088144647788,-107.42430674991324224 40.1174448122981957,-107.42430674991324224 40.1174448122981957,-107.42631019589980212 40.11971708125970082))"
geom_poly = ogr.CreateGeometryFromWkt(poly_Wkt)
``````
``````# Create 4 square polygons
geom_poly_envelope = geom_poly.GetEnvelope()
minX = geom_poly_envelope[0]
minY = geom_poly_envelope[2]
maxX = geom_poly_envelope[1]
maxY = geom_poly_envelope[3]

'''
coord0----coord1----coord2
|           |           |
coord3----coord4----coord5
|           |           |
coord6----coord7----coord8
'''
coord0 = minX, maxY
coord1 = minX+(maxX-minX)/2, maxY
coord2 = maxX, maxY
coord3 = minX, minY+(maxY-minY)/2
coord4 = minX+(maxX-minX)/2, minY+(maxY-minY)/2
coord5 = maxX, minY+(maxY-minY)/2
coord6 = minX, minY
coord7 = minX+(maxX-minX)/2, minY
coord8 = maxX, minY

ringTopLeft = ogr.Geometry(ogr.wkbLinearRing)
polyTopLeft = ogr.Geometry(ogr.wkbPolygon)

ringTopRight = ogr.Geometry(ogr.wkbLinearRing)
polyTopRight = ogr.Geometry(ogr.wkbPolygon)

ringBottomLeft = ogr.Geometry(ogr.wkbLinearRing)
polyBottomLeft = ogr.Geometry(ogr.wkbPolygon)

ringBottomRight = ogr.Geometry(ogr.wkbLinearRing)
polyBottomRight = ogr.Geometry(ogr.wkbPolygon)
``````
``````# Intersect 4 squares polygons with test polygon
quaterPolyTopLeft = polyTopLeft.Intersection(geom_poly)
quaterPolyTopRight =  polyTopRight.Intersection(geom_poly)
quaterPolyBottomLeft =  polyBottomLeft.Intersection(geom_poly)
quaterPolyBottomRight =  polyBottomRight.Intersection(geom_poly)
``````
``````# Create centroids of each intersected polygon
centroidTopLeft = quaterPolyTopLeft.Centroid()
centroidTopRight =  quaterPolyTopRight.Centroid()
centroidBottomLeft =  quaterPolyBottomLeft.Centroid()
centroidBottomRight =  quaterPolyBottomRight.Centroid()
``````