GDAL讀寫向量檔案——Python
GDAL讀寫向量檔案——Python
在Python中使用OGR時,先要匯入OGR庫,如果需要對中文的支援,還需要匯入GDAL庫,具體程式碼如下。Python建立的shp結果如圖1所示。
圖1 Python建立向量結果
-
#-*- coding: cp936 -*-
-
try:
-
from osgeo import gdal
-
from osgeo import ogr
-
exceptImportError:
-
import gdal
-
import ogr
1.讀取向量
-
#-*- 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("資料集關閉!")
執行上面的程式碼,如果不設定屬性過濾,輸出內容如圖3‑9上半部分所示,如過設定了屬性過濾,輸出內容如圖3‑9下半部分所示。(Python輸出的中文轉為編碼了)。
圖2 OGR庫使用Python讀取向量示例
2.寫入向量
在使用Python建立向量圖形的時候,使用的WKT格式的字串來進行建立。也可以使用其他的方式進行建立。程式碼如下,寫出來的向量圖形和屬性表如圖1所示。
-
#-*- 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.向量資料管理
-
defVectorDelete(strVectorFile):
-
# 註冊所有的驅動
-
ogr.RegisterAll()
-
oDriver = None
-
#開啟向量
-
oDS = ogr.Open(strVectorFile, 0)
-
if oDS == None:
-
os.remove(strVectorFile)
-
return
-
oDriver = oDS.GetDriver()
-
if oDriver == None:
-
os.remove(strVectorFile)
-
return
-
ifoDriver.DeleteDataSource(strVectorFile) == ogr.OGRERR_NONE:
-
return
-
else:
-
os.remove(strVectorFile)
-
defVectorRename(strOldFile, strNewFile):
-
# 註冊所有的驅動
-
ogr.RegisterAll()
-
oDriver = None
-
#開啟向量
-
oDS = Ogr.Open(strOldFile, 0)
-
if oDS == None :
-
os.rename(strOldFile,strNewFile)
-
return
-
oDriver = oDS.GetDriver()
-
if oDriver == None:
-
os.rename(strOldFile,strNewFile)
-
return
-
oDDS = oDriver.CopyDataSource(oDS,strNewFile, None)
-
if oDDS == None:
-
os.rename(strOldFile,strNewFile)
-
if oDriver.DeleteDataSource(strOldFile)== ogr.OGRERR_NONE:
-
return
-
else :
-
os.rename(strOldFile,strNewFile)