QGISとOracleで逆ジオコーディング
今回は、QGISのPythonコンソールとオラクルを使って逆ジオコーディングをしてみたいと思います。
逆ジオコーディングとは
Wikipediaによるとジオコーディングは、以下のように定義されています。
ジオコーディング(英語: geocoding)は、各種情報に対して、関連する地理座標(典型的には緯度・経度)を付加すること、およびこれに関する技術やソフトウェアをいう。付加された地理座標のことをジオコードと称する。
出典元 Wikipedia ジオコーディング
一般的にジオコーディングは、住所の名称やコードから該当する場所の座標(主に緯度経度)を求めることを指します。
一方で、逆ジオコーディングは緯度経度などの座標から住所などの場所情報を求めることを指します。
逆ジオコーディングを提供しているサイトを以下にご紹介します。
[blogcard url=”https://www.finds.jp/rgeocode/index.html.ja”]
上記サイトでは、以下のように緯度経度を指定してリクエストすると、
https://www.finds.jp/ws/rgeocode.php?XML&lat=35.678813&lon=139.7631761
以下のようにXMLなどの形式で該当場所の住所が返却されます。
<?xml version="1.0" encoding="UTF-8"?>
<rgeocode xmlns="http://finds.jp/ts" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:xsd="http://www.w3.org/2001/XMLSchema">
<status>200</status>
<result>
<prefecture>
<pcode>13</pcode>
<pname>東京都</pname>
</prefecture>
<municipality>
<mname>千代田区</mname>
<mcode>13101</mcode>
</municipality>
<local>
<section>丸の内二丁目</section>
<homenumber>6</homenumber>
<distance>11.063310233226</distance>
<latitude>35.678785</latitude>
<longitude>139.763058</longitude>
</local>
</result>
<argument>
<latitude>35.678813</latitude>
<longitude>139.7631761</longitude>
<localradius>500</localradius>
<localmax>1</localmax>
</argument>
<meta name="thanks" content="このサービスは 国土交通省 提供 国土数値情報(行政区域) を利用しています" />
<meta name="thanks" content="このサービスは 国土交通省 提供 街区レベル位置参照情報および大字・町丁目レベル位置参照情報 を利用しています" />
</rgeocode>
まとめると、下図の赤線(地図から住所)が逆ジオコーディング、青線がジオコーディングとなります。
作成するプログラムと使用データ
今回は、QGISのPythonコンソールとオラクルを使って逆ジオコーディングを行うプログラムを作成します。
イメージは以下となります。
図を説明すると、以下のようになります。
①QGIS上にOpenStreetMapを表示し、Pythonコンソール上のプログラムにて座標を取得
②Pythonコンソール上のプログラムはオラクル上の空間データベースへクエリを発行
③検索結果を画面に表示
使用するデータなどは、以下の通りです。
[table id=4 /]
プログラムの作成
では作成作業に入りましょう。作業手順は以下の通りとなります。
1.cx_oracleの導入
2.データの確認
3.Pythonコンソールの起動
4.プログラムの作成
5.動作確認
1.cx_oracleの導入
QGISに同梱されているPython3.7へオラクル接続用のプラグイン(cx_oracle)を導入します。
この作業は、ちょっとした注意が必要です。それは、QGIS配下に専用のPython環境が存在しているので、そちらにcx_oracleを導入しなければなりません。
まずは、OSGeo4W Shellを管理者として起動します。
Shellが立ち上がったら以下のように作業します。(Windowsの場合)
>cd D:\products\qgis\apps\Python37
>set path=%path%;D:\products\qgis\apps\Python37
>python -V
Python 3.7.0
>python -m pip install –upgrade pip setuptools
>python -m pip install cx_oracle
Collecting cx_oracle
Downloading https://files.pythonhosted.org/packages/37/79/2537adc3db3c82a5fe36f89bd8e1ad6a3d08b15260135f9f785d6b7f2ffe/cx_Oracle-7.3.0-cp37-cp37m-win32.whl (134kB)
|████████████████████████████████| 143kB 1.1MB/s
Installing collected packages: cx-oracle
Successfully installed cx-oracle-7.3.0
最初のコマンドで、QGISのPythonディレクトリ(以降、%QGIS_PY%)へ移動します。
%QGIS_PY%は、QGISのインストールディレクトリ(私の場合、 D:\products\qgis )の下のapps\Python37となります。
次に、PATH環境変数に%QGIS_PY%を追加しPyhhonのバージョンがQGIS同梱の3.7であることを確認します。
あとは、念のためpipをアップグレードして、cx_oracleをインストールします。
2.データの確認
今回は、2015年版国勢調査の小地域(町丁・字等別)データを使用します。
元データのダウンロード、オラクルへの取り込みと動作確認は、こちらをご確認ください。
3.Pythonコンソールの起動
まず、QGISを起動してOpenStreetMap(XYZ)を表示します。
次に、メニューから[プラグイン]、[Pythonコンソール]を選択し表示されたPythonコンソールから、アイコン”エディタの表示”を押下すると以下の画面となります。(画面右下がPthoonエディタ)
4.プログラムの作成
Pythonコンソール上に、下記プログラムを打ち込みます。
なお、本プログラムのキャンバスからの座標取得部分は下記URLを参考に作成しています。
https://gis.stackovernet.com/ja/q/68447
from qgis.gui import QgsMapToolEmitPoint
from qgis.gui import QgsMessageBar
from qgis.core import (QgsCoordinateReferenceSystem,
QgsCoordinateTransform,
QgsProject,
QgsPointXY)
import cx_Oracle as co
class PrintClickedPoint(QgsMapToolEmitPoint):
def __init__(self, canvas):
self.canvas = canvas
QgsMapToolEmitPoint.__init__(self, self.canvas)
def canvasPressEvent( self, e ):
point = self.toMapCoordinates(self.canvas.mouseLastXY())
# 座標変換器を作成し変換を実施
crsSrc = QgsCoordinateReferenceSystem(3857)
crsDst = QgsCoordinateReferenceSystem(4612)
trans = QgsCoordinateTransform(crsSrc, crsDst, QgsProject.instance())
dstPnt = trans.transform(point, QgsCoordinateTransform.ForwardTransform)
# Db接続
con = co.connect("hogeUser", "hogePwd", "hogeTns")
# カーソルの取得
cu = con.cursor()
# クエリの実行
cu.execute(
"SELECT city_name||s_name FROM kokusei k WHERE SDO_CONTAINS( k.geometry," + \
"sdo_geometry(2001, 4612, " + \
"sdo_point_type(" + str(dstPnt[0]) + "," + str(dstPnt[1]) + \
",NULL), NULL, NULL)) = 'TRUE'"
)
# フェッチ
rows = cu.fetchone()
print(rows[0])
# クローズ
cu.close()
print ('({:.4f}, {:.4f})'.format(dstPnt[0], dstPnt[1]))
# コンソール出力
iface.messageBar().pushMessage("住所", rows[0], level=Qgis.Info)
canvas_clicked = PrintClickedPoint( iface.mapCanvas() )
iface.mapCanvas().setMapTool( canvas_clicked )
プログラムの概要を説明します。
1~8行目は、必要なモジュールのインポートです。
44、45行目がこのスクリプトが起動された時に呼び出されるメイン処理になります。
44行目はmapcanvasの取得と自作クリックイベント処理(PrintClickedPoint)のインスタンス生成、45行目が作成したインスタンス(canvas_clicked)をQGISへ登録する処理になります。
9~41行目はクリックイベント処理(PrintClickedPointクラス)で、10~12行目はコンストラクタ、14行目以降はクリック時のイベント処理となります。
コンストラクタの12行目で、QgsMapToolEmitPointのコンストラクタを再度呼び出し、canvasを自クラスの変数に置き換え(フック)しています。
クリック時のイベント処理は、14~20行目で座標の変換(3857→4612)を行い、 22~37行目でオラクルへの問い合わせを行っています。(この部分はリソースの無駄遣いがありそうですので見直しの余地があります)
最後の41行目でQGISの画面(メッセージバー)へ表示しています。
5.動作確認
最後に動作確認をしましょう。
プログラムの起動は、エディタの上部”スクリプトの実行”ボタンを押下します。
その状態で、警視庁付近で地図上でマウスクリックすると、以下のようになるはずです。
多摩地区も確認してみましょう。
まとめ
ここまで、以下を学習しました。
- QGISにてPythonコンソールの起動
- キャンバスから取得した座標を任意の座標系へ変換する処理
- cx_Oracleの導入とPythonコンソールからオラクルへのリクエスト方法
GIS構築の実業務では、クリックした座標を起点に処理を組んでいくケースが多いので、今回の座標取得は勉強になりました。
※留意点
32ビットのQGISではオラクル接続時に下記2種類のエラーが発生しました。
ORA-28040: 一致する認証プロトコルがありません
ORA-01017: ユーザー名/パスワードが無効です。ログオンは拒否されました。
解決方法は、以下のURLに示されています。
http://www.intellilink.co.jp/article/column/ora-report20180925.html
私の環境では、下記2つの手順により解決しましたので、参考にしてください。
1.sqlnet.oraへの追記
sqlnet.oraに、”SQLNET.ALLOWED_LOGON_VERSION_SERVER=11”を追加し、リスナーの再起動を行います。
2.ユーザーのパスワード再作成
下記SQLにてPythonからログオンするユーザーのパスワードを変更(再設定)します。
SQL>alter user hogeUser identified by hogePwd;