python使用gdal对shp读取，新建和更新

python使用gdal对shp读取，新建和更新

昨天要处理一个shp文件，读取里面的信息，做个计算然后写到后面新建的field里面。先写个外面网上都能找到的新建和读取吧。

1.读取shp文件

#-*- coding: cp936 -*-
try:from osgeo import gdalfrom osgeo import ogr
exceptImportError:import gdalimport ogrdefReadVectorFile():# 为了支持中文路径，请添加下面这句代码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)returnprint("打开文件【%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()# 为了演示，只输出一个要素信息breakprint("数据集关闭！")

2.新建shp文件

#-*- coding: cp936 -*-
try:from osgeo import gdalfrom osgeo import ogr
exceptImportError:import gdalimport ogrdefWriteVectorFile():# 为了支持中文路径，请添加下面这句代码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的字符型属性，字符长度为50oFieldName =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'%it(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。网上的教程都是新建和读取，都没有提到这个，结果自己蠢到试了好久都没有发现问题在哪，以为是什么数据类型与设置字段属性不匹配，一头雾水哈哈哈。

最后求关注，求点赞，欢迎大家关注我的公众号

记录所学所用，包括但不限于遥感、地信、气象、生态环境，机器学习知识，相关文献阅读，编程代码实现。偶尔荒腔走板的聊聊其他。欢迎不同领域的朋友们加入进来，多多交流。