昨天要處理一個shp文件，讀取里面的信息，做個計算然后寫到后面新建的field里面。先寫個外面網(wǎng)上都能找到的新建和讀取吧。

1.讀取shp文件

#-*- coding: cp936 -*-try: from osgeo import gdal from osgeo import ogrexceptImportError: import gdal import ogr defReadVectorFile(): # 為了支持中文路徑，請?zhí)砑酉旅孢@句代碼 gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8','NO') # 為了使屬性表字段支持中文，請?zhí)砑酉旅孢@句 gdal.SetConfigOption('SHAPE_ENCODING','') strVectorFile ='E:DatumGDALCsTestDebugbeijing.shp' # 注冊所有的驅(qū)動 ogr.RegisterAll() #打開數(shù)據(jù) ds = ogr.Open(strVectorFile, 0) if ds == None: print('打開文件【%s】失敗！', strVectorFile) return print('打開文件【%s】成功！', strVectorFile) # 獲取該數(shù)據(jù)源中的圖層個數(shù)，一般shp數(shù)據(jù)圖層只有一個，如果是mdb、dxf等圖層就會有多個 iLayerCount = ds.GetLayerCount() # 獲取第一個圖層 oLayer = ds.GetLayerByIndex(0) if oLayer == None: print('獲取第%d個圖層失敗！n', 0) return # 對圖層進(jìn)行初始化，如果對圖層進(jìn)行了過濾操作，執(zhí)行這句后，之前的過濾全部清空 oLayer.ResetReading() # 通過屬性表的SQL語句對圖層中的要素進(jìn)行篩選，這部分詳細(xì)參考SQL查詢章節(jié)內(nèi)容 oLayer.SetAttributeFilter(''NAME99'LIKE '北京市市轄區(qū)'') # 通過指定的幾何對象對圖層中的要素進(jìn)行篩選 #oLayer.SetSpatialFilter() # 通過指定的四至范圍對圖層中的要素進(jìn)行篩選 #oLayer.SetSpatialFilterRect() # 獲取圖層中的屬性表表頭并輸出 print('屬性表結(jié)構(gòu)信息：') 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())) # 輸出圖層中的要素個數(shù) print('要素個數(shù) = %d', oLayer.GetFeatureCount(0)) oFeature = oLayer.GetNextFeature() # 下面開始遍歷圖層中的要素 while oFeature is not None: print('當(dāng)前處理第%d個: n屬性值：', oFeature.GetFID()) # 獲取要素中的屬性表內(nèi)容 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('數(shù)據(jù)集關(guān)閉！')

2.新建shp文件

#-*- coding: cp936 -*-try: from osgeo import gdal from osgeo import ogrexceptImportError: import gdal import ogr defWriteVectorFile(): # 為了支持中文路徑，請?zhí)砑酉旅孢@句代碼 gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8','NO') # 為了使屬性表字段支持中文，請?zhí)砑酉旅孢@句 gdal.SetConfigOption('SHAPE_ENCODING','') strVectorFile ='E:TestPolygon.shp' # 注冊所有的驅(qū)動 ogr.RegisterAll() # 創(chuàng)建數(shù)據(jù)，這里以創(chuàng)建ESRI的shp文件為例 strDriverName = 'ESRIShapefile' oDriver =ogr.GetDriverByName(strDriverName) if oDriver == None: print('%s 驅(qū)動不可用！n', strDriverName) return # 創(chuàng)建數(shù)據(jù)源 oDS =oDriver.CreateDataSource(strVectorFile) if oDS == None: print('創(chuàng)建文件【%s】失敗！', strVectorFile) return # 創(chuàng)建圖層，創(chuàng)建一個多邊形圖層，這里沒有指定空間參考，如果需要的話，需要在這里進(jìn)行指定 papszLCO = [] oLayer =oDS.CreateLayer('TestPolygon', None, ogr.wkbPolygon, papszLCO) if oLayer == None: print('圖層創(chuàng)建失敗！n') return # 下面創(chuàng)建屬性表 # 先創(chuàng)建一個叫FieldID的整型屬性 oFieldID =ogr.FieldDefn('FieldID', ogr.OFTInteger) oLayer.CreateField(oFieldID, 1) # 再創(chuàng)建一個叫FeatureName的字符型屬性，字符長度為50 oFieldName =ogr.FieldDefn('FieldName', ogr.OFTString) oFieldName.SetWidth(100) oLayer.CreateField(oFieldName, 1) oDefn = oLayer.GetLayerDefn() # 創(chuàng)建三角形要素 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) # 創(chuàng)建矩形要素 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) # 創(chuàng)建五角形要素 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('數(shù)據(jù)集創(chuàng)建完成！n')

3.更新

其實更新無非就是獲取到field然后設(shè)置新值就可以了

其實用SetField()方法就行

import os,sysfrom osgeo import gdalfrom osgeo import ogrfrom osgeo import osrimport numpyimport transformer# 為了支持中文路徑，請?zhí)砑酉旅孢@句代碼 pathname = sys.argv[1]choose = sys.argv[2]gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8', 'NO')# 為了使屬性表字段支持中文，請?zhí)砑酉旅孢@句gdal.SetConfigOption('SHAPE_ENCODING', '')# 注冊所有的驅(qū)動ogr.RegisterAll()# 數(shù)據(jù)格式的驅(qū)動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();# 輸出圖層中的要素個數(shù)print ’要素個數(shù)=%d’%(layer0.GetFeatureCount(0))print ’屬性表結(jié)構(gòu)信息’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: # 獲取要素中的屬性表內(nèi)容 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，已經(jīng)設(shè)置了，所以不需要set。

網(wǎng)上的教程都是新建和讀取，都沒有提到這個，結(jié)果自己蠢到試了好久都沒有發(fā)現(xiàn)問題在哪，以為是什么數(shù)據(jù)類型與設(shè)置字段屬性不匹配，一頭霧水哈哈哈。

補充知識：python使用GDAL生成shp文件

GDAL是一個開源的地理工具包，其支持基本所有的地理操作，其有python、java、c等語言包，是地理信息C端開發(fā)不可越過的工具，鑒于python語言的簡單性，這里使用python中GDAL包來進(jìn)行shp文件的生成，這里本質(zhì)是利用ogc地理標(biāo)準(zhǔn)的坐標(biāo)字符串來生成shp。

第一步：安裝GDAL環(huán)境，建議下載后，本地安裝，注意與python版本號要對應(yīng)，可參考網(wǎng)上教程。

第二部：代碼分析

引入GDAL工具包

import osgeo.ogr as ogrimport osgeo.osr as osr

注冊驅(qū)動，這里是ESRI Shapefile類型，并設(shè)置shp文件名稱

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

注入投影信息，這里使用的4326，表示經(jīng)緯度坐標(biāo)，根據(jù)情況可以自行更改

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

這里定義的是，生成的要素類型，包括點、線、面

#ogr.wkbPoint 點#ogr.wkbLineString 線#ogr.wkbMultiPolygon 面

這里的圖層名稱要與上面注冊驅(qū)動的shp名稱一致

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

這里設(shè)置要素的屬性字段，其中設(shè)置了兩個字段，分別是Name、data，其中ogr.OFTString表示字符串類型，其長度都是14字節(jié)，可自行設(shè)置寬度

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數(shù)據(jù)

其中各要素格式如下：

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)

需要注意的是，這里應(yīng)該與上面定義的生成要素的類型保持一致，最后是清空緩存，這里多說一句，字符串語法與postgis等開源gis一致，都遵循ogc國際標(biāo)準(zhǔn)

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

結(jié)果如下：

用arcgis打開

可以使用該方法，下載在線shp數(shù)據(jù)，只需要知道所需要素的geojson格式數(shù)據(jù)中坐標(biāo)串即可。或者圖像識別中獲取的矢量邊界賦予經(jīng)緯度。

以上這篇python使用gdal對shp讀取,新建和更新的實例就是小編分享給大家的全部內(nèi)容了，希望能給大家一個參考，也希望大家多多支持好吧啦網(wǎng)。