python使用gdal对shp读取,新建和更新的实例

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

1.读取shp文件

#-*- coding: cp936 -*-try: from osgeo import gdal from osgeo import ogrexceptImportError: import gdal import ogr defReadVectorFile(): # 为了支持中文路径，请添加下面这句代码 gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8','NO') # 为了使属性表字段支持中文，请添加下面这句 gdal.SetConfigOption('SHAPE_ENCODING','') strVectorFile ='E:DatumGDALCsTestDebugbeijing.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 ogrexceptImportError: 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,sysfrom osgeo import gdalfrom osgeo import ogrfrom osgeo import osrimport numpyimport 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)#新建fieldfieldDefn = 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。

网上的教程都是新建和读取，都没有提到这个，结果自己蠢到试了好久都没有发现问题在哪，以为是什么数据类型与设置字段属性不匹配，一头雾水哈哈哈。

补充知识：python使用GDAL生成shp文件

GDAL是一个开源的地理工具包，其支持基本所有的地理操作，其有python、java、c等语言包，是地理信息C端开发不可越过的工具，鉴于python语言的简单性，这里使用python中GDAL包来进行shp文件的生成，这里本质是利用ogc地理标准的坐标字符串来生成shp。

第一步：安装GDAL环境，建议下载后，本地安装，注意与python版本号要对应，可参考网上教程。

第二部：代码分析

引入GDAL工具包

import osgeo.ogr as ogrimport osgeo.osr as osr

注册驱动，这里是ESRI Shapefile类型，并设置shp文件名称

driver = ogr.GetDriverByName('ESRI Shapefile')data_source = driver.CreateDataSource('ceshi.shp')

注入投影信息，这里使用的4326，表示经纬度坐标，根据情况可以自行更改

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

这里定义的是，生成的要素类型，包括点、线、面

#ogr.wkbPoint 点#ogr.wkbLineString 线#ogr.wkbMultiPolygon 面

这里的图层名称要与上面注册驱动的shp名称一致

layer = data_source.CreateLayer('ceshi', srs, ogr.wkbLineString)

这里设置要素的属性字段，其中设置了两个字段，分别是Name、data，其中ogr.OFTString表示字符串类型，其长度都是14字节，可自行设置宽度

field_name = ogr.FieldDefn('Name', ogr.OFTString)field_name.SetWidth(14)layer.CreateField(field_name)

field_name = ogr.FieldDefn('data', ogr.OFTString)field_name.SetWidth(14)layer.CreateField(field_name)

在生成的字段名中插入要素值，即属性表中每行的值

feature = ogr.Feature(layer.GetLayerDefn())feature.SetField('Name', 'ceshi')feature.SetField('data', '1.2')

核心部分，生成line数据

其中各要素格式如下：

POINT(6 10)LINESTRING(3 4,10 50,20 25)POLYGON((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2))MULTIPOINT(3.5 5.6, 4.8 10.5)MULTILINESTRING((3 4,10 50,20 25),(-5 -8,-10 -8,-15 -4))MULTIPOLYGON(((1 1,5 1,5 5,1 5,1 1),(2 2,2 3,3 3,3 2,2 2)),((6 3,9 2,9 4,6 3)))GEOMETRYCOLLECTION(POINT(4 6),LINESTRING(4 6,7 10))POINT ZM (1 1 5 60)POINT M (1 1 80)

需要注意的是，这里应该与上面定义的生成要素的类型保持一致，最后是清空缓存，这里多说一句，字符串语法与postgis等开源gis一致，都遵循ogc国际标准

wkt = ’LINESTRING(3 4,10 50,20 25)’line = ogr.CreateGeometryFromWkt(wkt)feature.SetGeometry(line)layer.CreateFeature(feature)feature = Nonedata_source = None

结果如下：

用arcgis打开

可以使用该方法，下载在线shp数据，只需要知道所需要素的geojson格式数据中坐标串即可。或者图像识别中获取的矢量边界赋予经纬度。

以上这篇python使用gdal对shp读取,新建和更新的实例就是小编分享给大家的全部内容了，希望能给大家一个参考，也希望大家多多支持好吧啦网。