隣接ブロックを考慮したジオハッシュによるリバースジオコーディング

これまでに作成した関数を組み合わせて,Yahoo!やGoogleなどWebサービスのAPIを用いない自前のリバースジオコーディング関数を作成します。
流れは以下のようになります。以前作成したリバースジオコーディング関数は隣接ブロックを考慮していないため精度に劣りますが,今回は隣接ブロックも考慮して精度を向上します。

  • 位置参照テーブルを用意する。
  • ジオハッシュ変換関数を作成する。
  • 位置参照テーブルにジオハッシュを追加する。
  • 隣接ブロックのジオハッシュ関数を作成する。
  • 緯度経度から住所を変換する関数を作成する。
  • アドインとして保存する。
  • 作成するツールから呼び出して使う。

これまでのエントリーを元に以下にまとめていきます。

位置参照テーブルを用意する

リバースジオコーディングの範囲を広げすぎると,動作が重くなるので,用途に応じて範囲を決定すべきと思います。
例えば自治体の施設点検ツールに活用するのであれば,その都道府県の範囲のみを対象とすれば良いでしょう。
以下の例では東京都を対象とします。

位置参照情報をダウンロードする

国土交通省の「位置参照情報ダウンロードサービス」から対象とする位置参照情報をダウンロードします。東京都全域を作成するので,まず「都道府県単位」を選択します。
rgca01
次の画面で都道府県からは対象とする「東京都」を選択し,街区レベル,大字・町丁目レベルの両方を使うので「すべて」を選択します。またデータ整備年度は最新の「平成24年」を選択します。
rgca02
次の画面でダウンロードするデータを選択します。前の画面で東京都を選択しているので,すべてにチェックします。利用規約に同意した後,ダウンロードして下さい。
rgca03

位置参照テーブルを作成する

ダウンロードした2つのファイルを解凍します。中のCSVファイルは同じファイル名なので,まず「大字・町丁目レベル」のデータから開きます。
rgca04
rgca06
適当な名前をつけて,Excelマクロ有効ブック(拡張子xlsm)で保存します。
rgca05
不要な列は削除して,都道府県名,市町村名,大字町丁目名を結合した住所データを作成します。
rgca07
新たに作成した住所の列は「値として貼り付け」し,最終的に住所・緯度・経度のみのデータとします。
rgca08
次にもう一つの「街区レベル」のCSVファイルを開きます。
rgca09
同様に番地までをつなげて住所データを作成します。
rgca10
rgca11
しかしながら,国土交通省のデータ量が十分でなく,複数の番地で同じ座標が登録されている場合があります。
rgca12
リバースジオコーディングでは一つの座標から一つの住所を返せば良いので,代表値として最初の行だけ残して重複は削除することとします。
rgca13
rgca14
こちらの「街区レベル」データを,最初の「大字・町丁目レベル」データに統合します。
rgca15
統合した位置参照情報の範囲をテーブルに変換します。
rgca16

(参考エントリー「ジオハッシュ関数でWebサービスのAPIを使わずに簡易に逆ジオコーディング」)

ジオハッシュ変換関数を作成する

ジオハッシュは階層的な空間データ構造で,任意の精度で表現でき(桁数が多いと精度が高い)近隣の2地点を表すコードは、似たような文字列から構成されるハッシュ値です。より詳しい説明はこのページに記載しました。
VBEの画面に移動し,標準モジュールを追加します。
rgca17
10進数の緯度・経度をジオハッシュに変換する関数 CnvGeohash,その中で用いる10進数の緯度・経度をそれぞれ2進数に変換する関数を標準モジュールに追加します。
コードの詳細はこちらのページに記載してあります。
rgca18

位置参照テーブルにジオハッシュを割り当てる

シートに戻り,D列に「ジオハッシュ」列を作成し,先ほど作成したジオハッシュ変換関数により,各住所の位置情報に対応するジオハッシュを求めます。
rgca19
東京都の位置参照情報は重複を除いて大字・町丁目レベルと街区レベルの合計で約21万件ありましたので,ジオハッシュの計算に約80秒かかりました。(MacBook Pro Retina Core-i7 2.3GHzの場合。1件あたり約0.0004秒)
計算が終わったら,再計算が起こらないように,値として貼り付けしておきます。
rgca20

隣接ブロックのジオハッシュを求める関数を作成する。

ジオハッシュは原則,似ている文字列は近くなるのですが,隣接ブロックを考慮する必要があります。その理由はこのページに記載しました。また,隣接ブロックの求め方はこのページに,そのコードはこのページに記載しました。
標準モジュールに,この隣接ブロックのジオハッシュを求める関数を作成します。(省略された部分を補う必要があります。各ブロックの隣接ブロックを調べて下さい。)
rgca21

最後にこれまでの関数を使って,リバースジオコーディングを行う関数を作ります。
標準モジュールに以下のプロシージャを追加します。

ErrExit1:

Exit2:

End Function
[/crayon]
ソースコード内に緑色のコメントで説明を入れていますが,以下にポイントを説明します。
所々,確認用に Debug.Print を入れていますが,実際の運用時は不要です。
〔14〜16行目〕
まず,リバースジオコーディングしたい緯度・経度をジオハッシュに変換し,その桁数を求めておきます。
〔20〜27行目のループ〕
変換したジオハッシュで始まる文字列が,用意した変換テーブル内に見つかるかフィルタをかけて調べます。見つからなかった場合,桁を短くして繰り返します。桁が短くなってから見つかっても,知りたい位置から遠くなってしまうので,6桁までとしています。
〔34〜51行目のループ〕
上記のループで知りたい位置が含まれるジオハッシュが変換テーブル内に見つかった場合,その周りの隣接ブロックにもっと近い住所がないか調べる必要があります。
見つかったジオハッシュを中心に周りの8ブロックのジオハッシュを求め,その文字列で始まるジオハッシュを変換テーブル内から探します。
〔43〜39行目のネストされたループ〕
各隣接ブロックで見つかったジオハッシュはあとで最も近い位置を調べるため,いったん配列に行番号を蓄積しておきます。配列数はその都度増やしていくためRedim Preserveを用います。
〔56〜58行目〕
みつかったジオハッシュが含まれる行の緯度経度と知りたい位置の距離を計算します。
住所は候補として格納しておきます。
〔60〜69行目のループ〕
複数見つかった場合は,順次比較して,当初の位置と最も近い位置を探します。
最も近い位置の住所を戻り値として返します。
〔76行目〕
30行目から飛んできますが,ジオハッシュを6桁まで減らしても見つからない場合は,エラー値を返します。

この関数のテストを実施してみます。都庁付近で緯度35.69,経度139.692(緑の矢印のところ)をGoogleMapでリバースジオコーディングを行うと,「東京都西新宿2丁目8-1」と出ました。(赤Aのバルーン)
rgca22
イミディエイトウィンドウでこの緯度・経度を作成した関数で変換してみます。
rgca23
Debug.Printにより確認用の結果が出ていますが,
1行目は緯度35.69,経度139.692のジオハッシュがxn774c0t7d0yであることを示しています。
2行目は20万件の住所変換テーブルからこのジオハッシュを含むジオハッシュを探したところ,xn774c0と7桁まで短くしたら見つかったことを示しています。
3行目〜11行目はこのジオハッシュを中心とした隣接ブロックのジオハッシュです。
rgca24
12行目の”8″はこれらの9つのジオハッシュを含むジオハッシュを位置参照テーブルから探したところ8件見つかったことを示しています。
8件のうち,最も近い住所データは「東京都新宿区西新宿二丁目」でした。
実は「東京都新宿区西新宿二丁目8」のデータもあったのですが,最も近いデータはたまたま「東京都新宿区西新宿二丁目」の方だったようです。
rgca25

なお,今回は東京都全域で約20万件の位置参照テーブルとしたため,変換に数秒かかります。このスピードで実用的でない場合は,位置参照テーブルを分割して,地域またはジオハッシュの最初の桁などでシートを分割して検索範囲を狭めるのが良いでしょう。
実際,私がある県で16地域に分割し位置参照テーブルを1万件程度以下までに抑えたところ,変換は瞬時でした。このときはExcel 2003に対応させるため,少なくとも6万5千件以下にする必要もあったためですが。
この関数を利用するファイルすべてに位置参照テーブルを入れると重くなるので,実際はアドインとして作成し,各マクロから呼び出して使うのが良いと思われます。

関連記事