オラクルからのポリゴンデータ抽出

今回は、オラクルのジオメトリカラム(SDO_GEOMETRY)をWKTテキストファイルに出力します。

この作業に至るまで

以前の記事にて整備したオラクル11XEの環境を使って、国勢調査の小地域ポリゴンデータをWKTファイルとして出力するPythonスクリプトを作成しようと思っていたところ、思わぬトラップが待ち受けていました。以降、順次説明します。

なお、WKT(Well Known Text)とは、図形を表現するテキストファイルの一種で、GeoJSONなどと同じようなものです。ただし、表現できるのは図形のみで属性を付与することはできません。

Wikipediaより引用

WKTの詳細は以下を確認してください。

https://ja.wikipedia.org/wiki/Well-known_text

テーブル

今回使用するテーブルを確認します。

  • データはOracle Map Builderにて作成
  • データは東京都全域(島しょ部を含む)のH27国勢調査小地域ポリゴンデータ
  • テーブル名はAREA、ジオメトリカラム(SDO_GEOMETRY)の名称はgeometry
  • レコード件数は6,010件

以下に、SQL*Plusにて表示したカラム構成を示します。

AREAテーブルのカラム構成

空間関数の動作確認

次に、AREAテーブルに対して空間関数の動作確認をしてみます。
例えば、逆ジオコードに使う点を内包するポリゴンの抽出は問題なく動作します。

SQL> SELECT city_name||s_name FROM area k WHERE SDO_CONTAINS( k.geometry,sdo_geometry(2001, 4612, sdo_point_type(139.259829,35.631703,NULL), NULL, NULL)) = ‘TRUE’;

CITY_NAME||S_NAME
——————————————————————————–
八王子市高尾町

同様に、点に近いポリゴンの取得も問題なく動作します。

SQL> SELECT city_name||s_name ,sdo_nn_distance(1) dist FROM area g WHERE  SDO_NN(g.geometry,sdo_geometry(2001, 4612, sdo_point_type(139.7670944,35.6806302,NULL), NULL, NULL),’sdo_num_res=5′ ,1) = ‘TRUE’ order by dist;

CITY_NAME||S_NAME DIST
——————————————————–
千代田区丸の内1丁目 0
中央区八重洲2丁目 142.799695
中央区八重洲1丁目 153.273322
千代田区丸の内2丁目 181.262214
中央区京橋1丁目 306.666962

WKT出力関数はNG

しかし、以下のようにWKTを取得するSQLを2種類発行しましたが、いずれもNULLが返却されます。


SQL>select a.geometry.Get_Wkt() from area a where rownum <10;
A.GEOMETRY.GET_WKT()
——————————————————————-

9行が選択されました。

SQL>select SDO_UTIL.TO_WKTGEOMETRY(geometry) from area where rownum <10;
SDO_UTIL.TO_WKTGEOMETRY(GEOMETRY)
——————————————————————-

9行が選択されました。

この現象は、Oracle11XEに(データベースのインスタンスに)JavaVMが組み込まれていないことが原因のようです。具体的には、WKT出力処理(SDO_GEOMETRY.Get_WKTメソッド、SDO_UTIL.TO_WKTGEOMETRYメソッド)は、オラクルインスタンスに組み込まれたJavaVM上のアプリケーションとして実装されているようなのですが、XEでは肝心のJavaVMが組み込まれておらず動作しません。
この情報源は以下(英語)です。10XEについて書かれていますが、11XEも現象は同じです。

https://stackoverflow.com/questions/44832223/oracle-converting-sdo-geometry-to-wkt

なお、18XEは未確認です。

オラクルのジオメトリ型

こうなったら、オラクルのジオメトリカラムの座標を直接参照してWKTファイルを作成するしかなさそうです。
参照対象のジオメトリカラムの詳細を確認します。

SDO_GEOMETRY

まずは、オラクルのジオメトリ型(SDO_GEOMETRY)を確認します。

オラクルテーブル上に点やポリゴンを収容するカラムは、SDO_GEOMETRY型です。
SDO_GEOMETRY型は、VARCHAR2やNUMBER型(プリミティブな型と言ってもよいと思います)とは異なるオブジェクト型で、ユーザーが任意の型を定義することができます。SDO_GEOMETRY型は、オラクルが定義したオブジェクト型となります。
SQL*Plusを使ってSDO_GEOMETRY型を確認してみましょう。

SQL> desc SDO_GEOMETRY
名前         型
—————————– ————————————————————
SDO_GTYPE NUMBER
SDO_SRID NUMBER
SDO_POINT MDSYS.SDO_POINT_TYPE
SDO_ELEM_INFO MDSYS.SDO_ELEM_INFO_ARRAY
SDO_ORDINATES MDSYS.SDO_ORDINATE_ARRAY
METHOD
—————————–
MEMBER FUNCTION GET_GTYPE RETURNS NUMBER
MEMBER FUNCTION GET_DIMS RETURNS NUMBER
MEMBER FUNCTION GET_LRS_DIM RETURNS NUMBER
MEMBER FUNCTION GET_WKB RETURNS BLOB
MEMBER FUNCTION GET_WKT RETURNS CLOB

定義された型とメソッドが確認できます。
オラクルは、オブジェクト型を使って空間データを定義し、基本的にはストアードファンクションを使って座標変換や幾何計算などの空間関数を実現しています。
また、オラクルの場合はRDBMS上での空間データの扱い方標準(SQL/MM)が定義される前から製品化されているため(Oracle7のSpatial Data Option)SQL/MMで定義された”ST_XXXX”といった形の関数とは異なる関数体系となっています。

SDO_GEOMETRYの詳細は、以下を確認してください。

https://docs.oracle.com/cd/E16338_01/appdev.112/b72087/sdo_objrelschema.htm#i1004087

SDO_GTYPE、SDO_SRID、SDO_ELEM_INFO

次に、SDO_GEOMETRYの内部に存在する、SDO_GTYPE、SDO_SRID、SDO_ELEM_INFOを確認します、

SQL> select a.geometry.sdo_gtype sdo_gtype ,a.geometry.sdo_srid sdo_srid ,a.geometry.sdo_elem_info sdo_ellem_info from area a where rownum<2;
SDO_GTYPE  SDO_SRID  SDO_ELEM_INFO
—————–  ————–  ————————————————–
2003     4612     SDO_ELEM_INFO_ARRAY(1, 1003, 1)

SDO_GTYPEはデータ型を意味しており、”2003”は2次元のポリゴンという意味になります。
SDO_SRIDは、空間参照系(=EPSGコード)のことです。4612はJGD2000です。
SDO_ELEM_INFOは、少しわかりにくいのですが、3項目で1つの意味を持ちます。上記の場合の1,1003,1は、”SDO_ORDINATE_ARRAYの1点目から始まるポリゴン外枠”という意味になります。

SDO_ELEM_INFOが以下の場合は、”SDO_ORDINATE_ARRAYの1点目から176座標目(X、Yを別々にカウント)までのポリゴン外枠と、それ以降の中抜きポリゴン”を意味します。

[1, 1003, 1, 177, 2003, 1]

このように、SDO_ELEM_INFOにはSDO_ORDINATE_ARRAY(構成点座標点列)の解釈方法が書かれています。

SDO_ORDINATE_ARRAY

次に、SDO_ORDINATE_ARRAYを確認します。

SQL> select a.geometry.sdo_ordinates from area a where rownum<2;
GEOMETRY.SDO_ORDINATES
———————————————————————————————————
SDO_ORDINATE_ARRAY(139.739265, 35.6830937, 139.740285, 35.6832588, 139.740227, 35.6836813, 139.74021, 35.6837541, 139.740131, 35.6841094, 139.74012, 35.6841646,139.739546, 35.6841158, 139.739436, 35.6846023, 139.739377, 35.6848627, 139.73926, 35.6853755, 139.739256, 35.6853949, 139.739195, 35.6857349, 139.739111, 35.6858594, 139.738954, 35.6860923, 139.73891, 35.6860468, 139.738058, 35.6851638, 139.737979, 35.6851373, 139.737923, 35.6851269, 139.737265, 35.6849505, 139.737545, 35.6842969, 139.737726, 35.6839121, 139.738147, 35.6839315, 139.73815, 35.6838702, 139.738169, 35.6834578, 139.738195, 35.6829611, 139.738379, 35.6829874, 139.738767, 35.6830305, 139.739224, 35.6830885, 139.739265, 35.6830937)

こちらは、構成点座標点列ですが、前述のSDO_ELEM_INFOの2つ目の項目値により決まる回転方向ルールを守らなければなりません。オラクルのリファレンスによると、以下のようになります。

オラクルのリファレンス 2.2.4 SDO_ELEM_INFO より

回転ルールはWKTと同じです。
この回転方向が守られないと、幾何関数で結果が異常になったりしますので、注意してください。

これらの詳細は、以下のリファレンス確認してください。

https://docs.oracle.com/cd/E57425_01/121/SPATL/sdo_geometry-object-type.htm#GUID-683FF8C5-A773-4018-932D-2AF6EC8BC119

実装

さて、いよいよ実装です。
PythonからSDO_GEOMETRYのSDO_ORDINATE_ARRAYの座標を取得しながらWKTファイルを作成します。
今回作成したプログラムの条件は以下となります。

  • 変換対象はポリゴンのみ(WKTのPOLYGON)
  • 点(POINT)やマルチポリゴン(MULTIPOLYGON)はデータが存在しないため未考慮
  • SDO_ELEM_INFOに従って外枠と内枠のポリゴンを判定してWKT化

PythonからSDO_GEOMETRY型へのアクセス方法

調べたところ、Pythonのオラクル接続ライブラリのcx_oracleを使えば、ジオメトリを直接操作できそうです。

https://github.com/oracle/python-cx_Oracle/blob/master/samples/InsertGeometry.py

上記のGitHubのページから一部を抜粋すると、以下のようにPythonにてSDO_GEOMETRYオブジェクトを生成、操作するようです。

GitHubの oracle/python-cx_Oracle より抜粋

環境整備

では、Python側の環境を整備します。
pip(pip3)にてcx_oracleを導入します。

% pip3 install cx_oracle
Collecting cx_oracle
Downloading cx_Oracle-8.0.1.tar.gz (325 kB)
|████████████████████████████████| 325 kB 693 kB/s
Using legacy setup.py install for cx-oracle, since package ‘wheel’ is not installed.
Installing collected packages: cx-oracle
Running setup.py install for cx-oracle … done
Successfully installed cx-oracle-8.0.1

MacOSXにてPythonは3.9を使いましたが、cx_oracleのバージョンは8.0.1でした。

Pythonスクリプト

作成したPythonスクリプトは以下のようになります。

import cx_Oracle as co

if __name__ == '__main__':
    # Db接続
    co.init_oracle_client(lib_dir="/Users/takamotokeiji/products/oracle/instantclient_12_2")
    con = co.connect("user/pwd@moselle")
    # カーソルの取得
    cu = con.cursor()
    # クエリの実行
    cu.execute("select pref||city||s_area ,city_name||s_name ,a.geometry from area a")

    # ファイルのオープン
    wktfd = open("./outfile.wkt" ,'w' ,encoding='utf-8')

    # フェッチ
    rows = cu.fetchall()
    # レコードの取得と表示
    for row in rows:
        elinfo = row[2].SDO_ELEM_INFO.aslist()
        if len(elinfo) > 3 :
            # 座標リストの取得
            xylist = row[2].SDO_ORDINATES.aslist()
            print(xylist)
            xylist_s2 = [xylist[idx:idx + 2] for idx in range(0,len(xylist), 2)]
            # 要素情報の取得
            elinfo_s3 = [elinfo[idx:idx + 3] for idx in range(0,len(elinfo), 3)]
            elinfo_numlist = elinfo[3::3]
            for (elinfo3,elinfo_nump) in zip(elinfo_s3 ,elinfo_numlist):
                if elinfo3[1]==1003:
                    # 外枠ポリゴン:先頭かつ必ず1つの前提
                    outBuff = "POLYGON(("
                    for xy in xylist_s2[elinfo3[0]-1:int(elinfo_nump/2)]:
                        outBuff += str(str(round(xy[0],7))+" "+str(round(xy[1],7))) + ","
                    outBuff = outBuff[:-1]
                    outBuff += ")"
                elif elinfo3[1]==2003:
                    # 穴あき部分
                    outBuff += ",("
                    for xy in xylist_s2[int(elinfo3[0]/2):int(elinfo_nump/2)]:
                        outBuff += str(str(round(xy[0],7))+" "+str(round(xy[1],7))) + ","
                    outBuff = outBuff[:-1]
                    outBuff += ")"
            
            # 終端の設定、ファイル書き込み
            wktfd.write( outBuff + ")" + "\n")
        else:
            outBuff = "POLYGON(("
            # シングルポリゴン
            xylist = row[2].SDO_ORDINATES.aslist()
            xylist_s2 = [xylist[idx:idx + 2] for idx in range(0,len(xylist), 2)]
            for xy in xylist_s2:
                outBuff += str(str(round(xy[0],7))+" "+str(round(xy[1],7))) + ","
            # 終端の設定、ファイル書き込み
            outBuff = outBuff[:-1]
            wktfd.write( outBuff + "))" + "\n")
    # カーソルのクローズ
    cu.close()
    # ファイルのクローズ
    wktfd.close()

プログラムの説明は以下となります。

行数内容
5Mac+InstantClientの環境ではオラクルライブラリのパスを明示的に指定しないとランタイムエラーが発生します。
19SDO_ELEM_INFOをリストとして取得します。
20SDO_ELEM_INFOの構成数で単一ポリゴンか中抜けポリゴンかを判断
22(中抜けPGN)SDO_ORDINATESをリストとして取得します。
24(中抜けPGN)SDO_ORDINATESを2次元配列化します。
27(中抜けPGN)SDO_ELEM_INFOを3次元配列化します。
28(中抜けPGN)3次元配列化したSDO_ELEM_INFOを使ってLoop
31〜35(中抜けPGN)外枠ポリゴンのWKTテキスト出力
35〜43(中抜けPGN)穴あき部分のWKTテキスト出力
47〜55単一ポリゴンのWKT出力

出力ファイルの検証

プログラムを実行すると、6,010行、12.6MBのWKTファイルができました。
ファイルの内容は以下のようになります。

完成したWKTファイルをOpenJUMPで表示すると、以下のようになります。

まとめ

ちょっと苦労しましたが、SDO_GEOMETRYオブジェクトのの内容が確認でき、PythonからSDO_GEOMETRYへのアクセスも確認できました。

正直、”select geometry.Get_Wkt() from area;”をファイル出力すれば終了だと思っていたので、思わぬ展開でした。