Python+GDAL/OGR向量資料讀寫__Python

       常見的向量資料格式有Shapefile、GeoJSON、CSV，及檔案資料庫gdb和空間資料庫PostGIS，不論是何種格式的資料或如何儲存，一旦開啟資料來源、擷取向量圖層後（詳情參考OGR操作向量資料的類結構圖），對資料的操作都一樣。下面對向量資料的讀寫進行詳細的介紹。     一、開啟不同的向量檔案

    1、定義開啟資料來源的函數，並遍曆所有的圖層，輸出他們的名字和圖層

'''    輸出資料來源中的圖層    參數：fn 資料來源的路徑         is_write 開啟資料來源的模式，0 表示唯讀模式，1 表示讀寫入模式'''def print_layers(fn, is_write):    ds = ogr.Open(fn, is_write)    if ds is None:        #raise OSError('Could not open %s', fn)        raise OSError('Could not open {}'.format(fn))    for i in range(ds.GetLayerCount()):        lyr = ds.GetLayer(i)        print('{0}:{1}'.format(i, lyr.GetName()))

    2、開啟PostGIS資料來源

    利用上面開啟資料來源並輸出圖層的方法，開啟本地的Post GIS資料庫 postgis_24_sample，輸出圖層如下：

#開啟postgis資料庫    print_layers('PG:user=postgres password=postgis dbname=postgis_24_sample', 0)

0:tiger.county1:tiger.state2:tiger.place3:tiger.cousub4:tiger.edges5:tiger.addrfeat6:tiger.faces7:tiger.zcta58:tiger.tract9:tiger.tabblock10:tiger.bg

    3、開啟檔案夾資料來源（shapefiles和CSV)

    開啟shapefiles檔案夾，輸出的結果如下所示：

#開啟shapefile檔案夾    shp_fn = r'E:\HZZ_GXDP_Spatial_Data'    print_layers(shp_fn, 0)

0:CommunityRiver1:CountyRiver2:StreetRiver3:Video4:WRY

    二、讀取向量資料

        讀取向量資料的步驟包括：

        1、開啟資料來源；

        2、擷取圖層；

        3、擷取圖層中的要素；

        4、擷取要素的屬性和幾何資訊。

        下面是讀取向量資料的程式碼片段：

import sysfrom osgeo import ogrimport ospybook as pbfn = r'D:\soft\geoprocessing-with-python\china_basic_map'ds = ogr.Open(fn, 0)if ds is None:    sys.exit('Could not open {0}.'.format(fn))'''索引擷取資料來源中的圖層方式'''lyr = ds.GetLayer(0)print(lyr.GetName())i = 0for fea in lyr:    pt = fea.geometry()    x = pt.GetX()    y = pt.GetY()    '''欄位名稱 擷取屬性值'''    code = fea.GetField('GBCODE')    '''對象方式擷取屬性值'''    length = fea.LENGTH    '''對象方式擷取屬性值2'''    lpoly_ = fea['LPOLY_']    '''索引方式擷取屬性值'''    fnode_ = fea.GetField(0)    '''擷取特定類型的屬性值'''    str_len = fea.GetFieldAsString('LENGTH')    print(code, length, fnode_, lpoly_, str_len, x, y)    i += 1    if i == 20:        break'''名稱擷取資料來源中的圖層'''lyr2 = ds.GetLayer('國家')print(lyr2.GetName())

    三、擷取資料中繼資料

        向量資料的中繼資料，包括資料範圍、幾何類型、空間參考及圖層對象的概要資訊等，下面的程式碼片段展示如何擷取這些中繼資料資訊：

import sysfrom osgeo import  ogrfn = r'D:\soft\geoprocessing-with-python\china_basic_map'ds = ogr.Open(fn, 0)if ds is None:    sys.exit('Could not open {0}.'.format(fn))lyr = ds.GetLayer(0)'''擷取圖層範圍'''extent = lyr.GetExtent()print(extent)'''(73.44696044921875, 135.08583068847656, 3.408477306365967, 53.557926177978516)''''''(min_x, max_x, min_y, max_y)''''''擷取圖層集合類型，傳回值為數字''''''1:點， 2：線， 3：面'''geom_type = lyr.GetGeomType()print(geom_type) #2print(geom_type == ogr.wkbPoint) #Falseprint(geom_type == ogr.wkbLineString) #Trueprint(geom_type == ogr.wkbPolygon) #False'''通過要素擷取幾何類型'''fea = lyr.GetFeature(0)print(fea.geometry().GetGeometryType()) #2print(fea.geometry().GetGeometryName()) #LINESTRING'''擷取圖層的空間參考'''print(lyr.GetSpatialRef())# GEOGCS["GCS_WGS_1984",#     DATUM["WGS_1984",#         SPHEROID["WGS_84",6378137.0,298.257223563]],#     PRIMEM["Greenwich",0.0],#     UNIT["Degree",0.0174532925199433],#     AUTHORITY["EPSG","4326"]]print(lyr.schema)'''擷取圖層屬性名稱及資料類型'''for field in lyr.schema:    print(field.name, field.GetTypeName())# FNODE_ Integer64# TNODE_ Integer64# LPOLY_ Integer64# RPOLY_ Integer64# LENGTH Real# BOU2_4M_ Integer64# BOU2_4M_ID Integer64# GBCODE Integer

    四、向量資料寫入

        1、向量資料寫入的思路如下：

        （1）以讀寫入模式開啟（建立）資料來源；

        （2）擷取待添加要素的圖層（或新建立圖層以添加要素）；

        （3）建立空要素，並為幾何對象和屬性賦值；

        （4）將建立的要素插入圖層中。

        下面是利用一個圖層中的要素建立一個新圖層的詳細代碼：

# -*- coding:utf-8 -*-import sysfrom osgeo import ogr'''建立一個圖層（根據一個圖層中的要素建立）'''fn = r'D:\soft\geoprocessing-with-python\china_basic_map''''以讀寫入模式開啟資料來源'''ds = ogr.Open(fn, 1)if ds is None:    sys.exit('Could not open {0}.'.format(fn))in_lyr = ds.GetLayer('省會城市')# lyr = ds.GetLayer('國家')## from ospybook.vectorplotter import VectorPlotter# vp = VectorPlotter(False)# vp.plot(lyr, fill=False)# vp.draw()'''若存在同名圖層，則先刪除再創鍵'''if ds.GetLayer('capital_city'):    ds.DeleteLayer('capital_city')out_lyr = ds.CreateLayer('capital_city',                         in_lyr.GetSpatialRef(),                         ogr.wkbPoint)'''根據元圖層欄位對象建立新圖層的欄位'''out_lyr.CreateFields(in_lyr.schema)'''根據圖層定義建立一個空的要素'''out_defn = out_lyr.GetLayerDefn()out_fea = ogr.Feature(out_defn)for in_fea in in_lyr:    geom = in_fea.geometry()    out_fea.SetGeometry(geom)    for i in range(in_fea.GetFieldCount()):        value = in_fea.GetField(i)        if i != 5:            out_fea.SetField(i, value)        else:            print(value)    out_lyr.CreateFeature(out_fea)

        2、建立新的資料來源

        建立新資料來源的關鍵在使用正確的驅動程式，每種驅動程式只處理操作一種類型的向量資料。有兩種方式擷取驅動程式：（1）從一個已經開啟的資料集中擷取，這將允許建立一個新的資料來源，它和已存在的資料來源向量資料格式一樣；（2）使用OGR中的GetDriverByName函數，傳遞給它驅動程式的簡稱，驅動程式的名稱在OGR網站有介紹。下面是程式碼片段：

''''''''''建立新的資料來源'''''''''import sysfrom osgeo import ogrfn = r'D:\soft\geoprocessing-with-python\china_basic_map''''擷取資料驅動的第一種方式：根據已有資料來源擷取'''ds = ogr.Open(fn, 0)if ds is None:    sys.exit('Could not open {0}.'.format(fn))driver = ds.GetDriver()print(driver.GetName())'''擷取資料來源的第二種方式：GetDriverByName'''json_driver = ogr.GetDriverByName('GeoJSON')print(json_driver.GetName())'''建立GeoJson資料來源，資料來源路徑應到具體檔案名稱'''json_fn = r'D:\soft\geoprocessing-with-python\china_basic_map\json_fn.json'json_ds = json_driver.CreateDataSource(json_fn)if json_ds is None:    sys.exit('Could not create {0}'.format(json_fn))print(json_ds)

        3、建立屬性欄位

        要將一個欄位添加到圖層中，需要一個包含欄位名稱、資料類型、欄位和精度等重要訊息的FieldDefn對象，下面是建立屬性欄位的程式碼片段：

'''UseExceptions'''ogr.UseExceptions()json_fn = r'D:\soft\geoprocessing-with-python\china_basic_map\json_fn4.json'json_driver = ogr.GetDriverByName('GeoJSON')print('start')try:    json_ds = json_driver.CreateDataSource(json_fn)    '''建立屬性資料欄位'''    lyr = json_ds.CreateLayer('layer')    coord_fld = ogr.FieldDefn('X', ogr.OFTReal)    coord_fld.SetWidth(8)    coord_fld.SetPrecision(3)    lyr.CreateField(coord_fld)    coord_fld.SetName('Y')    lyr.CreateField(coord_fld)    #建立要素    fea = ogr.Feature(lyr.GetLayerDefn())    fea.SetField('X', 12.34)    fea.SetFID(1)    lyr.CreateFeature(fea)    #結果同步至硬碟（儲存到檔案中）    json_ds.SyncToDisk()except RuntimeError as e:    print(e)print('end')

    五、更新現有資料

        在進行向量資料處理時，有時需要更新現有資料，而不是建立一個全新的資料集。是否能夠更新以及支援的編輯操作，取決於資料格式。常見的更新資料操作包括 改變圖層定義（屬性欄位）、要素添加、更新和刪除等。下面分別進行說明：

        1、改變圖層定義

        包括修改屬性欄位的定義、新增屬性欄位及刪除屬性欄位，具體程式碼片段如下：

import sysfrom osgeo import ogr'''''''更新現有資料'''''''fn = r'D:\soft\geoprocessing-with-python\china_basic_map'ds = ogr.Open(fn, 1)if ds is None:    sys.exit('Could not open {0}.'.format(fn))'''更改圖層定義(新增、刪除、修改欄位)'''lyr = ds.GetLayer(0)lyrdefn = lyr.GetLayerDefn()i = lyrdefn.GetFieldIndex('GBCODE')'''擷取欄位的資料類型'''fld_type = lyrdefn.GetFieldDefn(i).GetType()print(fld_type)'''定義新的屬性欄位'''fld_defn = ogr.FieldDefn('GBCODECODE', fld_type)fld_defn2 = ogr.FieldDefn('NEWFIELD4', ogr.OFTInteger)'''使用ALTER_NAME_FLAG更改屬性欄位名稱，必須保證欄位的資料類型一致'''lyr.AlterFieldDefn(i, fld_defn, ogr.ALTER_NAME_FLAG)'''建立屬性欄位'''lyr.CreateField(fld_defn2)'''刪除欄位(參數為欄位索引值)'''lyr.DeleteField(lyr.FindFieldIndex('NEWFIELD', 0))lyr.DeleteField(lyrdefn.GetFieldIndex('NEWFIELD2'))

        2、要素添加、更新和刪除

'''增加、刪除和更新要素'''fn = r'D:\soft\geoprocessing-with-python\china_basic_map'ds = ogr.Open(fn, 1)if ds is None:    sys.exit('Could not open {0}.'.format(fn))#更新要素（新增屬性欄位）lyr = ds.GetLayer(0)lyr.CreateField(ogr.FieldDefn('NEWFIELD', ogr.OFTInteger))n = 1for fea in lyr:    fea.SetField('NEWFIELD', n*n)    lyr.SetFeature(fea)    n += 1#刪除要素for fea in lyr:    print(fea.GetField('NEWFIELD'))    if fea.GetField('NEWFIELD') == 3179089:        lyr.DeleteFeature(fea.GetFID())print(lyr.GetFeatureCount())

        對資料進行更新操作後，一般需要對資料庫進行壓縮、重組或重新計算資料的空間範圍等，如下：

ds.ExecuteSQL('REPACK ' + lyr.GetName()) #壓縮資料庫ds.ExecuteSQL('VACUUM') #重組功能ds.ExecuteSQL('RECOMPUTE EXTENT ON ' + lyr.GetName()) #重新計算圖層的空間範圍print(lyr.GetExtent())