GDAL_OGR简单的使用教程

检测安装

from osgeo import gdal

查看版本

gdal.VersionInfo('VERSION_NUM')

# '2040100'

开启python异常

默认情况下，发生错误时，GDAL/OGR Python绑定不会引发异常。相反，它们返回错误值（例如None），并将错误消息写入sys.stdout。你可以通过调用UseExceptions()函数来开启异常：

from osgeo import gdal

# 开启异常
gdal.UseExceptions()

# 打开不存在的数据集
ds = gdal.Open('test.tif')

# 开启异常前
ERROR 4: test.tif: No such file or directory
# 开启异常后
RuntimeError                              Traceback (most recent call last)
<ipython-input-5-6ef000fdc647> in <module>
----> 1 gdal.Open("test.tif")

c:\program files\python37\lib\site-packages\osgeo\gdal.py in Open(*args)
   3114 def Open(*args):
   3115     """Open(char const * utf8_path, GDALAccess eAccess) -> Dataset"""
-> 3116     return _gdal.Open(*args)
   3117
   3118 def OpenEx(*args, **kwargs):

RuntimeError: test.tif: No such file or directory

你可以在运行的任何时候使用以下命令禁用GDAL/OGR异常：

gdal.DontUseExceptions()

安装GDAL OGR错误处理

安装GDAL错误处理程序功能，以捕获GDAL错误、类和消息。仅适用于GDAL 1.10以上的版本。

from osgeo import ogr, osr, gdal

# GDAL错误处理方法
def gdal_error_handler(err_class, err_num, err_msg):
    errtype = {
            gdal.CE_None:'None',
            gdal.CE_Debug:'Debug',
            gdal.CE_Warning:'Warning',
            gdal.CE_Failure:'Failure',
            gdal.CE_Fatal:'Fatal'
    }
    err_msg = err_msg.replace('\n',' ')
    err_class = errtype.get(err_class, 'None')
    print('Error Number: %s' % (err_num))
    print('Error Type: %s' % (err_class))
    print('Error Message: %s' % (err_msg))

if __name__=='__main__':

    # 安装错误处理
    gdal.PushErrorHandler(gdal_error_handler)

    # 抛出一个假的错误
    gdal.Error(1, 2, 'test error')

    # 卸载错误处理
    gdal.PopErrorHandler()

栅格图层

关闭栅格数据

from osgeo import gdal

# 打开数据
ds = gdal.Open('test.tif')

# 关闭数据（非必须）
ds = None

获取元数据

from osgeo import gdal
gtif = gdal.Open( "merge.tif" )
print (gtif.GetMetadata())

波段融合

from osgeo import gdal
# 打开数据
ds1=gdal.Open("LM51130321998341HAJ00_B1.TIF")
ds2=gdal.Open("LM51130321998341HAJ00_B2.TIF")
ds3=gdal.Open("LM51130321998341HAJ00_B3.TIF")
ds4=gdal.Open("LM51130321998341HAJ00_B4.TIF")
# 获取波段
b1=ds1.GetRasterBand(1)
b2=ds1.GetRasterBand(1)
b3=ds1.GetRasterBand(1)
b4=ds1.GetRasterBand(1)
# 创建数据
driver=gdal.GetDriverByName("GTiff")
out_ds=driver.Create('merge.tif',ds1.RasterXSize,ds1.RasterYSize,4,gdal.GDT_Byte)
# 设置坐标和投影
out_ds.SetGeoTransform(ds1.GetGeoTransform())
out_ds.SetProjection(ds1.GetProjection())

ob1=out_ds.GetRasterBand(1)
ob2=out_ds.GetRasterBand(2)
ob3=out_ds.GetRasterBand(3)
ob4=out_ds.GetRasterBand(4)
# 写入数据
ob1.WriteArray(b1.ReadAsArray())
ob2.WriteArray(b2.ReadAsArray())
ob3.WriteArray(b3.ReadAsArray())
ob4.WriteArray(b4.ReadAsArray())

out_ds=None

获取波段

from osgeo import gdal
import sys
# 允许python异常
gdal.UseExceptions()

try:
    src_ds = gdal.Open( "merge.tif" )
except RuntimeError as e:
    print(e)
    sys.exit(1)

try:
    srcband = src_ds.GetRasterBand(1)
except RuntimeError as e:
    # 试一下 GetRasterBand(10)
    print(e)
    sys.exit(1)

遍历波段

from osgeo import gdal
import sys

src_ds = gdal.Open( "merge.tif" )

print ("波段个数: ", src_ds.RasterCount)
for band in range( src_ds.RasterCount ):
    band += 1
    print ("获取波段: ", band)
    srcband = src_ds.GetRasterBand(band)
    if srcband is None:
        continue

    stats = srcband.GetStatistics( True, True )
    if stats is None:
        continue

    print ("[ STATS ] =  Minimum=%.3f, Maximum=%.3f, Mean=%.3f, StdDev=%.3f" % ( \
                stats[0], stats[1], stats[2], stats[3] ))

获取波段信息

from osgeo import gdal
import sys
gdal.UseExceptions()


def main( band_num, input_file ):
    src_ds = gdal.Open( input_file )

    try:
        srcband = src_ds.GetRasterBand(band_num)
    except RuntimeError as e:
        print( e)
        sys.exit(1)


    print ("[ NO DATA VALUE ] = ", srcband.GetNoDataValue())
    print ("[ MIN ] = ", srcband.GetMinimum())
    print ("[ MAX ] = ", srcband.GetMaximum())
    print ("[ SCALE ] = ", srcband.GetScale())
    print( "[ UNIT TYPE ] = ", srcband.GetUnitType())
    ctable = srcband.GetColorTable()

    if ctable is None:
        print ('No ColorTable found')
        sys.exit(1)

    print ("[ COLOR TABLE COUNT ] = ", ctable.GetCount())
    for i in range( 0, ctable.GetCount() ):
        entry = ctable.GetColorEntry( i )
        if not entry:
            continue
        print ("[ COLOR ENTRY RGB ] = ", ctable.GetColorEntryAsRGB( i, entry ))

main( 1,"merge.tif" )

多边形化栅格波段

from osgeo import gdal, ogr
import sys

gdal.UseExceptions()

src_ds = gdal.Open( "test.tif" )

srcband = src_ds.GetRasterBand(3)

dst_layername = "POLYGONIZED_STUFF"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

栅格化矢量数据

from osgeo import gdal, ogr

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

vector_fn = 'test.shp'

raster_fn = 'test.tif'

# 打开数据并读取范围
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
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('GTiff').Create(raster_fn, x_res, y_res, 1, 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=[0])

用shapefile裁剪栅格

from osgeo import gdal, gdalnumeric, ogr, osr
from PIL import Image, ImageDraw
import os, sys
gdal.UseExceptions()


def imageToArray(i):
    """
    Image转换为数组.
    """
    a=gdalnumeric.fromstring(i.tostring(),'b')
    a.shape=i.im.size[1], i.im.size[0]
    return a

def arrayToImage(a):
    """
   数组转换为Image.
    """
    i=Image.fromstring('L',(a.shape[1],a.shape[0]),
            (a.astype('b')).tostring())
    return i

def world2Pixel(geoMatrix, x, y):
  """
  用地理仿射变换计算像素坐标
  """
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  rtnX = geoMatrix[2]
  rtnY = geoMatrix[4]
  pixel = int((x - ulX) / xDist)
  line = int((ulY - y) / xDist)
  return (pixel, line)

#
#  EDIT: this is basically an overloaded
#  version of the gdal_array.OpenArray passing in xoff, yoff explicitly
#  so we can pass these params off to CopyDatasetInfo
#
def OpenArray( array, prototype_ds = None, xoff=0, yoff=0 ):
    ds = gdal.Open( gdalnumeric.GetArrayFilename(array) )

    if ds is not None and prototype_ds is not None:
        if type(prototype_ds).__name__ == 'str':
            prototype_ds = gdal.Open( prototype_ds )
        if prototype_ds is not None:
            gdalnumeric.CopyDatasetInfo( prototype_ds, ds, xoff=xoff, yoff=yoff )
    return ds

def histogram(a, bins=range(0,256)):
  """
  Histogram function for multi-dimensional array.
  a = array
  bins = range of numbers to match
  """
  fa = a.flat
  n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)
  n = gdalnumeric.concatenate([n, [len(fa)]])
  hist = n[1:]-n[:-1]
  return hist

def stretch(a):
  """
  Performs a histogram stretch on a gdalnumeric array image.
  """
  hist = histogram(a)
  im = arrayToImage(a)
  lut = []
  for b in range(0, len(hist), 256):
    # step size
    step = reduce(operator.add, hist[b:b+256]) / 255
    # create equalization lookup table
    n = 0
    for i in range(256):
      lut.append(n / step)
      n = n + hist[i+b]
  im = im.point(lut)
  return imageToArray(im)

def main( shapefile_path, raster_path ):
    # 读取数据到数组
    srcArray = gdalnumeric.LoadFile(raster_path)

    # 获取地理仿射变换
    srcImage = gdal.Open(raster_path)
    geoTrans = srcImage.GetGeoTransform()

    # 获取矢量图层
    shapef = ogr.Open(shapefile_path)
    lyr = shapef.GetLayer( os.path.split( os.path.splitext( shapefile_path )[0] )[1] )
    poly = lyr.GetNextFeature()

    # 获取范围
    minX, maxX, minY, maxY = lyr.GetExtent()
    ulX, ulY = world2Pixel(geoTrans, minX, maxY)
    lrX, lrY = world2Pixel(geoTrans, maxX, minY)

    # 计算像素大小
    pxWidth = int(lrX - ulX)
    pxHeight = int(lrY - ulY)

    clip = srcArray[:, ulY:lrY, ulX:lrX]

    #
    # 像素偏移
    #
    xoffset =  ulX
    yoffset =  ulY
    print ("Xoffset, Yoffset = ( %f, %f )" % ( xoffset, yoffset ))

    # 创建新的仿射变换
    geoTrans = list(geoTrans)
    geoTrans[0] = minX
    geoTrans[3] = maxY

    # Map points to pixels for drawing the
    # boundary on a blank 8-bit,
    # black and white, mask image.
    points = []
    pixels = []
    geom = poly.GetGeometryRef()
    pts = geom.GetGeometryRef(0)
    for p in range(pts.GetPointCount()):
      points.append((pts.GetX(p), pts.GetY(p)))
    for p in points:
      pixels.append(world2Pixel(geoTrans, p[0], p[1]))
    rasterPoly = Image.new("L", (pxWidth, pxHeight), 1)
    rasterize = ImageDraw.Draw(rasterPoly)
    rasterize.polygon(pixels, 0)
    mask = imageToArray(rasterPoly)

    # Clip the image using the mask
    clip = gdalnumeric.choose(mask, \
        (clip, 0)).astype(gdalnumeric.uint8)

    # This image has 3 bands so we stretch each one to make them
    # visually brighter
    for i in range(3):
      clip[i,:,:] = stretch(clip[i,:,:])

    # Save new tiff
    #
    #  EDIT: instead of SaveArray, let's break all the
    #  SaveArray steps out more explicity so
    #  we can overwrite the offset of the destination
    #  raster
    #
    ### the old way using SaveArray
    #
    # gdalnumeric.SaveArray(clip, "OUTPUT.tif", format="GTiff", prototype=raster_path)
    #
    ###
    #
    gtiffDriver = gdal.GetDriverByName( 'GTiff' )
    if gtiffDriver is None:
        raise ValueError("Can't find GeoTiff Driver")
    gtiffDriver.CreateCopy( "OUTPUT.tif",
        OpenArray( clip, prototype_ds=raster_path, xoff=xoffset, yoff=yoffset )
    )

    # Save as an 8-bit jpeg for an easy, quick preview
    clip = clip.astype(gdalnumeric.uint8)
    gdalnumeric.SaveArray(clip, "OUTPUT.jpg", format="JPEG")

    gdal.ErrorReset()


if __name__ == '__main__':

    #
    # example run : $ python clip.py /<full-path>/<shapefile-name>.shp /<full-path>/<raster-name>.tif
    #
    if len( sys.argv ) < 2:
        print "[ ERROR ] you must two args. 1) the full shapefile path and 2) the full raster path"
        sys.exit( 1 )

    main( sys.argv[1], sys.argv[2] )

区域统计

import gdal, ogr, osr, numpy
import sys


def zonal_stats(feat, input_zone_polygon, input_value_raster):

    # Open data
    raster = gdal.Open(input_value_raster)
    shp = ogr.Open(input_zone_polygon)
    lyr = shp.GetLayer()

    # Get raster georeference info
    transform = raster.GetGeoTransform()
    xOrigin = transform[0]
    yOrigin = transform[3]
    pixelWidth = transform[1]
    pixelHeight = transform[5]

    # Reproject vector geometry to same projection as raster
    sourceSR = lyr.GetSpatialRef()
    targetSR = osr.SpatialReference()
    targetSR.ImportFromWkt(raster.GetProjectionRef())
    coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
    feat = lyr.GetNextFeature()
    geom = feat.GetGeometryRef()
    geom.Transform(coordTrans)

    # Get extent of feat
    geom = feat.GetGeometryRef()
    if (geom.GetGeometryName() == 'MULTIPOLYGON'):
        count = 0
        pointsX = []; pointsY = []
        for polygon in geom:
            geomInner = geom.GetGeometryRef(count)
            ring = geomInner.GetGeometryRef(0)
            numpoints = ring.GetPointCount()
            for p in range(numpoints):
                    lon, lat, z = ring.GetPoint(p)
                    pointsX.append(lon)
                    pointsY.append(lat)
            count += 1
    elif (geom.GetGeometryName() == 'POLYGON'):
        ring = geom.GetGeometryRef(0)
        numpoints = ring.GetPointCount()
        pointsX = []; pointsY = []
        for p in range(numpoints):
                lon, lat, z = ring.GetPoint(p)
                pointsX.append(lon)
                pointsY.append(lat)

    else:
        sys.exit("ERROR: Geometry needs to be either Polygon or Multipolygon")

    xmin = min(pointsX)
    xmax = max(pointsX)
    ymin = min(pointsY)
    ymax = max(pointsY)

    # Specify offset and rows and columns to read
    xoff = int((xmin - xOrigin)/pixelWidth)
    yoff = int((yOrigin - ymax)/pixelWidth)
    xcount = int((xmax - xmin)/pixelWidth)+1
    ycount = int((ymax - ymin)/pixelWidth)+1

    # Create memory target raster
    target_ds = gdal.GetDriverByName('MEM').Create('', xcount, ycount, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((
        xmin, pixelWidth, 0,
        ymax, 0, pixelHeight,
    ))

    # Create for target raster the same projection as for the value raster
    raster_srs = osr.SpatialReference()
    raster_srs.ImportFromWkt(raster.GetProjectionRef())
    target_ds.SetProjection(raster_srs.ExportToWkt())

    # Rasterize zone polygon to raster
    gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])

    # Read raster as arrays
    banddataraster = raster.GetRasterBand(1)
    dataraster = banddataraster.ReadAsArray(xoff, yoff, xcount, ycount).astype(numpy.float)

    bandmask = target_ds.GetRasterBand(1)
    datamask = bandmask.ReadAsArray(0, 0, xcount, ycount).astype(numpy.float)

    # Mask zone of raster
    zoneraster = numpy.ma.masked_array(dataraster,  numpy.logical_not(datamask))

    # Calculate statistics of zonal raster
    return numpy.average(zoneraster),numpy.mean(zoneraster),numpy.median(zoneraster),numpy.std(zoneraster),numpy.var(zoneraster)


def loop_zonal_stats(input_zone_polygon, input_value_raster):

    shp = ogr.Open(input_zone_polygon)
    lyr = shp.GetLayer()
    featList = range(lyr.GetFeatureCount())
    statDict = {}

    for FID in featList:
        feat = lyr.GetFeature(FID)
        meanValue = zonal_stats(feat, input_zone_polygon, input_value_raster)
        statDict[FID] = meanValue
    return statDict

def main(input_zone_polygon, input_value_raster):
    return loop_zonal_stats(input_zone_polygon, input_value_raster)


if __name__ == "__main__":

    #
    # Returns for each feature a dictionary item (FID) with the statistical values in the following order: Average, Mean, Medain, Standard Deviation, Variance
    #
    # example run : $ python grid.py <full-path><output-shapefile-name>.shp xmin xmax ymin ymax gridHeight gridWidth
    #

    if len( sys.argv ) != 3:
        print "[ ERROR ] you must supply two arguments: input-zone-shapefile-name.shp input-value-raster-name.tif "
        sys.exit( 1 )
    print 'Returns for each feature a dictionary item (FID) with the statistical values in the following order: Average, Mean, Medain, Standard Deviation, Variance'
    print main( sys.argv[1], sys.argv[2] )

从数组创建栅格

import gdal, ogr, os, osr
import numpy as np


def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):

    cols = array.shape[1]
    rows = array.shape[0]
    originX = rasterOrigin[0]
    originY = rasterOrigin[1]

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
    # 反转数组
    reversed_arr = array[::-1]
    # 数组转栅格
    array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) 


if __name__ == "__main__":
    rasterOrigin = (-123.25745,45.43013)
    pixelWidth = 10
    pixelHeight = 10
    newRasterfn = 'test.tif'
    array = np.array([[ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 1, 1],
                      [ 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
                      [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])


    main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)

替换无效值

import gdal, ogr, osr, os
import numpy as np

def raster2array(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.ReadAsArray()

def getNoDataValue(rasterfn):
    raster = gdal.Open(rasterfn)
    band = raster.GetRasterBand(1)
    return band.GetNoDataValue()

def array2raster(rasterfn,newRasterfn,array):
    raster = gdal.Open(rasterfn)
    geotransform = raster.GetGeoTransform()
    originX = geotransform[0]
    originY = geotransform[3]
    pixelWidth = geotransform[1]
    pixelHeight = geotransform[5]
    cols = raster.RasterXSize
    rows = raster.RasterYSize

    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    outRasterSRS.ImportFromWkt(raster.GetProjectionRef())
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


rasterfn = 'test.tif'
newValue = 0
newRasterfn = 'testNew.tif'

# 栅格转数组
rasterArray = raster2array(rasterfn)

# 获取无效值
noDataValue = getNoDataValue(rasterfn)

# 更新无效值
rasterArray[rasterArray == noDataValue] = newValue

# 数组转栅格
array2raster(rasterfn,newRasterfn,rasterArray)

投影

创建投影

from osgeo import osr
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326) #WGS84

重投影

from osgeo import ogr
from osgeo import osr

source = osr.SpatialReference()
source.ImportFromEPSG(2927)

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

transform = osr.CoordinateTransformation(source, target)

point = ogr.CreateGeometryFromWkt("POINT (1120351.57 741921.42)")
point.Transform(transform)

print (point.ExportToWkt())

获取投影

from osgeo import ogr, osr
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open(r'./data/test.shp')

# 图层
layer = dataset.GetLayer()
spatialRef = layer.GetSpatialRef()
# 几何
feature = layer.GetNextFeature()
geom = feature.GetGeometryRef()
spatialRef = geom.GetSpatialReference()

重投影图层

from osgeo import ogr, osr
import os

driver = ogr.GetDriverByName('ESRI Shapefile')

# 输入空间参考
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(2927)

# 输出空间参考
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(4326)

# 创建坐标转换
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# 输入图层
inDataSet = driver.Open(r'./data/test.shp')
inLayer = inDataSet.GetLayer()

# 创建输出图层
outputShapefile = r'./data/test_4326.shp'
if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)
outDataSet = driver.CreateDataSource(outputShapefile)
outLayer = outDataSet.CreateLayer(
    "test_4326", 
    geom_type=ogr.wkbMultiPolygon,
    srs=outSpatialRef # 输出prj文件
)

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

outLayerDefn = outLayer.GetLayerDefn()

# 遍历要素
inFeature = inLayer.GetNextFeature()
while inFeature:
    # 获取几何
    geom = inFeature.GetGeometryRef()
    # 重投影几何
    geom.Transform(coordTrans)
    # 创建要素
    outFeature = ogr.Feature(outLayerDefn)
    # 设置几何、属性
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # 添加要素
    outLayer.CreateFeature(outFeature)
    # 取消引用
    outFeature = None
    inFeature = inLayer.GetNextFeature()

# 保存并关闭
inDataSet = None
outDataSet = None

输出投影

from osgeo import ogr, osr
driver = ogr.GetDriverByName('ESRI Shapefile')
dataset = driver.Open(r'./data/test_4326.shp')
layer = dataset.GetLayer()
spatialRef = layer.GetSpatialRef()

spatialRef.ExportToWkt()
spatialRef.ExportToPrettyWkt()
spatialRef.ExportToPCI()
spatialRef.ExportToUSGS()
spatialRef.ExportToXML()

创建ESRI投影文件

from osgeo import ogr, osr

spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(26912)

spatialRef.MorphToESRI()
file = open('yourshpfile.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()

创建质心

centroidTopLeft = quaterPolyTopLeft.Centroid()
centroidTopRight = quaterPolyTopRight.Centroid()
centroidBottomLeft = quaterPolyBottomLeft.Centroid()
centroidBottomRight = quaterPolyBottomRight.Centroid()

矢量图层

删除文件

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')