Overview

緯度経度から符号化方式の Quadkey に変換して、近傍検索を試してみます。
実践編なので予備知識などは別途ググってください。

Quadkey is 何？

ざっくり言うと、領域を再帰的に4分割して絞り込んでいく手法で、桁1つは4進数表記です。
4分木をそのまま落とし込んだもので、2次元のZ階数曲線と同じものです。

詳しい説明はこの辺見てください。

可視化

GeoJson の可視化サイトが便利なので、確認にはあとでこちらを使います。
https://www.townsendjennings.com/geo

データを取ってくる

とりあえずサンプル代わりに、
国土交通省の位置参照情報ダウンロードサービスから、平成29年度の東京都の 市区町村全域 -> 街区のデータを取得します。

中身

例によって SJIS なので UTF-8 に変換します。

こんな感じで入ってるはず

$ nkf -w 13_2017.csv | head
"都道府県名","市区町村名","大字・丁目名","小字・通称名","街区符号・地番","座標系番号","Ｘ座標","Ｙ座標","緯度","経度","住居表示フラグ","代表フラグ","更新前履歴フラグ","更新後履歴フラグ"
"東京都","千代田区","麹町六丁目","","1","9","-34981.9","-9228.7","35.684649","139.731373","0","1","1","0"
"東京都","千代田区","麹町六丁目","","5","9","-34981.9","-9228.7","35.684649","139.731373","0","1","1","0"
"東京都","千代田区","六番町","","6","9","-34653.0","-9100.3","35.687614","139.732787","0","1","1","0"
"東京都","千代田区","六番町","","8","9","-34653.0","-9100.3","35.687614","139.732787","0","1","1","0"
"東京都","千代田区","六番町","","10","9","-34653.0","-9100.3","35.687614","139.732787","0","1","1","0"
"東京都","千代田区","六番町","","12","9","-34653.0","-9100.3","35.687614","139.732787","0","1","1","0"
"東京都","千代田区","六番町","","14","9","-34653.0","-9100.3","35.687614","139.732787","0","1","1","0"
"東京都","千代田区","麹町六丁目","","1","9","-34968.4","-9134.2","35.684770","139.732417","0","0","1","0"
"東京都","千代田区","紀尾井町","","5","9","-35262.7","-9173.3","35.682118","139.731988","1","1","1","0"

DB に入れる

このままだと扱いづらいので sqlite3 に入れます。

$ sqlite3 --version
3.32.2 2020-06-04 12:58:43 ec02243ea6ce33b090870ae55ab8aa2534b54d216d45c4aa2fdbb00e86861e8c

投入

$ nkf -w 13_2017.csv | sqlite3 -csv -header a.db ".import /dev/stdin gis"

確認

$ sqlite3 -header a.db "select * from gis limit 3"
都道府県名|市区町村名|大字・丁目名|小字・通称名|街区符号・地番|座標系番号|Ｘ座標|Ｙ座標|緯度|経度|住居 表示フラグ|代表フラグ|更新前履歴フラグ|更新後履歴フラグ
東京都|千代田区|麹町六丁目||1|9|-34981.9|-9228.7|35.684649|139.731373|0|1|1|0
東京都|千代田区|麹町六丁目||5|9|-34981.9|-9228.7|35.684649|139.731373|0|1|1|0
東京都|千代田区|六番町||6|9|-34653.0|-9100.3|35.687614|139.732787|0|1|1|0

緯度経度から Quadkey を算出

なんかのライブラリを使ってもいいんですけど、勉強がてらに書いたので今回はこれを使います。

$ g++ quadkey.cpp -o quadkey

投入

ズームレベル 20 で算出してみます。

# !/bin/bash

F1=$(mktemp)
F2=$(mktemp)

echo 'select distinct "緯度","経度" from gis' \
    | sqlite3 -separator $'\t' a.db \
    | tee ${F2} | ./quadkey -z 20 >${F1}

{
    printf "%s\t%s\t%s\n" lat lon quadkey
    paste ${F2} ${F1}
} | sqlite3 -separator $'\t' a.db ".import /dev/stdin quadkey"

確認

$ sqlite3 -header a.db "select * from quadkey limit 3"
lat|lon|quadkey
35.684649|139.731373|13300211230133032303
35.687614|139.732787|13300211230133031023
35.684770|139.732417|13300211230133033202

ここまで準備。

近傍検索してみる

試しに JR新宿駅の周辺を調べてみます。
こちらによると、JR新宿駅の緯度経度は

| 駅名 | 緯度 | 経度 |
| 新宿 | 35.690921 | 139.70025799999996 |

だそうなので、quadkey に変換すると

$ echo 35.690921 139.70025799999996 | ./quadkey
13300211230123111232

可視化してみると

zoom level 20

quadkey は↑で算出した 13300211230123111232 です。

狭ｗ
ズームレベル20だと狭いので、もう少し範囲を広げてみます。

zoom level 17

quadkey は 13300211230123111 (末尾から3桁落とします)。

多少広くなりました。

この範囲内の住所を近傍検索してみます。

quadkey で前方一致検索

さて、ここで再び sql に戻ってきます。
上記の zoom level 17 の 1330021123012311 で検索してみると、

a.sql

SELECT
  gis."都道府県名"
  , gis."市区町村名"
  , gis."大字・丁目名"
  , gis."小字・通称名"
  , gis."街区符号・地番"
  , gis."緯度"
  , gis."経度"
  , quadkey.quadkey
FROM gis
JOIN quadkey
  on gis."緯度"  = quadkey.lat
  AND gis."経度" = quadkey.lon
WHERE
  quadkey.quadkey LIKE '1330021123012311%'

$ sqlite3 -header a.db <a.sql
都道府県名|市区町村名|大字・丁目名|小字・通称名|街区符号・地番|緯度|経度|quadkey
東京都|新宿区|西新宿一丁目||8|35.690594|139.696837|13300211230123112000
東京都|新宿区|西新宿一丁目||7|35.691535|139.697283|13300211230123110203
東京都|新宿区|西新宿一丁目||13|35.689844|139.696710|13300211230123112022
東京都|新宿区|西新宿一丁目||13|35.689493|139.696739|13300211230123112200
東京都|新宿区|西新宿一丁目||14|35.688810|139.696712|13300211230123112222
東京都|新宿区|西新宿一丁目||6|35.692378|139.697332|13300211230123110021
東京都|新宿区|西新宿一丁目||14|35.689054|139.696938|13300211230123112220
東京都|新宿区|西新宿一丁目||12|35.689722|139.697257|13300211230123112023
東京都|新宿区|西新宿一丁目||15|35.688997|139.697313|13300211230123112221
東京都|新宿区|西新宿一丁目||9|35.690639|139.697777|13300211230123112011
東京都|新宿区|西新宿一丁目||5|35.692557|139.698498|13300211230123110103
東京都|新宿区|西新宿一丁目||11|35.689749|139.697791|13300211230123112033
東京都|新宿区|西新宿一丁目||16|35.689262|139.697834|13300211230123112213
東京都|新宿区|西新宿一丁目||16|35.688964|139.697835|13300211230123112231
東京都|新宿区|西新宿一丁目||18|35.688570|139.698069|13300211230123112322
東京都|新宿区|西新宿一丁目||10|35.689781|139.698410|13300211230123112123
東京都|新宿区|西新宿一丁目||17|35.689240|139.698448|13300211230123112303
東京都|新宿区|西新宿一丁目||18|35.688804|139.698626|13300211230123112323
東京都|新宿区|西新宿一丁目||4|35.692912|139.698873|13300211230123110110
東京都|新宿区|西新宿一丁目||1|35.688581|139.699425|13300211230123113222
東京都|新宿区|西新宿一丁目||2|35.692873|139.699397|13300211230123110111
東京都|新宿区|西新宿一丁目||2|35.692862|139.699638|13300211230123111000
東京都|新宿区|新宿三丁目||23|35.692761|139.700709|13300211230123111011
東京都|新宿区|新宿三丁目||24|35.692649|139.701359|13300211230123111103
東京都|新宿区|新宿三丁目||37|35.690187|139.701458|13300211230123113121
東京都|新宿区|新宿三丁目||26|35.691574|139.701713|13300211230123111312
東京都|新宿区|新宿三丁目||36|35.690308|139.701923|13300211230123113113
東京都|新宿区|新宿三丁目||27|35.691036|139.701941|13300211230123111333
東京都|新宿区|新宿三丁目||26|35.691737|139.702148|13300211230123111311
東京都|新宿区|新宿三丁目||25|35.692480|139.701982|13300211230123111113
東京都|新宿区|西新宿一丁目||1|35.690402|139.699375|13300211230123112113
東京都|新宿区|新宿三丁目||37|35.689556|139.701733|13300211230123113310
東京都|新宿区|新宿三丁目||38|35.690743|139.700583|13300211230123113011

こんな感じになります。

と言ってもいまいちよくわからないので、

可視化してみます

c.sql

SELECT
  json_group_array(json_array(lon, lat))
FROM quadkey
WHERE
  quadkey.quadkey LIKE '1330021123012311%'

a.jq

{
  "type": "MultiPoint",
  "crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
    }
  },
  "coordinates": .
}

# !/bin/bash

printf "https://www.townsendjennings.com/geo/?geojson=%s" $(sqlite3 a.db <c.sql | jq -c -f a.jq)

とかやって URL を生成。
これをブラウザで開くと、

緯度経度で Point がマークされて可視化してくれます。
先程の範囲内に収まっているのが確認できました。

本編ここまで。

隣接領域の話

さて、quadkey は領域なので、緯度経度と検索半径によってははみ出てしまうため、
周辺8タイルを含めた合計9タイル（将棋の王将の動く範囲と同じ）を探索領域にすると思います。

引用元

上図で例えると、このあたりとかですね。

この、隣接領域の quadkey を算出するロジックについての説明です。
＃わかりやすく説明してるサイトが見付からなかったので自分で書くことにした

こちらをベースに説明します。
https://github.com/arc279/calculate-quadkey/blob/master/src/python/quadkey_test.py

結論

先に結論を述べると、quadkey(4進数) を 2進数表記した時に、

LSB 側から1始まりで数えて、
奇数ビット目が横のタイル座標
偶数ビット目が縦のタイル座標

を表しています。
2進数2桁 で 4進数1桁 を表してるのが要点です。

と言っても、ビット演算に馴染みがないと何を言っているのかわからねーと思うので、以下で説明します。

横の隣接

level 4 横1周 で説明します。

'0200'
'0201'
'0210'
'0211'
'0300'
'0301'
'0310'
'0311'
'1200'
'1201'
'1210'
'1211'
'1300'
'1301'
'1310'
'1311'

は、横に1タイルずつずらしていって、端から端まで1周するまでの quadkey で、4進数です。
level 4 なので、 2^4 = 16 タイルで1周します。

ここで、4進数 → 2進数に直してみます。

quadkey binstr
===============
0200 0b00100000
0201 0b00100001
0210 0b00100100
0211 0b00100101
0300 0b00110000
0301 0b00110001
0310 0b00110100
0311 0b00110101
1200 0b01100000
1201 0b01100001
1210 0b01100100
1211 0b01100101
1300 0b01110000
1301 0b01110001
1310 0b01110100
1311 0b01110101

2進数表記をよく見ると、LSB側から数えて第1bit、第3bit、第5bit... と1つ飛ばしで繰り上がっているのがわかります。

なので、偶数ビットを除去してLSB側に詰めた値に、横方向のオフセットを加算すれば、正しく繰り上がりが処理できます。
ループするのではみ出さないようにビットマスクも必要です。
https://github.com/arc279/calculate-quadkey/blob/master/src/python/quadkey.py#L38

縦の隣接

level 4 縦1周 で説明しますが、横の隣接と要領は同じです.
こちらはLSB側から数えて第2bit、第4bit、第6bit... と偶数ビットに注目します。

quadkey binstr
===============
1010 0b01000100
1012 0b01000110
1030 0b01001100
1032 0b01001110
1210 0b01100100
1212 0b01100110
1230 0b01101100
1232 0b01101110
3010 0b11000100
3012 0b11000110
3030 0b11001100
3032 0b11001110
3210 0b11100100
3212 0b11100110
3230 0b11101100
3232 0b11101110

ビットマスクが必要なのも同様です。
https://github.com/arc279/calculate-quadkey/blob/master/src/python/quadkey.py#L39

タイル座標移動した後は

奇数ビットと偶数ビットで分割したのと逆の作業で統合します。
https://github.com/arc279/calculate-quadkey/blob/master/src/python/quadkey.py#L44

何故か python に itoa が無いので、基数変換で文字列に直すのに numpy.base_repr 使ってます。
c++ とかで書くなら itoa で radix 指定すれば4進数に変換できます。

おわり。
効率的なビット演算のやり方とかあったら誰か教えてください。

追記

と思ったら、Z階数曲線の wikipedia にビット演算のやり方載ってた…
https://ja.wikipedia.org/wiki/Z%E9%9A%8E%E6%95%B0%E6%9B%B2%E7%B7%9A#%E3%83%87%E3%83%BC%E3%82%BF%E6%A7%8B%E9%80%A0%E3%81%A8%E5%BA%A7%E6%A8%99%E5%80%A4%E7%BE%A4%E3%81%AE%E5%8F%96%E3%82%8A%E6%89%B1%E3%81%84