GDAL讀寫向量檔案——Python

在Python中使用OGR時，先要匯入OGR庫，如果需要對中文的支援，還需要匯入GDAL庫，具體程式碼如下。Python建立的shp結果如圖1所示。

圖1 Python建立向量結果

1. #-*- coding: cp936 -*-

2. try:

3. from osgeo import gdal

4. from osgeo import ogr

5. exceptImportError:

6. import gdal

7. import ogr

1．讀取向量

1. #-*- coding: cp936 -*-

2. try:

3. from osgeo import gdal

4. from osgeo import ogr

5. exceptImportError:

6. import gdal

7. import ogr

8.  

9. defReadVectorFile():

10. # 為了支援中文路徑，請新增下面這句程式碼

11. gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")

12. # 為了使屬性表字段支援中文，請新增下面這句

13. gdal.SetConfigOption("SHAPE_ENCODING","")

14.  

15. strVectorFile ="E:\\Datum\\GDALCsTest\\Debug\\beijing.shp"

16.  

17. # 註冊所有的驅動

18. ogr.RegisterAll()

19.  

20. #開啟資料

21. ds = ogr.Open(strVectorFile, 0)

22. if ds == None:

23. print("開啟檔案【%s】失敗！", strVectorFile)

24. return

25.  

26. print("開啟檔案【%s】成功！", strVectorFile)

27.  

28. # 獲取該資料來源中的圖層個數，一般shp資料圖層只有一個，如果是mdb、dxf等圖層就會有多個

29. iLayerCount = ds.GetLayerCount()

30.  

31. # 獲取第一個圖層

32. oLayer = ds.GetLayerByIndex(0)

33. if oLayer == None:

34. print("獲取第%d個圖層失敗！\n", 0)

35. return

36.  

37. # 對圖層進行初始化，如果對圖層進行了過濾操作，執行這句後，之前的過濾全部清空

38. oLayer.ResetReading()

39.  

40. # 通過屬性表的SQL語句對圖層中的要素進行篩選，這部分詳細參考SQL查詢章節內容

41. oLayer.SetAttributeFilter("\"NAME99\"LIKE \"北京市市轄區\"")

42.  

43. # 通過指定的幾何物件對圖層中的要素進行篩選

44. #oLayer.SetSpatialFilter()

45.  

46. # 通過指定的四至範圍對圖層中的要素進行篩選

47. #oLayer.SetSpatialFilterRect()

48.  

49. # 獲取圖層中的屬性表表頭並輸出

50. print("屬性表結構資訊：")

51. oDefn = oLayer.GetLayerDefn()

52. iFieldCount = oDefn.GetFieldCount()

53. for iAttr in range(iFieldCount):

54. oField =oDefn.GetFieldDefn(iAttr)

55. print( "%s: %s(%d.%d)" % ( \

56. oField.GetNameRef(),\

57. oField.GetFieldTypeName(oField.GetType() ), \

58. oField.GetWidth(),\

59. oField.GetPrecision()))

60.  

61. # 輸出圖層中的要素個數

62. print("要素個數 = %d", oLayer.GetFeatureCount(0))

63.  

64. oFeature = oLayer.GetNextFeature()

65. # 下面開始遍歷圖層中的要素

66. while oFeature is not None:

67. print("當前處理第%d個: \n屬性值：", oFeature.GetFID())

68. # 獲取要素中的屬性表內容

69. for iField inrange(iFieldCount):

70. oFieldDefn =oDefn.GetFieldDefn(iField)

71. line = " %s (%s) = " % ( \

72. oFieldDefn.GetNameRef(),\

73. ogr.GetFieldTypeName(oFieldDefn.GetType()))

74.  

75. ifoFeature.IsFieldSet( iField ):

76. line = line+ "%s" % (oFeature.GetFieldAsString( iField ) )

77. else:

78. line = line+ "(null)"

79.  

80. print(line)

81.  

82. # 獲取要素中的幾何體

83. oGeometry =oFeature.GetGeometryRef()

84.  

85. # 為了演示，只輸出一個要素資訊

86. break

87.  

88. print("資料集關閉！")

執行上面的程式碼，如果不設定屬性過濾，輸出內容如圖3‑9上半部分所示，如過設定了屬性過濾，輸出內容如圖3‑9下半部分所示。（Python輸出的中文轉為編碼了）。

圖2 OGR庫使用Python讀取向量示例

2．寫入向量

在使用Python建立向量圖形的時候，使用的WKT格式的字串來進行建立。也可以使用其他的方式進行建立。程式碼如下，寫出來的向量圖形和屬性表如圖1所示。

1. #-*- coding: cp936 -*-

2. try:

3. from osgeo import gdal

4. from osgeo import ogr

5. exceptImportError:

6. import gdal

7. import ogr

8.  

9. defWriteVectorFile():

10. # 為了支援中文路徑，請新增下面這句程式碼

11. gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8","NO")

12. # 為了使屬性表字段支援中文，請新增下面這句

13. gdal.SetConfigOption("SHAPE_ENCODING","")

14.  

15. strVectorFile ="E:\\TestPolygon.shp"

16.  

17. # 註冊所有的驅動

18. ogr.RegisterAll()

19.  

20. # 建立資料，這裡以建立ESRI的shp檔案為例

21. strDriverName = "ESRIShapefile"

22. oDriver =ogr.GetDriverByName(strDriverName)

23. if oDriver == None:

24. print("%s 驅動不可用！\n", strDriverName)

25. return

26.  

27. # 建立資料來源

28. oDS =oDriver.CreateDataSource(strVectorFile)

29. if oDS == None:

30. print("建立檔案【%s】失敗！", strVectorFile)

31. return

32.  

33. # 建立圖層，建立一個多邊形圖層，這裡沒有指定空間參考，如果需要的話，需要在這裡進行指定

34. papszLCO = []

35. oLayer =oDS.CreateLayer("TestPolygon", None, ogr.wkbPolygon, papszLCO)

36. if oLayer == None:

37. print("圖層建立失敗！\n")

38. return

39.  

40. # 下面建立屬性表

41. # 先建立一個叫FieldID的整型屬性

42. oFieldID =ogr.FieldDefn("FieldID", ogr.OFTInteger)

43. oLayer.CreateField(oFieldID, 1)

44.  

45. # 再建立一個叫FeatureName的字元型屬性，字元長度為50

46. oFieldName =ogr.FieldDefn("FieldName", ogr.OFTString)

47. oFieldName.SetWidth(100)

48. oLayer.CreateField(oFieldName, 1)

49.  

50. oDefn = oLayer.GetLayerDefn()

51.  

52. # 建立三角形要素

53. oFeatureTriangle = ogr.Feature(oDefn)

54. oFeatureTriangle.SetField(0, 0)

55. oFeatureTriangle.SetField(1, "三角形")

56. geomTriangle =ogr.CreateGeometryFromWkt("POLYGON ((0 0,20 0,10 15,0 0))")

57. oFeatureTriangle.SetGeometry(geomTriangle)

58. oLayer.CreateFeature(oFeatureTriangle)

59.  

60. # 建立矩形要素

61. oFeatureRectangle = ogr.Feature(oDefn)

62. oFeatureRectangle.SetField(0, 1)

63. oFeatureRectangle.SetField(1, "矩形")

64. geomRectangle =ogr.CreateGeometryFromWkt("POLYGON ((30 0,60 0,60 30,30 30,30 0))")

65. oFeatureRectangle.SetGeometry(geomRectangle)

66. oLayer.CreateFeature(oFeatureRectangle)

67.  

68. # 建立五角形要素

69. oFeaturePentagon = ogr.Feature(oDefn)

70. oFeaturePentagon.SetField(0, 2)

71. oFeaturePentagon.SetField(1, "五角形")

72. geomPentagon =ogr.CreateGeometryFromWkt("POLYGON ((70 0,85 0,90 15,80 30,65 15,700))")

73. oFeaturePentagon.SetGeometry(geomPentagon)

74. oLayer.CreateFeature(oFeaturePentagon)

75.  

76. oDS.Destroy()

77. print("資料集建立完成！\n")

3．向量資料管理

1. defVectorDelete(strVectorFile):

2. # 註冊所有的驅動

3. ogr.RegisterAll()

4.  

5. oDriver = None

6. #開啟向量

7. oDS = ogr.Open(strVectorFile, 0)

8. if oDS == None:

9. os.remove(strVectorFile)

10. return

11.  

12. oDriver = oDS.GetDriver()

13. if oDriver == None:

14. os.remove(strVectorFile)

15. return

16.  

17. ifoDriver.DeleteDataSource(strVectorFile) == ogr.OGRERR_NONE:

18. return

19. else:

20. os.remove(strVectorFile)

21.  

22. defVectorRename(strOldFile, strNewFile):

23. # 註冊所有的驅動

24. ogr.RegisterAll()

25.  

26. oDriver = None

27. #開啟向量

28. oDS = Ogr.Open(strOldFile, 0)

29. if oDS == None :

30. os.rename(strOldFile,strNewFile)

31. return

32.  

33. oDriver = oDS.GetDriver()

34. if oDriver == None:

35. os.rename(strOldFile,strNewFile)

36. return

37.  

38. oDDS = oDriver.CopyDataSource(oDS,strNewFile, None)

39. if oDDS == None:

40. os.rename(strOldFile,strNewFile)

41.  

42. if oDriver.DeleteDataSource(strOldFile)== ogr.OGRERR_NONE:

43. return

44. else :

45. os.rename(strOldFile,strNewFile)