Python+GDAL/OGRベクトルデータ読み書き

一般的なベクトルデータフォーマットは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.county
1:tiger.state
2:tiger.place
3:tiger.cousub
4:tiger.edges
5:tiger.addrfeat
6:tiger.faces
7:tiger.zcta5
8:tiger.tract
9:tiger.tabblock
10:tiger.bg

3、フォルダデータソース(shapefilesとCSV)を開く
shapefilesフォルダを開き、次のように出力します.

#  shapefile   
    shp_fn = r'E:\HZZ_GXDP_Spatial_Data'
    print_layers(shp_fn, 0)

0:CommunityRiver
1:CountyRiver
2:StreetRiver
3:Video
4:WRY

二、ベクトルデータを読み取る
ベクトルデータを読み込むには、次の手順に従います.
1、データソースを開く;
2、レイヤーを取得する;
3、レイヤーの要素を取得する;
4、フィーチャのプロパティとジオメトリ情報を取得します.
次に、ベクトルデータを読み込むコードクリップを示します.

import sys
from osgeo import ogr
import ospybook as pb

fn = 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 = 0
for 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 sys
from osgeo import  ogr

fn = 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) #2
print(geom_type == ogr.wkbPoint) #False
print(geom_type == ogr.wkbLineString) #True
print(geom_type == ogr.wkbPolygon) #False

'''          '''
fea = lyr.GetFeature(0)
print(fea.geometry().GetGeometryType()) #2
print(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 sys
from 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つのタイプのベクトルデータのみを処理します.ドライバを取得するには、(1)開いているデータセットから取得します.これにより、既存のデータソースベクトルデータフォーマットと同じ新しいデータソースを作成できます.(2)OGRのGetDriverByName関数を使用して、そのドライバの略称を渡します.ドライバの名前はOGRのウェブサイトで紹介されています.次はコードクリップです.

''''''''''       '''''''''
import sys
from osgeo import ogr

fn = 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 sys
from 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 = 1
for 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())