python使用gdal对shp读取，新建和更新
昨天要处理一个shp文件，读取里面的信息，做个计算然后写到后面新建的field里面。先写个外面网上都能找到的新建和读取吧。
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.更新
其实更新无非就是获取到field然后设置新值就可以了
其实用SetField()方法就行
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()，就是你更新好了field的feature一定要重新set一下，不然是根本起不到任何改变的。新建的时候有createfeature，已经设置了，所以不需要set。网上的教程都是新建和读取，都没有提到这个，结果自己蠢到试了好久都没有发现问题在哪，以为是什么数据类型与设置字段属性不匹配，一头雾水哈哈哈。
转载自：https://blog.csdn.net/lionhenryzxxy/article/details/86714511
