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

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

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

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

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

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

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

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

国土交通省の「位置参照情報ダウンロードサービス」から対象とする位置参照情報をダウンロードします。東京都全域を作成するので，まず「都道府県単位」を選択します。

次の画面で都道府県からは対象とする「東京都」を選択し，街区レベル，大字・町丁目レベルの両方を使うので「すべて」を選択します。またデータ整備年度は最新の「平成24年」を選択します。

次の画面でダウンロードするデータを選択します。前の画面で東京都を選択しているので，すべてにチェックします。利用規約に同意した後，ダウンロードして下さい。

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

ダウンロードした2つのファイルを解凍します。中のCSVファイルは同じファイル名なので，まず「大字・町丁目レベル」のデータから開きます。

適当な名前をつけて，Excelマクロ有効ブック（拡張子xlsm）で保存します。

不要な列は削除して，都道府県名，市町村名，大字町丁目名を結合した住所データを作成します。

新たに作成した住所の列は「値として貼り付け」し，最終的に住所・緯度・経度のみのデータとします。

次にもう一つの「街区レベル」のCSVファイルを開きます。

同様に番地までをつなげて住所データを作成します。

しかしながら，国土交通省のデータ量が十分でなく，複数の番地で同じ座標が登録されている場合があります。

リバースジオコーディングでは一つの座標から一つの住所を返せば良いので，代表値として最初の行だけ残して重複は削除することとします。

こちらの「街区レベル」データを，最初の「大字・町丁目レベル」データに統合します。

統合した位置参照情報の範囲をテーブルに変換します。

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

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

ジオハッシュは階層的な空間データ構造で，任意の精度で表現でき（桁数が多いと精度が高い）近隣の2地点を表すコードは、似たような文字列から構成されるハッシュ値です。より詳しい説明はこのページに記載しました。
VBEの画面に移動し，標準モジュールを追加します。

10進数の緯度・経度をジオハッシュに変換する関数 CnvGeohash，その中で用いる10進数の緯度・経度をそれぞれ2進数に変換する関数を標準モジュールに追加します。
コードの詳細はこちらのページに記載してあります。

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

シートに戻り，D列に「ジオハッシュ」列を作成し，先ほど作成したジオハッシュ変換関数により，各住所の位置情報に対応するジオハッシュを求めます。

東京都の位置参照情報は重複を除いて大字・町丁目レベルと街区レベルの合計で約21万件ありましたので，ジオハッシュの計算に約80秒かかりました。（MacBook Pro Retina Core-i7 2.3GHzの場合。1件あたり約0.0004秒）
計算が終わったら，再計算が起こらないように，値として貼り付けしておきます。

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

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

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

'隣接ブロックを考慮したジオハッシュによるリバースジオコーディング
Public Function ReverseGeoCoding(ByVal Lat As Double, ByVal Lng As Double) As String

[crayon-5d319f67d6114598233561 ]<code>Dim Geohash As String
Dim Dist As Double, Dist1 As Double
Dim i As Integer, Digits As Integer
Dim R() As Long
Dim c As Variant
Dim k As Long
Dim NBhash As String

Application.ScreenUpdating = False

Geohash = CnvGeoHash(Lat, Lng) '緯度経度をジオハッシュに変換
Debug.Print Geohash
Digits = Len(Geohash) '最初のジオハッシュの桁数

With ThisWorkbook.Worksheets("Sheet1").ListObjects("テーブル1")
    With .ListColumns("ジオハッシュ")
        For i = Digits To 6 Step -1 'ジオハッシュの最大桁数から6桁まで
            '桁数を短くしていく
            Geohash = Left(Geohash, i)
            'そのジオハッシュで始まるジオハッシュで絞り込み
            .Parent.Range.AutoFilter field:=.Index, Criteria1:=Geohash &amp; "*"
            '見つかった場合ループを抜け出す
            If .Range.SpecialCells(xlCellTypeVisible).Count &gt; 1 Then Exit For
        Next

        '見つからなかった場合エラー値を返して終了
        If i = 5 Then GoTo ErrExit1

        Debug.Print Geohash

        For i = 0 To 8 '見つかった場合，そのジオハッシュと隣接8ブロックに対して
            '順次絞り込みを行う
            NBhash = NeighborBlock(Geohash, i)
            Debug.Print NBhash
            .Parent.Range.AutoFilter _
                field:=.Index, Criteria1:=NBhash &amp; "*"
            '各ジオハッシュ内に住所が見つかった場合
            If .Range.SpecialCells(xlCellTypeVisible).Count &gt; 1 Then
                '見つかったセルを順次調べる
                For Each c In .Range.SpecialCells(xlCellTypeVisible)
                    If c.Row &gt; 1 Then '1行目はテーブル見出しなので無視する
                        k = k + 1
                        ReDim Preserve R(k) '配列を見つかった数まで拡張する
                        R(k) = c.Row '配列に見つかった住所の行番号を格納する
                    End If
                Next
            End If
        Next
    End With

    Debug.Print k
    '少なくとも1つは見つかっているので，1つめの位置と当初の位置の距離を求める。
    Dist = ((Lat - .ListColumns("緯度").Range(R(1))) ^ 2 + _
            (Lng - .ListColumns("経度").Range(R(1))) ^ 2) ^ 0.5
    ReverseGeoCoding = .ListColumns("住所").Range(R(1)) '1つめの住所を返す

    For i = 2 To k '複数見つかった場合
        '各位置と当初の位置の距離を求める。
        Dist1 = ((Lat - .ListColumns("緯度").Range(R(i))) ^ 2 + _
                 (Lng - .ListColumns("経度").Range(R(i))) ^ 2) ^ 0.5
        '1つめの位置よりも近い場合
        If Dist1 &lt; Dist Then
            Dist = Dist1 'それを採用する
            ReverseGeoCoding = .ListColumns("住所").Range(R(i)) '住所を更新する。
        End If
    Next
End With

GoTo Exit2

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

Public Function ReverseGeoCoding(ByVal Lat As Double, ByVal Lng As Double) As String

[crayon-5d319f67d6114598233561 ]<code>Dim Geohash As String

Dim Dist As Double, Dist1 As Double

Dim i As Integer, Digits As Integer

Dim R() As Long

Dim c As Variant

Dim k As Long

Dim NBhash As String

Application.ScreenUpdating = False

Geohash = CnvGeoHash(Lat, Lng) '緯度経度をジオハッシュに変換

Debug.Print Geohash

Digits = Len(Geohash) '最初のジオハッシュの桁数

With ThisWorkbook.Worksheets("Sheet1").ListObjects("テーブル1")

With .ListColumns("ジオハッシュ")

For i = Digits To 6 Step -1 'ジオハッシュの最大桁数から6桁まで

'桁数を短くしていく

Geohash = Left(Geohash, i)

'そのジオハッシュで始まるジオハッシュで絞り込み

.Parent.Range.AutoFilter field:=.Index, Criteria1:=Geohash & "*"

'見つかった場合ループを抜け出す

If .Range.SpecialCells(xlCellTypeVisible).Count > 1 Then Exit For

'見つからなかった場合エラー値を返して終了

If i = 5 Then GoTo ErrExit1

Debug.Print Geohash

For i = 0 To 8 '見つかった場合，そのジオハッシュと隣接8ブロックに対して

'順次絞り込みを行う

NBhash = NeighborBlock(Geohash, i)

Debug.Print NBhash

.Parent.Range.AutoFilter _

field:=.Index, Criteria1:=NBhash & "*"

'各ジオハッシュ内に住所が見つかった場合

If .Range.SpecialCells(xlCellTypeVisible).Count > 1 Then

'見つかったセルを順次調べる

For Each c In .Range.SpecialCells(xlCellTypeVisible)

If c.Row > 1 Then '1行目はテーブル見出しなので無視する

k = k + 1

ReDim Preserve R(k) '配列を見つかった数まで拡張する

R(k) = c.Row '配列に見つかった住所の行番号を格納する

End If

End With

Debug.Print k

'少なくとも1つは見つかっているので，1つめの位置と当初の位置の距離を求める。

Dist = ((Lat - .ListColumns("緯度").Range(R(1))) ^ 2 + _

(Lng - .ListColumns("経度").Range(R(1))) ^ 2) ^ 0.5

ReverseGeoCoding = .ListColumns("住所").Range(R(1)) '1つめの住所を返す

For i = 2 To k '複数見つかった場合

'各位置と当初の位置の距離を求める。

Dist1 = ((Lat - .ListColumns("緯度").Range(R(i))) ^ 2 + _

(Lng - .ListColumns("経度").Range(R(i))) ^ 2) ^ 0.5

'1つめの位置よりも近い場合

If Dist1 < Dist Then

Dist = Dist1 'それを採用する

ReverseGeoCoding = .ListColumns("住所").Range(R(i)) '住所を更新する。

End If

End With

GoTo Exit2

ErrExit1:

ReverseGeoCoding = "#N/A" 'エラー値を返す

1 2	ReverseGeoCoding = "#N/A" 'エラー値を返す

Exit2:

Application.ScreenUpdating = True

1 2	Application.ScreenUpdating = True

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でリバースジオコーディングを行うと，「東京都西新宿２丁目８-１」と出ました。（赤Aのバルーン）

イミディエイトウィンドウでこの緯度・経度を作成した関数で変換してみます。

Debug.Printにより確認用の結果が出ていますが，
1行目は緯度35.69，経度139.692のジオハッシュがxn774c0t7d0yであることを示しています。
2行目は20万件の住所変換テーブルからこのジオハッシュを含むジオハッシュを探したところ，xn774c0と7桁まで短くしたら見つかったことを示しています。
3行目〜11行目はこのジオハッシュを中心とした隣接ブロックのジオハッシュです。

12行目の”8″はこれらの9つのジオハッシュを含むジオハッシュを位置参照テーブルから探したところ8件見つかったことを示しています。
8件のうち，最も近い住所データは「東京都新宿区西新宿二丁目」でした。
実は「東京都新宿区西新宿二丁目8」のデータもあったのですが，最も近いデータはたまたま「東京都新宿区西新宿二丁目」の方だったようです。

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