GDAL讀寫向量檔案——Python

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

圖1 Python建立向量結果

#-*- coding: cp936 -*-try:         from osgeo import gdal         from osgeo import ogrexceptImportError:         import gdal         import ogr

1．讀取向量

#-*- coding: cp936 -*-try:         from osgeo import gdal         from osgeo import ogrexceptImportError:         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 ogrexceptImportError:         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)