gdal python | 对shp读取，新建和更新

1.读取shp文件

#-*- coding: cp936 -*-
try:
         from osgeo import gdal
         from osgeo import ogr
exceptImportError:
         import gdal
         import ogr
 
defReadVectorFile():
         # 为了支持中文路径，请添加下面这句代码
         gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
         # 为了使属性表字段支持中文，请添加下面这句
         gdal.SetConfigOption("SHAPE_ENCODING","")
 
         strVectorFile ="E:\\Datum\\GDALCsTest\\Debug\\beijing.shp"
 
         # 注册所有的驱动
         ogr.RegisterAll()
 
         #打开数据
         ds = ogr.Open(strVectorFile, 0)
         if ds == None:
                   print("打开文件【%s】失败！", strVectorFile)
                   return
 
         print("打开文件【%s】成功！", strVectorFile)
 
         # 获取该数据源中的图层个数，一般shp数据图层只有一个，如果是mdb、dxf等图层就会有多个
         iLayerCount = ds.GetLayerCount()
 
         # 获取第一个图层
         oLayer = ds.GetLayerByIndex(0)
         if oLayer == None:
                   print("获取第%d个图层失败！\n", 0)
                   return
 
         # 对图层进行初始化，如果对图层进行了过滤操作，执行这句后，之前的过滤全部清空
         oLayer.ResetReading()
 
         # 通过属性表的SQL语句对图层中的要素进行筛选，这部分详细参考SQL查询章节内容
         oLayer.SetAttributeFilter("\"NAME99\"LIKE \"北京市市辖区\"")
 
         # 通过指定的几何对象对图层中的要素进行筛选
         #oLayer.SetSpatialFilter()
 
         # 通过指定的四至范围对图层中的要素进行筛选
         #oLayer.SetSpatialFilterRect()
 
         # 获取图层中的属性表表头并输出
         print("属性表结构信息：")
         oDefn = oLayer.GetLayerDefn()
         iFieldCount = oDefn.GetFieldCount()
         for iAttr in range(iFieldCount):
                   oField =oDefn.GetFieldDefn(iAttr)
                   print( "%s: %s(%d.%d)" % ( \
                                     oField.GetNameRef(),\
                                     oField.GetFieldTypeName(oField.GetType() ), \
                                     oField.GetWidth(),\
                                     oField.GetPrecision()))
 
         # 输出图层中的要素个数
         print("要素个数 = %d", oLayer.GetFeatureCount(0))
 
         oFeature = oLayer.GetNextFeature()
         # 下面开始遍历图层中的要素
         while oFeature is not None:
                   print("当前处理第%d个: \n属性值：", oFeature.GetFID())
                   # 获取要素中的属性表内容
                   for iField inrange(iFieldCount):
                            oFieldDefn =oDefn.GetFieldDefn(iField)
                            line =  " %s (%s) = " % ( \
                                               oFieldDefn.GetNameRef(),\
                                               ogr.GetFieldTypeName(oFieldDefn.GetType()))
 
                            ifoFeature.IsFieldSet( iField ):
                                     line = line+ "%s" % (oFeature.GetFieldAsString( iField ) )
                            else:
                                     line = line+ "(null)"
 
                            print(line)
        
                   # 获取要素中的几何体
                   oGeometry =oFeature.GetGeometryRef()
 
                   # 为了演示，只输出一个要素信息
                   break
 
         print("数据集关闭！")

2.新建shp文件

#-*- coding: cp936 -*-
try:
         from osgeo import gdal
         from osgeo import ogr
exceptImportError:
         import gdal
         import ogr
 
defWriteVectorFile():
         # 为了支持中文路径，请添加下面这句代码
         gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")
         # 为了使属性表字段支持中文，请添加下面这句
         gdal.SetConfigOption("SHAPE_ENCODING","")
 
         strVectorFile ="E:\\TestPolygon.shp"
 
         # 注册所有的驱动
         ogr.RegisterAll()
 
         # 创建数据，这里以创建ESRI的shp文件为例
         strDriverName = "ESRIShapefile"
         oDriver =ogr.GetDriverByName(strDriverName)
         if oDriver == None:
                   print("%s 驱动不可用！\n", strDriverName)
                   return
        
         # 创建数据源
         oDS =oDriver.CreateDataSource(strVectorFile)
         if oDS == None:
                   print("创建文件【%s】失败！", strVectorFile)
                   return
 
         # 创建图层，创建一个多边形图层，这里没有指定空间参考，如果需要的话，需要在这里进行指定
         papszLCO = []
         oLayer =oDS.CreateLayer("TestPolygon", None, ogr.wkbPolygon, papszLCO)
         if oLayer == None:
                   print("图层创建失败！\n")
                   return
 
         # 下面创建属性表
         # 先创建一个叫FieldID的整型属性
         oFieldID =ogr.FieldDefn("FieldID", ogr.OFTInteger)
         oLayer.CreateField(oFieldID, 1)
 
         # 再创建一个叫FeatureName的字符型属性，字符长度为50
         oFieldName =ogr.FieldDefn("FieldName", ogr.OFTString)
         oFieldName.SetWidth(100)
         oLayer.CreateField(oFieldName, 1)
 
         oDefn = oLayer.GetLayerDefn()
 
         # 创建三角形要素
         oFeatureTriangle = ogr.Feature(oDefn)
         oFeatureTriangle.SetField(0, 0)
         oFeatureTriangle.SetField(1, "三角形")
         geomTriangle =ogr.CreateGeometryFromWkt("POLYGON ((0 0,20 0,10 15,0 0))")
         oFeatureTriangle.SetGeometry(geomTriangle)
         oLayer.CreateFeature(oFeatureTriangle)
 
         # 创建矩形要素
         oFeatureRectangle = ogr.Feature(oDefn)
         oFeatureRectangle.SetField(0, 1)
         oFeatureRectangle.SetField(1, "矩形")
         geomRectangle =ogr.CreateGeometryFromWkt("POLYGON ((30 0,60 0,60 30,30 30,30 0))")
         oFeatureRectangle.SetGeometry(geomRectangle)
         oLayer.CreateFeature(oFeatureRectangle)
 
         # 创建五角形要素
         oFeaturePentagon = ogr.Feature(oDefn)
         oFeaturePentagon.SetField(0, 2)
         oFeaturePentagon.SetField(1, "五角形")
         geomPentagon =ogr.CreateGeometryFromWkt("POLYGON ((70 0,85 0,90 15,80 30,65 15,700))")
         oFeaturePentagon.SetGeometry(geomPentagon)
         oLayer.CreateFeature(oFeaturePentagon)
 
         oDS.Destroy()
         print("数据集创建完成！\n")

3.更新

import os,sys
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
import numpy
import transformer
# 为了支持中文路径，请添加下面这句代码
 
 
pathname = sys.argv[1]
choose = sys.argv[2]
gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "NO")
# 为了使属性表字段支持中文，请添加下面这句
gdal.SetConfigOption("SHAPE_ENCODING", "")
# 注册所有的驱动
ogr.RegisterAll()
# 数据格式的驱动
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.Open(pathname, update=1)
if ds is None:
    print 'Could not open %s'%pathname
    sys.exit(1)
# 获取第0个图层
layer0 = ds.GetLayerByIndex(0);
# 投影
spatialRef = layer0.GetSpatialRef();
# 输出图层中的要素个数
print '要素个数=%d'%(layer0.GetFeatureCount(0))
print '属性表结构信息'
defn = layer0.GetLayerDefn()
fieldindex = defn.GetFieldIndex('x')
xfield = defn.GetFieldDefn(fieldindex)
#新建field
fieldDefn = ogr.FieldDefn('newx', xfield.GetType())
fieldDefn.SetWidth(32)
fieldDefn.SetPrecision(6)
layer0.CreateField(fieldDefn,1)
fieldDefn = ogr.FieldDefn('newy', xfield.GetType())
fieldDefn.SetWidth(32)
fieldDefn.SetPrecision(6)
layer0.CreateField(fieldDefn,1)
feature = layer0.GetNextFeature()
# 下面开始遍历图层中的要素
while feature is not None:
    # 获取要素中的属性表内容
    x = feature.GetFieldAsDouble('x')
    y = feature.GetFieldAsDouble('y')
    newx, newy = transformer.begintrans(choose, x, y)
    feature.SetField('newx', newx)
    feature.SetField('newy', newy)
    layer0.SetFeature(feature)
    feature = layer0.GetNextFeature()
feature.Destroy()
ds.Destroy()

   更新之后，必须SetFeature()，不然数据不会更新。新建的时候有createfeature，已经设置了，所以不需要set。