昨天要處理一個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","")
 
   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',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 ogr
import 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,20 25),(-5 -8,-10 -8,-15 -4))
MULTIPOLYGON(((1 1,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,20 25)'
line = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometry(line)
layer.CreateFeature(feature)

feature = None
data_source = None

結果如下：

用arcgis開啟

可以使用該方法，下載線上shp資料，只需要知道所需要素的geojson格式資料中座標串即可。或者影象識別中獲取的向量邊界賦予經緯度。

以上這篇python使用gdal對shp讀取,新建和更新的例項就是小編分享給大家的全部內容了，希望能給大家一個參考，也希望大家多多支援我們。