python中如何使用gdal對shp讀取、新建和更新

這篇文章主要介紹了python中如何使用gdal對shp讀取、新建和更新，具有一定借鑒價值，感興趣的朋友可以參考下，希望大家閱讀完這篇文章之后大有收獲，下面讓小編帶著大家一起了解一下。

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。

網上的教程都是新建和讀取，都沒有提到這個，結果自己蠢到試了好久都沒有發現問題在哪，以為是什么數據類型與設置字段屬性不匹配，一頭霧水哈哈哈。

補充知識：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,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 = None
data_source = None

結果如下：

用arcgis打開

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

感謝你能夠認真閱讀完這篇文章，希望小編分享的“python中如何使用gdal對shp讀取、新建和更新”這篇文章對大家有幫助，同時也希望大家多多支持億速云，關注億速云行業資訊頻道，更多相關知識等著你來學習!