【Python と GIS】GDAL を使ってファイルジオデータベースを操作してみよう！

さて、本日も GDAL について書いてみようと思います。前回はGDAL で Shapeファイルの読み込みと属性検索について書きましたが、今回はGDAL でのファイルジオデータベースの操作について書いてみようと思います。興味のある方は前回のエントリーもぜひ読んでみてください。

www.gis-py.com

GDALとは

GDAL はフリーの GIS ライブラリです。GDAL — GDAL documentation
これがあればShapeの読み込みからジオメトリの作成まで基本的なことは何でもできるかと思います。

ファイルジオデータベースとは

以下エントリーでも紹介していますが、ファイルジオデータベースとは ArcGIS で使用されるデータフォーマットです。

www.gis-py.com

インストール

公式サイトの手順に従ってインストールをお願いします。https://gdal.org/download.html

実行環境

Windows10(64bit)
Python3.6.5
Jupyter Notebook

サンプルデータ

ArcMap で以下を作成しました。

GDB・・・Sample.gdb
フィーチャクラス・・・GDAL_TEST
カラム・・・テスト1、テスト2

f:id:sanvarie:20190713163137p:plain

f:id:sanvarie:20190713163224p:plain

サンプルコード

GDAL_TEST フィーチャクラスのテスト1、テスト2カラムの値を取得するサンプルです。

from osgeo import ogr

# GDBオープン
driver = ogr.GetDriverByName("OpenFileGDB")
gdb = driver.Open(r"D:\python\data\Sample.gdb", 0)

list1, list2 = [], []

# GDBに格納されているフィーチャクラス数分ループ
for featsClass_idx in range(gdb.GetLayerCount()):

    # フィーチャクラスを取得
    featureClass = gdb.GetLayerByIndex(featsClass_idx)

    layerDefinition = featureClass.GetLayerDefn()

    # GDAL_TESTフィーチャクラスの場合のみ処理を実行
    if featureClass.GetName() == "GDAL_TEST":

        while featureClass:

            feat = featureClass.GetNextFeature()

            if feat is None:
                break

            else:
                for column_num in range(layerDefinition.GetFieldCount()):

                    # カラム名取得
                    column_name = layerDefinition.GetFieldDefn(column_num).GetName()

                    #テスト1、テスト2カラムのときのみ属性を取得する
                    if column_name == "テスト1":
                        list1.append(feat.GetField(column_name))
                    elif column_name == "テスト2":
                        list2.append(feat.GetField(column_name))

print(list1)
print(list2)

del gdb

結果を見ると、想定通りテスト1、テスト2カラムの値を取得できたことがわかりました。
f:id:sanvarie:20190713163525p:plain

処理のポイント

以下にポイントとなる箇所を記述します。

GDBオープン

# GDBオープン
driver = ogr.GetDriverByName("OpenFileGDB")
gdb = driver.Open(r"D:\python\data\Sample.gdb", 0)

フィーチャクラス取得

GetLayerByIndex(index) を使用することでフィーチャクラスを取得することができます。

# GDBに格納されているフィーチャクラス数分ループ
for featsClass_idx in range(gdb.GetLayerCount()):

    # フィーチャクラスを取得
    featureClass = gdb.GetLayerByIndex(featsClass_idx)

フィーチャクラスの走査

GetNextFeature() を使用することで、フィーチャクラス内のフィーチャにアクセスできます。

while featureClass:
    feat = featureClass.GetNextFeature()

属性の取得

基本的には前回のエントリー(https://www.gis-py.com/entry/gdal-shape)と同じですが、今回はリスト内包表記は使用していません。ぱっと見はこちらの方がわかりやすいかもしれないですね。

for column_num in range(layerDefinition.GetFieldCount()):

    # カラム名取得
    column_name = layerDefinition.GetFieldDefn(column_num).GetName()

    #テスト1、テスト2カラムのときのみ属性を取得する
    if column_name == "テスト1":
        list1.append(feat.GetField(column_name))
    elif column_name == "テスト2":
        list2.append(feat.GetField(column_name))

本日は以上になります。また GDAL について書いてみようと思いますが、今度はもう少し突っ込んだ内容にしてみようかと思います。

GIS奮闘記

元GISエンジニアの技術紹介ブログ。主にPythonを使用。