写真のジオタグなど位置情報から住所へ変換する逆ジオコーディングについては，WebサービスのAPIを利用する方法が考えられます。しかしながら，非公開アプリケーションとして用いる場合は，ライセンス料がネックとなることもあります。
そこで，逆ジオコーディングを自前で構築する方法はないものかと，情報収集を行ったところ，国土交通省が提供している「位置情報参照情報ダウンロードサービス」が使えそうだということが分かりました。また，位置情報を特定するには「ジオハッシュ」が有効だということも分かりました。
さらにいろいろ調べましたが，Web上にはまだ詳しい情報が少なくそのまま適用できそうなものも公開されていなかったため，自分でジオハッシュ関数を作成してみることにしました。
ジオハッシュとは
Geohash（ジオハッシュ）は、Gustavo Niemeyerがgeohash.orgというWebサービスを作成中に発明した経緯度に基づくジオコーディング方法の一つである。パブリックドメインになっている。階層的な空間データ構造であり、空間を分割していくことによって表現する。
ジオハッシュは、任意の精度で表現できる、文字列の末尾を削っていくと徐々に精度が落ちる、といった特徴がある。
そのため、近隣の2地点を表すコードは、似たような文字列から構成されることが多い。同時に、より多くの文字列が一致すれば、当該2点がより近いことを表す。
利用例としては，地点を特定したり，例えばデータベースなどに地点情報を格納したり，ジオタグへの適用も提案されているそうです。
アルゴリズム
Wikipediaのページではジオハッシュから緯度経度へ復号（デコード）する考え方が載っていますので，これを逆に読み取り，緯度経度からジオハッシュにエンコードする方法を考えてみます。まず，緯度については−90°〜0（南緯）と0〜＋90°（北緯）の2つに分割した場合，対象位置情報の緯度がどちらに入るかで，2進数のビットを決めます。大区分（0〜90°＝北緯）に入る場合は”1″，小区分（−90〜0°＝南緯）に入る場合は”0″になります。例えば日本は北緯なので”1″になります。これを左側の1ビット目とします。
次に，0~90°を再度2分割し，小区分（0〜45°），大区分（45°〜90°）のどちらに入るかを見ます。例えば東京の北緯はだいたい35°なので，小区分に入ります。したがって，左から2ビット目は”0″となります。さらに二分割すると，3ビット目は”1″になります。
これを繰り返すことにより，東京（＝都庁の位置）の緯度35.68952は，101100101100001000101000010111と変換されました。
経度についても大区分（0〜180°＝東経），小区分（−180〜0°＝西経）からスタートし，同様に2進数に変換すると，東京の経度139.69170は11100011010101100001100100010001となりました。
次にこれらから1つの2進数に変換しますが，左から順に奇数ビットに経度，偶数ビットに緯度を混ぜながら結合します。桁が不足するときは0で補います。上記の場合は，
経度：11100011010101100001100100010001
緯度：101100101100001000101000010111
結合：1110110100001110011100100010110000000110110000100001001101010010
この結合したものを5ビットずつに分割し，下表の文字列に変換します。（base32エンコード）
| 2進数 | 10進数 | base32 |
| 00000 | 0 | 0 |
| 00001 | 1 | 1 |
| 00010 | 2 | 2 |
| 00011 | 3 | 3 |
| 00100 | 4 | 4 |
| 00101 | 5 | 5 |
| 00110 | 6 | 6 |
| 00111 | 7 | 7 |
| 01000 | 8 | 8 |
| 01001 | 9 | 9 |
| 01010 | 10 | b |
| 01011 | 11 | c |
| 01100 | 12 | d |
| 01101 | 13 | e |
| 01110 | 14 | f |
| 01111 | 15 | g |
| 10000 | 16 | h |
| 10001 | 17 | j |
| 10010 | 18 | k |
| 10011 | 19 | m |
| 10100 | 20 | n |
| 10101 | 21 | p |
| 10110 | 22 | q |
| 10111 | 23 | r |
| 11000 | 24 | s |
| 11001 | 25 | t |
| 11010 | 26 | u |
| 11011 | 27 | v |
| 11100 | 28 | w |
| 11101 | 29 | x |
| 11110 | 30 | y |
| 11111 | 31 | z |
先ほど結合した東京の位置情報の2進数を変換すると下表の様になりました。
| 11101 | 10100 | 00111 | 00111 | 00100 | 01011 | 00000 | 00110 | 11000 | 01000 | 01001 | 10101 | 0010 |
| x | n | 7 | 7 | 4 | c | 0 | 6 | s | 8 | 9 | p | 切り捨て |
コーディング
上記のようなアルゴリズムを確認しましたので，VBAで関数を作成してみます。
冗長な部分もあるかもしれませんが，可読性を優先しており，十分なスピードも得られています。
まずは，2進数への変換部分です。標準モジュールに以下のコードを作成します。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 |
Option Explicit '緯度を2進数へ変換 Function CnvLat10to2(ByVal Lat10 As Double) As String Dim b As Integer, L As Double If Lat10 > 0 Then CnvLat10to2 = "1" '北緯の場合1ビット目は"1" Else CnvLat10to2 = "0" '南緯の場合1ビット目は"0" Lat10 = Lat10 + 90 End If '十分な精度が得られるまで繰り返す Do Until Lat10 < 1e-07 'コンマ1秒程度の精度を考慮 b = b + 1 'ビットをインクリメント L = 90 / 2 ^ b '区分を順次半分へ If Lat10 < L Then '小区分の場合 CnvLat10to2 = CnvLat10to2 & "0" Else '大区分の場合 CnvLat10to2 = CnvLat10to2 & "1" Lat10 = Lat10 - L '評価した分を差し引く End If Loop End Function '経度を2進数へ変換 Function CnvLng10to2(ByVal Lng10 As Double) As String Dim b As Integer, L As Double If Lng10 > 0 Then CnvLng10to2 = "1" '東経の場合1ビット目は"1" Else CnvLng10to2 = "0" '西経の場合1ビット目は"0" Lng10 = Lng10 + 180 End If '十分な精度が得られるまで繰り返す Do Until Lng10 < 1e-07 'コンマ1秒程度の精度を考慮 b = b + 1 'ビットをインクリメント L = 180 / 2 ^ b '区分を順次半分へ If Lng10 < L Then '小区分の場合 CnvLng10to2 = CnvLng10to2 & "0" Else '大区分の場合 CnvLng10to2 = CnvLng10to2 & "1" Lng10 = Lng10 - L '評価した分を差し引く End If Loop End Function |
次に，ジオハッシュへ変換します。引数は10進数の緯度経度で，上の関数を呼び出して2進数に変換します。戻り値はジオハッシュの文字列となります。
|
54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 |
'緯度経度からジオハッシュへ変換 Function CnvGeoHash(ByVal Lat10 As Double, _ ByVal Lng10 As Double) As String Dim Lat2 As String, Lng2 As String Dim LonLat2 As String, Hash As String Dim Digits As Integer, i As Integer Lat2 = CnvLat10to2(Lat10) '緯度を2進数へ変換 Lng2 = CnvLng10to2(Lng10) '経度を2進数へ変換 '少ない方の桁に足りない分だけ"0"を付けて桁を揃える Digits = Len(Lat2) - Len(Lng2) '桁数の差 If Digits > 0 Then '経度の方が少ない場合 Lng2 = Lng2 & String(Digits, 48) '少ない桁数だけ"0"を追加 Else '緯度の方が少ない場合 Lat2 = Lat2 & String(-Digits, 48) End If Digits = Len(Lat2) '揃えた桁数 For i = 1 To Digits '奇数ビットに経度、偶数ビットに緯度を入れる LonLat2 = LonLat2 & Mid$(Lng2, i, 1) & Mid(Lat2, i, 1) Next i '5桁ごとにBase32でエンコード Digits = Len(LonLat2) For i = 1 To Digits Step 5 Hash = Mid$(LonLat2, i, 5) '順次5桁毎に取り出す Select Case Hash '各桁のハッシュ値に変換 Case "00000": Hash = "0" Case "00001": Hash = "1" Case "00010": Hash = "2" Case "00011": Hash = "3" Case "00100": Hash = "4" Case "00101": Hash = "5" Case "00110": Hash = "6" Case "00111": Hash = "7" Case "01000": Hash = "8" Case "01001": Hash = "9" Case "01010": Hash = "b" Case "01011": Hash = "c" Case "01100": Hash = "d" Case "01101": Hash = "e" Case "01110": Hash = "f" Case "01111": Hash = "g" Case "10000": Hash = "h" Case "10001": Hash = "j" Case "10010": Hash = "k" Case "10011": Hash = "m" Case "10100": Hash = "n" Case "10101": Hash = "p" Case "10110": Hash = "q" Case "10111": Hash = "r" Case "11000": Hash = "s" Case "11001": Hash = "t" Case "11010": Hash = "u" Case "11011": Hash = "v" Case "11100": Hash = "w" Case "11101": Hash = "x" Case "11110": Hash = "y" Case "11111": Hash = "z" Case Else: Hash = "" '5桁未満の場合は切り捨て End Select CnvGeoHash = CnvGeoHash & Hash '各桁のハッシュ値を結合 Next End Function |
これでジオハッシュ関数は完成です。
次のプロシージャでテストしてみます。
|
128 129 130 131 132 133 |
Sub test() Debug.Print CnvGeoHash(35.68952, 139.6917) End Sub |
変換結果がイミディエイトウィンドウに表示されました。
念のため，下記のサイトで確認してみると，同じ結果が得られています。
