オラクルからのポリゴンデータ抽出
今回は、オラクルのジオメトリカラム（SDO_GEOMETRY）をWKTテキストファイルに出力します。
この作業に至るまで
以前の記事にて整備したオラクル11XEの環境を使って、国勢調査の小地域ポリゴンデータをWKTファイルとして出力するPythonスクリプトを作成しようと思っていたところ、思わぬトラップが待ち受けていました。以降、順次説明します。
なお、WKT（Well Known Text）とは、図形を表現するテキストファイルの一種で、GeoJSONなどと同じようなものです。ただし、表現できるのは図形のみで属性を付与することはできません。
WKTの詳細は以下を確認してください。
https://ja.wikipedia.org/wiki/Well-known_text
テーブル
今回使用するテーブルを確認します。
- データはOracle Map Builderにて作成
- データは東京都全域（島しょ部を含む）のH27国勢調査小地域ポリゴンデータ
- テーブル名はAREA、ジオメトリカラム（SDO_GEOMETRY）の名称はgeometry
- レコード件数は6,010件
以下に、SQL*Plusにて表示したカラム構成を示します。
空間関数の動作確認
次に、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
——————————————————–
千代田区丸の内１丁目 0
中央区八重洲２丁目 142.799695
中央区八重洲１丁目 153.273322
千代田区丸の内２丁目 181.262214
中央区京橋１丁目 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つ目の項目値により決まる回転方向ルールを守らなければなりません。オラクルのリファレンスによると、以下のようになります。
回転ルールはWKTと同じです。
この回転方向が守られないと、幾何関数で結果が異常になったりしますので、注意してください。
これらの詳細は、以下のリファレンス確認してください。
実装
さて、いよいよ実装です。
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オブジェクトを生成、操作するようです。
環境整備
では、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()
プログラムの説明は以下となります。
[table id=27 column_widths=”10%|90%”/]
出力ファイルの検証
プログラムを実行すると、6,010行、12.6MBのWKTファイルができました。
ファイルの内容は以下のようになります。
完成したWKTファイルをOpenJUMPで表示すると、以下のようになります。
まとめ
ちょっと苦労しましたが、SDO_GEOMETRYオブジェクトのの内容が確認でき、PythonからSDO_GEOMETRYへのアクセスも確認できました。
正直、”select geometry.Get_Wkt() from area;”をファイル出力すれば終了だと思っていたので、思わぬ展開でした。