矢量图层

删除文件

from osgeo import ogr
import os

DriverName = "ESRI Shapefile"      # e.g.: GeoJSON, ESRI Shapefile
FileName = 'test.shp'
driver = ogr.GetDriverByName(DriverName)
if os.path.exists(FileName):
     driver.DeleteDataSource(FileName)

获取OGR驱动列表

import ogr
cnt = ogr.GetDriverCount()
formatsList = []  # Empty List

for i in range(cnt):
    driver = ogr.GetDriver(i)
    driverName = driver.GetName()
    if not driverName in formatsList:
        formatsList.append(driverName)

formatsList.sort() # 排序

print(formatsList)

驱动是否可用

from osgeo import ogr

## Shapefile 是否可用?
driverName = "ESRI Shapefile"
drv = ogr.GetDriverByName( driverName )
if drv is None:
    print("%s 驱动不可用.\n" % driverName)
else:
    print  ("%s 驱动可用.\n" % driverName)

## PostgreSQL 是否可用?
driverName = "PostgreSQL"
drv = ogr.GetDriverByName( driverName )
if drv is None:
    print ("%s 驱动不可用.\n" % driverName)
else:
    print  ("%s 驱动可用.\n" % driverName)

## File GeoDatabase 是否可用?
driverName = "FileGDB"
drv = ogr.GetDriverByName( driverName )
if drv is None:
    print ("%s 驱动不可用.\n" % driverName)
else:
    print  ("%s 驱动可用.\n" % driverName)

## SDE 是否可用?
driverName = "SDE"
drv = ogr.GetDriverByName( driverName )
if drv is None:
    print ("%s 驱动不可用.\n" % driverName)
else:
    print  ("%s 驱动可用.\n" % driverName)

获取shapefile要素个数

import os
from osgeo import ogr

daShapefile = r"/test.shp"

driver = ogr.GetDriverByName('ESRI Shapefile')

dataSource = driver.Open(daShapefile, 0) # 0 只读. 1 读写.

# 检查数据源是否有效.
if dataSource is None:
    print ('不能打开 %s' % (daShapefile))
else:
    print ('打开 %s' % (daShapefile))
    layer = dataSource.GetLayer()
    featureCount = layer.GetFeatureCount()
    print ("%s 要素个数: %d" % (os.path.basename(daShapefile),featureCount))

获取PostGIS图层

from osgeo import ogr

databaseServer = "localhost"
databaseName = "test2020"
databaseUser = "postgres"
databasePW = "123456"


connString = "PG: host=%s dbname=%s user=%s password=%s" %(databaseServer,databaseName,databaseUser,databasePW)

conn = ogr.Open(connString)

layerList = []
for i in conn:
    daLayer = i.GetName()
    if not daLayer in layerList:
        layerList.append(daLayer)

layerList.sort()

for j in layerList:
    print(j)

# 关闭连接
conn = None

获取PostGIS图层的要素

from osgeo import ogr
import sys

databaseServer = "localhost"
databaseName = "test2020"
databaseUser = "postgres"
databasePW = "123456"
connString = "PG: host=%s dbname=%s user=%s password=%s" % (databaseServer,databaseName,databaseUser,databasePW)

def GetPGLayer( lyr_name ):
    conn = ogr.Open(connString)

    lyr = conn.GetLayer( lyr_name )
    if lyr is None:
        sys.exit( 1 )

    featureCount = lyr.GetFeatureCount()
    print ("%s 要素个数 : %d" % ( lyr_name , featureCount ))

    # 关闭连接
    conn = None


if __name__ == '__main__':
    GetPGLayer( "test" )

获取ESRI GDB图层

import sys
from osgeo import ogr

ogr.UseExceptions()

driver = ogr.GetDriverByName("OpenFileGDB")

# opening the FileGDB
try:
    gdb = driver.Open("/sparse.gdb", 0)
except Exception as e:
    print(e)
    sys.exit()

featsClassList = []

# 获取图层
for featsClass_idx in range(gdb.GetLayerCount()):
    featsClass = gdb.GetLayerByIndex(featsClass_idx)
    featsClassList.append(featsClass.GetName())

featsClassList.sort()

for featsClass in featsClassList:
    print (featsClass)

加载数据到内存

from osgeo import ogr

# 打开输入数据源
indriver=ogr.GetDriverByName('SQLite')
srcdb = indriver.Open('/poly_spatialite.sqlite',0)

# 创建内存输出数据源
outdriver=ogr.GetDriverByName('MEMORY')
source=outdriver.CreateDataSource('memData')

# 打开内存数据
tmp=outdriver.Open('memData',1)

# 复制图层到内存
poly_mem=source.CopyLayer(srcdb.GetLayer('poly'),'poly',['OVERWRITE=YES'])

# 新的图层可直接被访问，poly_mem 或者 source.GetLayer('poly'):
layer=source.GetLayer('poly')
for feature in layer:
    feature.SetField('area',1)

遍历要素

from osgeo import ogr
import os

shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()

for feature in layer:
    print( feature.GetField("area"))
layer.ResetReading()

遍历要素几何

from osgeo import ogr
import os

shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()

for feature in layer:
    geom = feature.GetGeometryRef()
    # 获取质心
    print (geom.Centroid().ExportToWkt())

过滤属性

from osgeo import ogr
import os

shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()

layer.SetAttributeFilter("area = 5268.813")

for feature in layer:
    print (feature.GetField("area"))

空间过滤

from osgeo import ogr
import os

shapefile = "/test.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shapefile, 0)
layer = dataSource.GetLayer()

wkt = "POLYGON ((479386 4764749,481098 4764226,480772 4763114,478681 4763159,479386 4764749))"
layer.SetSpatialFilter(ogr.CreateGeometryFromWkt(wkt))

for feature in layer:
    print (feature.GetField("area"))

获取要素字段

from osgeo import ogr

daShapefile = r"/test.shp"

dataSource = ogr.Open(daShapefile)
daLayer = dataSource.GetLayer(0)
layerDefinition = daLayer.GetLayerDefn()


for i in range(layerDefinition.GetFieldCount()):
    fieldName =  layerDefinition.GetFieldDefn(i).GetName()
    fieldTypeCode = layerDefinition.GetFieldDefn(i).GetType()
    fieldType = layerDefinition.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
    fieldWidth = layerDefinition.GetFieldDefn(i).GetWidth()
    GetPrecision = layerDefinition.GetFieldDefn(i).GetPrecision()

    print (fieldName + " - " + fieldType+ " " + str(fieldWidth) + " " + str(GetPrecision))

获取PostGIS图层字段

from osgeo import ogr
import sys

databaseServer = "localhost"
databaseName = "test2020"
databaseUser = "postgres"
databasePW = "123456"
connString = "PG: host=%s dbname=%s user=%s password=%s" %(databaseServer,databaseName,databaseUser,databasePW)


def GetPGLayerFields( lyr_name ):
    conn = ogr.Open(connString)

    lyr = conn.GetLayer( lyr_name )
    if lyr is None:
        sys.exit( 1 )

    lyrDefn = lyr.GetLayerDefn()


    for i in range( lyrDefn.GetFieldCount() ):
        fieldName =  lyrDefn.GetFieldDefn(i).GetName()
        fieldTypeCode = lyrDefn.GetFieldDefn(i).GetType()
        fieldType = lyrDefn.GetFieldDefn(i).GetFieldTypeName(fieldTypeCode)
        fieldWidth = lyrDefn.GetFieldDefn(i).GetWidth()
        GetPrecision = lyrDefn.GetFieldDefn(i).GetPrecision()

        print (fieldName + " - " + fieldType+ " " + str(fieldWidth) + " " + str(GetPrecision))


if __name__ == '__main__':
    GetPGLayerFields( "test" )

获取图层能力

from osgeo import ogr

ds = ogr.Open("/test.shp",0)
layer = ds.GetLayer()
capabilities = [
    ogr.OLCRandomRead,
    ogr.OLCSequentialWrite,
    ogr.OLCRandomWrite,
    ogr.OLCFastSpatialFilter,
    ogr.OLCFastFeatureCount,
    ogr.OLCFastGetExtent,
    ogr.OLCCreateField,
    ogr.OLCDeleteField,
    ogr.OLCReorderFields,
    ogr.OLCAlterFieldDefn,
    ogr.OLCTransactions,
    ogr.OLCDeleteFeature,
    ogr.OLCFastSetNextByIndex,
    ogr.OLCStringsAsUTF8,
    ogr.OLCIgnoreFields
]

print("图层能力:")
for cap in capabilities:
    print("  %s = %s" % (cap, layer.TestCapability(cap)))

WFS图层和遍历要素

import sys

from osgeo import ogr, osr, gdal

# 获取WFS驱动
wfs_drv = ogr.GetDriverByName('WFS')

# 加快查询多图层WFS服务
gdal.SetConfigOption('OGR_WFS_LOAD_MULTIPLE_LAYER_DEFN', 'NO')

# 设置分页的配置。适用于WFS 2.0服务以及WFS 1.0和1.1以及其他一些服务。
gdal.SetConfigOption('OGR_WFS_PAGING_ALLOWED', 'YES')
gdal.SetConfigOption('OGR_WFS_PAGE_SIZE', '10000')

url = 'http://sampleserver6.arcgisonline.com/arcgis/services/SampleWorldCities/MapServer/WFSServer'
wfs_ds = wfs_drv.Open('WFS:' + url)
if not wfs_ds:
    sys.exit('错误: 不能打开 WFS 数据源')
else:
    pass

# 遍历图层
for i in range(wfs_ds.GetLayerCount()):
    layer = wfs_ds.GetLayerByIndex(i)
    srs = layer.GetSpatialRef()
    print ('Layer: %s, Features: %s, SR: %s...' % (layer.GetName(), layer.GetFeatureCount(), srs.ExportToWkt()[0:50]))

    # 遍历要素
    feat = layer.GetNextFeature()
    while feat is not None:
        feat = layer.GetNextFeature()
        # do something more..
    feat = None

# 获取指定图层
layer = wfs_ds.GetLayerByName("esri:World")
if not layer:
    sys.exit('错误：不能找到图层：esri:World')
else:
    pass

设置HTTP代理

import sys

from osgeo import ogr, osr, gdal

server = 'proxy.example.com'
port = 3128

# 设置代理
gdal.SetConfigOption('GDAL_HTTP_PROXY', server + ':' + port)

# 没有用户名或密码的NTLM设置代理身份验证选项，因此单点登录有效
gdal.SetConfigOption('GDAL_PROXY_AUTH', 'NTLM')
gdal.SetConfigOption('GDAL_HTTP_PROXYUSERPWD', ' : ')

ds = ogr.Open('http://featureserver/cities/.geojson')
if not ds:
    sys.exit('ERROR: can not open GeoJSON datasource')
lyr = ds.GetLayer('OGRGeoJSON')
for feat in lyr:
    geom = feat.GetGeometryRef()
    print( geom.ExportToWkt())

读取CSV经纬度作为OGRVRTLayer

GDAL/OGR具有虚拟格式规范，该规范允许你从诸如CSV之类的平面表派生图层——它的功能远不止于此，因此请继续阅读。在下面的示例中，我们正在读取带有X、Y列和值的CSV。该CSV文件由XML文件包装，该XML文件将其描述为OGR层。以下是所有必要的部分和一个脚本，该脚本读取XML文件并打印出点的几何形状。

CSV文件：

ID,X,Y
1,-127.234343,47.234325
2,-127.003243,46.234343
3,-127.345646,45.234324
4,-126.234324,44.324234

XML文件

<OGRVRTDataSource>
    <OGRVRTLayer name="example">
        <SrcDataSource>example.csv</SrcDataSource>
        <SrcLayer>example</SrcLayer>
        <GeometryType>wkbPoint</GeometryType>
            <LayerSRS>WGS84</LayerSRS>
        <GeometryField encoding="PointFromColumns" x="X" y="Y"/>
    </OGRVRTLayer>
</OGRVRTDataSource>

from osgeo import ogr
ogr.UseExceptions()

inDataSource = ogr.Open("example_wrapper.vrt")
lyr = inDataSource.GetLayer('example')
for feat in lyr:
    geom = feat.GetGeometryRef()
    print (geom.ExportToWkt())

计算范围

from osgeo import ogr
import os

# Get a Layer's Extent
inShapefile = "states.shp"
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(inShapefile, 0)
inLayer = inDataSource.GetLayer()
extent = inLayer.GetExtent()

# 创建多边形
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(extent[0],extent[2])
ring.AddPoint(extent[1], extent[2])
ring.AddPoint(extent[1], extent[3])
ring.AddPoint(extent[0], extent[3])
ring.AddPoint(extent[0],extent[2])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

# 保存到新的shp文件
outShapefile = "new.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")

# 如果存在，先删除
if os.path.exists(outShapefile):
    outDriver.DeleteDataSource(outShapefile)

# 创建数据源
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("new", geom_type=ogr.wkbPolygon)

# 添加ID字段
idField = ogr.FieldDefn("id", ogr.OFTInteger)
outLayer.CreateField(idField)

# 创建要素
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(poly)
feature.SetField("id", 1)
outLayer.CreateFeature(feature)
feature = None

# 保存并关闭
inDataSource = None
outDataSource = None

计算凸包

from osgeo import ogr
import os

# 获得图层
inShapefile = "test.shp"
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(inShapefile, 0)
inLayer = inDataSource.GetLayer()

# 几何集合
geomcol = ogr.Geometry(ogr.wkbGeometryCollection)
for feature in inLayer:
    geomcol.AddGeometry(feature.GetGeometryRef())

# 计算凸包
convexhull = geomcol.ConvexHull()

# 保存
outShapefile = "test_convexhull.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")

# 如果存在，先删除
if os.path.exists(outShapefile):
    outDriver.DeleteDataSource(outShapefile)

# 输出
outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("test_convexhull", geom_type=ogr.wkbPolygon)

# 添加ID字段
idField = ogr.FieldDefn("id", ogr.OFTInteger)
outLayer.CreateField(idField)

# 创建要素
featureDefn = outLayer.GetLayerDefn()
feature = ogr.Feature(featureDefn)
feature.SetGeometry(convexhull)
feature.SetField("id", 1)
outLayer.CreateFeature(feature)
feature = None

# 保存并关闭
inDataSource = None
outDataSource = None

计算质心

from osgeo import ogr
import os

ogr.UseExceptions()

# 输入图层
inShapefile = "test.shp"
inDriver = ogr.GetDriverByName("ESRI Shapefile")
inDataSource = inDriver.Open(inShapefile, 0)
inLayer = inDataSource.GetLayer()

# 输出图层
outShapefile = "test_centroids.shp"
outDriver = ogr.GetDriverByName("ESRI Shapefile")

# 如果存在，先删除
if os.path.exists(outShapefile):
    outDriver.DeleteDataSource(outShapefile)

outDataSource = outDriver.CreateDataSource(outShapefile)
outLayer = outDataSource.CreateLayer("test_centroids", geom_type=ogr.wkbPoint)

# 添加字段
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)

# 获得要素定义
outLayerDefn = outLayer.GetLayerDefn()

# 添加要素
for i in range(0, inLayer.GetFeatureCount()):
    # 输入要素
    inFeature = inLayer.GetFeature(i)
    # 输出要素
    outFeature = ogr.Feature(outLayerDefn)
    # 设置字段值
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # 设置质心几何
    geom = inFeature.GetGeometryRef()
    inFeature = None
    centroid = geom.Centroid()
    outFeature.SetGeometry(centroid)
    # 添加新要素
    outLayer.CreateFeature(outFeature)
    outFeature = None

# 保存并关闭
inDataSource = None
outDataSource = None

创建新的shp数据并添加数据

import osgeo.ogr as ogr
import osgeo.osr as osr
import csv,os

# 使用字典读取数据
reader = csv.DictReader(open("volcano_data.txt","r"),
    delimiter='\t',
    quoting=csv.QUOTE_NONE)

# 驱动
outShapefile="volcanoes.shp"
driver = ogr.GetDriverByName("ESRI Shapefile")

if os.path.exists(outShapefile):
    driver.DeleteDataSource(outShapefile)

# 创建数据源
data_source = driver.CreateDataSource(outShapefile)

# 创建空间参考 WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

# 创建图层
layer = data_source.CreateLayer("volcanoes", srs, ogr.wkbPoint)

# 添加字段
field_name = ogr.FieldDefn("Name", ogr.OFTString)
field_name.SetWidth(24)
layer.CreateField(field_name)
field_region = ogr.FieldDefn("Region", ogr.OFTString)
field_region.SetWidth(24)
layer.CreateField(field_region)
layer.CreateField(ogr.FieldDefn("Latitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Longitude", ogr.OFTReal))
layer.CreateField(ogr.FieldDefn("Elevation", ogr.OFTInteger))

# 处理文本
for row in reader:
  # 创建要素
  feature = ogr.Feature(layer.GetLayerDefn())
  # 设置属性字段
  feature.SetField("Name", row['Name'])
  feature.SetField("Region", row['Region'])
  feature.SetField("Latitude", row['Latitude'])
  feature.SetField("Longitude", row['Longitude'])
  feature.SetField("Elevation", row['Elevation'])

  # 创建WKT
  wkt = "POINT(%f %f)" %  (float(row['Longitude']) , float(row['Latitude']))

  # 创建点
  point = ogr.CreateGeometryFromWkt(wkt)

  # 设置几何
  feature.SetGeometry(point)
  # 添加要素
  layer.CreateFeature(feature)
  # 删除引用
  feature = None

# 保存关闭
data_source = None

从WKT创建PostGIS表

import ogr, osr

database = 'test2020'
usr = 'postgres'
pw = '123456'
table = 'testtest'

wkt = "POINT (1120351.5712494177 741921.4223245403)"
point = ogr.CreateGeometryFromWkt(wkt)

connectionString = "PG:dbname='%s' user='%s' password='%s'" % (database,usr,pw)
ogrds = ogr.Open(connectionString)

srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)

layer = ogrds.CreateLayer(table, srs, ogr.wkbPoint, ['OVERWRITE=YES'] )

layerDefn = layer.GetLayerDefn()

feature = ogr.Feature(layerDefn)
feature.SetGeometry(point)

layer.StartTransaction()
layer.CreateFeature(feature)
feature = None
layer.CommitTransaction()

过滤和选择

ogr2ogr -f "ESRI Shapefile" junkmob.shp -select area -where "area = 5268.813" test.shp

# 该命令读取parcel_address.shp并生成junkmob.shp，area=5268.813输出area列

from osgeo import ogr
import os, sys

def main( field_name_target ):
    # 输入图层
    inShapefile = "test.shp"
    inDriver = ogr.GetDriverByName("ESRI Shapefile")
    inDataSource = inDriver.Open(inShapefile, 0)
    inLayer = inDataSource.GetLayer()
    inLayer.SetAttributeFilter("area = 5268.813")

    # 创建输出图层
    outShapefile = os.path.join( os.path.split( inShapefile )[0], "junkmob.shp" )
    outDriver = ogr.GetDriverByName("ESRI Shapefile")

    # 存在，先删除
    if os.path.exists(outShapefile):
        outDriver.DeleteDataSource(outShapefile)

    # 创建输出shp
    outDataSource = outDriver.CreateDataSource(outShapefile)
    out_lyr_name = os.path.splitext( os.path.split( outShapefile )[1] )[0]
    outLayer = outDataSource.CreateLayer( out_lyr_name, geom_type=ogr.wkbMultiPolygon )

    # 添加字段
    inLayerDefn = inLayer.GetLayerDefn()
    for i in range(0, inLayerDefn.GetFieldCount()):
        fieldDefn = inLayerDefn.GetFieldDefn(i)
        fieldName = fieldDefn.GetName()
        if fieldName not in field_name_target:
            continue
        outLayer.CreateField(fieldDefn)

    # 要素定义
    outLayerDefn = outLayer.GetLayerDefn()

    # 添加要素
    for inFeature in inLayer:
        # 创建要素
        outFeature = ogr.Feature(outLayerDefn)

        # 添加字段
        for i in range(0, outLayerDefn.GetFieldCount()):
            fieldDefn = outLayerDefn.GetFieldDefn(i)
            fieldName = fieldDefn.GetName()
            if fieldName not in field_name_target:
                continue

            outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(),
                inFeature.GetField(i))

        # 设置几何
        geom = inFeature.GetGeometryRef()
        outFeature.SetGeometry(geom.Clone())
        # 创建要素
        outLayer.CreateFeature(outFeature)
        outFeature = None

    # 保存关闭
    inDataSource = None
    outDataSource = None


main( ["AREA","EAS_ID"])

合并图层

import os, ogr, osr

outputMergefn = 'merge.shp'
directory = "/Users/UserName/Downloads/"
fileStartsWith = 'test'
fileEndsWith = '.shp'
driverName = 'ESRI Shapefile'
geometryType = ogr.wkbPolygon

out_driver = ogr.GetDriverByName( driverName )
if os.path.exists(outputMergefn):
    out_driver.DeleteDataSource(outputMergefn)
out_ds = out_driver.CreateDataSource(outputMergefn)
out_layer = out_ds.CreateLayer(outputMergefn, geom_type=geometryType)

fileList = os.listdir(directory)

for file in fileList:
    if file.startswith(fileStartsWith) and file.endswith(fileEndsWith):
        print file
        ds = ogr.Open(directory+file)
        lyr = ds.GetLayer()
        for feat in lyr:
            out_feat = ogr.Feature(out_layer.GetLayerDefn())
            out_feat.SetGeometry(feat.GetGeometryRef().Clone())
            out_layer.CreateFeature(out_feat)
            out_feat = None
            out_layer.SyncToDisk()

获取OSM街道名称

TODO：测试

import ogr

ds = ogr.Open('test.osm')
layer = ds.GetLayer()

nameList = []
for feature in layer:
    if feature.GetField("highway") != None:
        name = feature.GetField("name")
        if name != None and name not in nameList:
            nameList.append(name)

print (nameList)

创建鱼网

import os, sys
import ogr
from math import ceil


def main(outputGridfn,xmin,xmax,ymin,ymax,gridHeight,gridWidth):
    xmin = float(xmin)
    xmax = float(xmax)
    ymin = float(ymin)
    ymax = float(ymax)
    gridWidth = float(gridWidth)
    gridHeight = float(gridHeight)

    # get rows
    rows = ceil((ymax-ymin)/gridHeight)
    # get columns
    cols = ceil((xmax-xmin)/gridWidth)

    # start grid cell envelope
    ringXleftOrigin = xmin
    ringXrightOrigin = xmin + gridWidth
    ringYtopOrigin = ymax
    ringYbottomOrigin = ymax-gridHeight

    # create output file
    outDriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputGridfn):
        os.remove(outputGridfn)
    outDataSource = outDriver.CreateDataSource(outputGridfn)
    outLayer = outDataSource.CreateLayer(outputGridfn,geom_type=ogr.wkbPolygon )
    featureDefn = outLayer.GetLayerDefn()

    # create grid cells
    countcols = 0
    while countcols < cols:
        countcols += 1

        # reset envelope for rows
        ringYtop = ringYtopOrigin
        ringYbottom =ringYbottomOrigin
        countrows = 0

        while countrows < rows:
            countrows += 1
            ring = ogr.Geometry(ogr.wkbLinearRing)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYtop)
            ring.AddPoint(ringXrightOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYbottom)
            ring.AddPoint(ringXleftOrigin, ringYtop)
            poly = ogr.Geometry(ogr.wkbPolygon)
            poly.AddGeometry(ring)

            # add new geom to layer
            outFeature = ogr.Feature(featureDefn)
            outFeature.SetGeometry(poly)
            outLayer.CreateFeature(outFeature)
            outFeature = None

            # new envelope for next poly
            ringYtop = ringYtop - gridHeight
            ringYbottom = ringYbottom - gridHeight

        # new envelope for next poly
        ringXleftOrigin = ringXleftOrigin + gridWidth
        ringXrightOrigin = ringXrightOrigin + gridWidth

    # Save and close DataSources
    outDataSource = None


if __name__ == "__main__":

    #
    # example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth
    #

    if len( sys.argv ) != 8:
        print "[ ERROR ] you must supply seven arguments: output-shapefile-name.shp xmin xmax ymin ymax gridHeight gridWidth"
        sys.exit( 1 )

    main( sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7] )

面转线

import ogr, os

def poly2line(input_poly,output_line):

    source_ds = ogr.Open(input_poly)
    source_layer = source_ds.GetLayer()

    # polygon2geometryCollection
    geomcol =  ogr.Geometry(ogr.wkbGeometryCollection)
    for feat in source_layer:
        geom = feat.GetGeometryRef()
        ring = geom.GetGeometryRef(0)
        geomcol.AddGeometry(ring)

    # geometryCollection2shp
    shpDriver = ogr.GetDriverByName("ESRI Shapefile")
    if os.path.exists(output_line):
            shpDriver.DeleteDataSource(output_line)
    outDataSource = shpDriver.CreateDataSource(output_line)
    outLayer = outDataSource.CreateLayer(output_line, geom_type=ogr.wkbMultiLineString)
    featureDefn = outLayer.GetLayerDefn()
    outFeature = ogr.Feature(featureDefn)
    outFeature.SetGeometry(geomcol)
    outLayer.CreateFeature(outFeature)
    outFeature = None

def main(input_poly,output_line):
    poly2line(input_poly,output_line)

if __name__ == "__main__":
    input_poly = 'test_polygon.shp'
    output_line = 'test_line.shp'

    main(input_poly,output_line)

创建缓冲区

import ogr, os

def createBuffer(inputfn, outputBufferfn, bufferDist):
    inputds = ogr.Open(inputfn)
    inputlyr = inputds.GetLayer()

    shpdriver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists(outputBufferfn):
        shpdriver.DeleteDataSource(outputBufferfn)
    outputBufferds = shpdriver.CreateDataSource(outputBufferfn)
    bufferlyr = outputBufferds.CreateLayer(outputBufferfn, geom_type=ogr.wkbPolygon)
    featureDefn = bufferlyr.GetLayerDefn()

    for feature in inputlyr:
        ingeom = feature.GetGeometryRef()
        geomBuffer = ingeom.Buffer(bufferDist)

        outFeature = ogr.Feature(featureDefn)
        outFeature.SetGeometry(geomBuffer)
        bufferlyr.CreateFeature(outFeature)
        outFeature = None

def main(inputfn, outputBufferfn, bufferDist):
    createBuffer(inputfn, outputBufferfn, bufferDist)


if __name__ == "__main__":
    inputfn = 'test.shp'
    outputBufferfn = 'testBuffer.shp'
    bufferDist = 10.0

    main(inputfn, outputBufferfn, bufferDist)

栅格化矢量图层

import ogr, gdal

vector_fn = 'test.shp'

# 定义像素大小和无效值
pixel_size = 25
NoData_value = 255

# 打开数据源，读取数据范围
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
source_srs = source_layer.GetSpatialRef()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# 创建目标数据源
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('MEM').Create('', x_res, y_res, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# 栅格化
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])

# 读取为数组
array = band.ReadAsArray()
print (array)

面转点

TODO：测试

import ogr, gdal
import numpy as np
import os

polygon_fn = 'test.shp'

# Define pixel_size which equals distance betweens points
pixel_size = 10

# Open the data source and read in the extent
source_ds = ogr.Open(polygon_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()

# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create('temp.tif', x_res, y_res, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(255)

# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[1])

# Read as array
array = band.ReadAsArray()

raster = gdal.Open('temp.tif')
geotransform = raster.GetGeoTransform()

# Convert array to point coordinates
count = 0
roadList = np.where(array == 1)
multipoint = ogr.Geometry(ogr.wkbMultiPoint)
for indexY in roadList[0]:
    indexX = roadList[1][count]
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    Xcoord = originX+pixelWidth*indexX
    Ycoord = originY+pixelHeight*indexY
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(Xcoord, Ycoord)
    multipoint.AddGeometry(point)
    count += 1

# Write point coordinates to Shapefile
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists('points.shp'):
    shpDriver.DeleteDataSource('points.shp')
outDataSource = shpDriver.CreateDataSource('points.shp')
outLayer = outDataSource.CreateLayer('points.shp', geom_type=ogr.wkbMultiPoint)
featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(multipoint)
outLayer.CreateFeature(outFeature)
outFeature = None

# Remove temporary files
os.remove('temp.tif')