QGISとPostGISで逆ジオコーディング

2020年3月15日

今回は、QGISのPythonコンソールとPostGISを使って逆ジオコーディングをしてみたいと思います。

逆ジオコーディングとは

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コンソールとPostGISを使って逆ジオコーディングを行うプログラムを作成します。
イメージは以下となります。

図を説明すると、以下のようになります。
①QGIS上にOpenStreetMapを表示し、Pythonコンソール上のプログラムにて座標を取得
②Pythonコンソール上のプログラムはPostGIS上の空間データベースへクエリを発行
③検索結果を画面に表示

使用するデータなどは、以下の通りです。

[table id=5 /]

プログラムの作成

ではプログラム作成に入りましょう。作業手順は以下の通りとなります。
1.psycopg2の導入
2.データの確認
3.Pythonコンソールの起動
4.プログラムの作成
5.動作確認

1.psycopg2の導入
QGISに同梱されているPython3.7へPostgreSQL接続用のプラグイン（psycopg2）を導入します。
この作業は、ちょっとした注意が必要です。それは、QGIS配下に専用のPython環境が存在しているので、そちらにpsycopg2を導入しなければなりません。
まずは、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
　Successfully installed pip-20.0.2 setuptools-46.0.0
>python -m pip install psycopg2
　Collecting psycopg2
　Downloading psycopg2-2.8.4-cp37-cp37m-win32.whl (983 kB)
|████████████████████████████████| 983 kB 731 kB/s
　Installing collected packages: psycopg2
　Successfully installed psycopg2-2.8.4

最初のコマンドで、QGISのPythonディレクトリ（以降、%QGIS_PY%）へ移動します。
%QGIS_PY%は、QGISのインストールディレクトリ（私の場合、 D:\products\qgis ）の下のapps\Python37となります。
次に、PATH環境変数に%QGIS_PY%を追加しPyhhonのバージョンがQGIS同梱の3.7であることを確認します。
あとは、念のためpipをアップグレードして、psycopg2をインストールします。

2.データの確認
今回は、2015年版国勢調査の小地域（町丁・字等別）データを使用します。
元データのダウンロード、PostGISへの取り込みと動作確認は、こちらをご確認ください。

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 psycopg2

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 = psycopg2.connect(host='sakura', port=5432,database='postgis', user='takamoto',password='hoge')

        # カーソルの取得
        cu = con.cursor()
        # クエリの実行
        cu.execute(
            "SELECT city_name||s_name as name FROM kokusei " + \
            "WHERE ST_Contains( geom ,ST_GeometryFromText('POINT(" + str(dstPnt[0]) + " " + str(dstPnt[1]) + ")',4612))"
        )
        # フェッチ
        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～7行目は、必要なモジュールのインポートです。
44、45行目がこのスクリプトが起動された時に呼び出されるメイン処理になります。
44行目はmapcanvasの取得と自作クリックイベント処理（PrintClickedPoint）のインスタンス生成、45行目が作成したインスタンス（canvas_clicked）をQGISへ登録する処理になります。
9～41行目はクリックイベント処理（PrintClickedPointクラス）で、10～12行目はコンストラクタ、14行目以降はクリック時のイベント処理となります。
コンストラクタの12行目で、QgsMapToolEmitPointのコンストラクタを再度呼び出し、canvasを自クラスの変数に置き換え（フック）しています。
クリック時のイベント処理は、18～21行目で座標の変換（3857→4612）を行い、 24～37行目でPostGISへの問い合わせを行っています。（この部分はリソースの無駄遣いがありそうですので見直しの余地があります）
最後の41行目でQGISの画面（メッセージバー）へ表示しています。

5.動作確認
最後に動作確認をしましょう。
プログラムの起動は、エディタの上部”スクリプトの実行”ボタンを押下します。
その状態で、警視庁付近で地図上でマウスクリックすると、以下のようになるはずです。