概要
- 本記事では、逆ジオコーディングの高速化と精度向上のための手法を解説します。
- 既存実装はGCSへの都度リクエストにより遅延（0.2～0.3秒）が発生していましたが、新実装ではデータ構造を圧縮し、メモリ上で処理可能にしました。
- geocellを用いたTrie構造を採用し、アルゴリズム最適化（末尾削減、エッジ補正、ノード統合、遷移圧縮）を行い、データサイズを2.44GiB→56.9MiBに縮小。
- これにより、行政区域の特定が高速化し、住所情報の解決がより効率化されました。
目次
イントロダクション
逆ジオコーディングは、ここでは特定の座標（緯度・経度）からその座標が属する行政区域（都道府県・市区町村）を特定するプロセスとします。また、緯度・経度としては北緯が正、南緯が負、東経が正、西経が負とします。たとえば、緯度=35, 経度=135とすると、"兵庫県西脇市"を返すという具合です。 SNSピリカでは、ごみを拾った投稿の座標から、逆ジオコーディングを行うことで所属する行政区画を特定し、その行政区画でごみが拾われたという記録として、統計や、見える化ページへの表示のために使っています。
このほど、既存の実装をフルスクラッチして新規に実装し、高速化を行いました。
より詳細な要件は以下になります。
また行政区域として、国土数値情報 | 行政区域データを採用します。
ここには所属未定地というものもあり、2024年現在では、たとえば新海面処分場(例: 35.587713, 139.819156)は所属未定地になります。これについてはNoneを返すこととします。
余談ですが、この行政区域データ、
- 124,325個以上、合計15,395,089頂点の多角形
- 1,904市区町村
となっていて、競技プログラマーなら心躍るサイズになっています。
用語の説明
geocell
geocellはgeomodelで提唱されている概念です。
緯度[-90,90], 経度[-180,180]の長方形領域を起点とし、等間隔に4x4で16分割を行って、指定座標を含むように長方形を再帰的に選びます。分割それぞれに対して16進数が割り当てられ、広い方から狭い方へ数が並べられます。たとえば、東京駅(35.6811124, 139.764516)を含む9桁のgeocellはde0d9eaf1となります。これは 緯度[35.68084716796875, 35.68153381347656] × 経度[139.76394653320312 139.76531982421875]の領域を表現します。
アイデアとしてはジオハッシュと同様です。
以下はgeomodelのソースコードにある16進数割り当ての図です。
+---+---+---+---+ (90, 180)
| a | b | e | f |
+---+---+---+---+
| 8 | 9 | c | d |
+---+---+---+---+
| 2 | 3 | 6 | 7 |
+---+---+---+---+
| 0 | 1 | 4 | 5 |
(-90,-180) +---+---+---+---+
行政区域文字列
行政区域を表現する文字列です。表現方法はいろいろありますが、SNSピリカでは都道府県名・郡名・市区町村名・政令指定都市の行政区名が含まれていれば良いとしています。
既存の実装と問題点
まず緯度・経度を含む9桁のgeocellを計算します。これを使って、GCSにある静的ファイルから行政区域文字列を得ます。
たとえば東京駅の緯度・経度であれば、パスde/de0dにあるファイルの16進数で9eaf1番目のレコードに行政区域文字列が格納されています。このレコードの部分だけ範囲リクエストして取得します。
各geocellに行政区域文字列をどのように割り当てたのでしょう？オリジナルのコードがないので推測ですが、geocellの中央の座標が属する行政区域を割り当てたと考えられます。行政区域自体は重複はしないことを利用すると、次のアルゴリズムで効率的に計算できます。
- 経度をx座標, 緯度をy座標として平面上で考える。
- geocellの中心点を並べると、同じy座標に多数の点が並んでいる。y座標を固定して
とする。
- それぞれの行政区域データについて、辺が直線
を通るなら、その交点のx座標をリスト化する。
- 辺が直線
の要求する最小単位は9桁の場合
であり、少なくとも
が整数にならないといけません。このとき
なので、座標の最小単位
では足りません。
- 辺が直線
- リストをソートし、0番目と1番目の間、2番目と3番目の間、…が行政区域データに含まれる。
この方法は金銭的コストもそこまでかからないですし良く見えますが、毎回GCSへのリクエストが走ってしまう点が問題です。1リクエストあたり0.2~0.3秒程度が追加でかかってしまいます。
新しい実装の詳細
この行政区域文字列を引くためのデータ構造はAPIを処理するインスタンスのメモリに乗るのではないかと考えました。 既存の実装において、中心点がいずれかの行政区域に属しているgeocellの個数は41,018,595個なので、そのまま1個ずつに行政区域文字列を置くと2.44GiB程度かかります。 ですが、行政区域は1,904種類しかないことと、多角形データはそこまで特異なデータではないため、もっと簡易的に表現できます。
geocellの文字列をtrie)で表すことを考えます。
16進数を文字列として扱ったので、遷移はたかだか16個になります。
また、trieの葉は行政区域文字列を持つノードにつなぐこととします。
実装は、build_trie関数の、for k, vの節が該当します。
東京駅のgeocellについて注目したtrieは以下のようになります。これは最終的な改善後に得られたtrieであることに注意してください。
できあがったtrieは2,785,790頂点になります。この時点で1569MiB程度になります。メモリの計測はmemory-profilerを使用しました。
trieを用いて逆ジオコーディングを行う時は、既存の実装と同様に、緯度・経度を含む9桁のgeocellを計算した後、trieのルートノードから1文字ずつ辿っていって到達した行政区域文字列が戻り値になります。
さらなる改良
自明な末尾を削る
当たり前のことですが、同じ行政区域は平面的には近い位置関係にあります。したがって、先頭8文字が同じgeocellのすべてが同じ行政区画に属することは非常によくあります。これを検出して、末尾1文字を削った文字列をたどった時点で当該行政区域のノードに行かせるようにします。これを、より小さい文字数でも行います。
これの効果は絶大で、頂点数が293,291になります。
実装は、shrink関数が該当します。
中心点が無所属かつ唯一の行政区域の辺がgeocellに含まれる場合、その行政区域をgeocellに割り当てる
geocellのなかには、陸地はあるが無所属と判定されるものもあります。 たとえば、この地図において、地図の中心がgeocellの中心、そこを中心として左右約107m, 上下約76mが9桁のgeocellに相当します。geocellの中心だけ見ると海であり、以前はこのgeocellは無所属と判定されていました。したがって宗谷岬付近でごみ拾いをしても、無所属として判定され、稚内市でのごみ拾いの記録には含まれなくなります。 このように、海岸でごみ拾いを投稿した場合、行政区域に所属しなくなることがあります。これをカバーします。
各多角形の辺について、一方の端点から他方の端点まで隣接geocellをたどっていくことを考えます。行政区域データには極端なデータはないため、現実的な時間で処理することができます。たどったgeocellにマーキングをして、最初の実装でgeocellが無所属で、かつマーキングが唯一の行政区域によってなされている場合、その行政区域をgeocellに割り当てます。
以下は処理のイメージ図です。赤い線分が多角形の辺、長方形がgeocell, 色のついた長方形が、マーキングされたgeocellで、番号はマークした順番を表します。
これにより299,713頂点となり、2%程度データは増えますが、無所属と判定される陸地を減らせるので良い改善だと思います。
実装は、plot_edge関数とその返り値を用いた処理に該当します。
無所属を許して自明な末尾を削る (不採用)
自明な末尾を削るで、単一の行政区域の他に、無所属も許して縮約します。
これにより204,907頂点になりますが、陸上でごみ拾いをしたのに沖合でごみ拾いをしたような投稿が、自治体の見える化ページに表示できてしまうということで却下になりました。
実装は、lift_unique_edges関数が該当します。
同一内容のノードをまとめる
ほとんどのノードは2行政区域の境界を表します。そうなると、全く同じ内容のノードが高頻度で現れます。これをまとめます。trieが木じゃなくなってきますがしょうがない。
229,902頂点になります。実装は、identify_same_rows関数が該当します。
遷移を表すリストを簡約化する
上記と同様の理由で、同じノードからの遷移先が同じになるケースが頻発します。これを次のようにまとめます。
[(リンクパターン), リンク先0, リンク先1, ... ]
リンク先の個数をtとしたとき、tを表現可能な最小ビット数をuとします。t→uは、たとえば2→1, 3→2, 4→2, 5→3, 8→3, 9→4, 16→4です。
リンクパターンは1桁の16通りのリンク先を、それぞれuビットで表し、それを間隔uビットで16個重ねた整数です。
たとえば、簡約前のノードの内容が [-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,1899,-1,8377]の場合、簡約後は[2214592512,-1,1899,8377]となります。この場合、t=3, u=2であり、2ビットごとなのでリンクパターンは22=4進数とみても差し支えないです。先頭要素を4進数で表すと
2214592512 = 2010000000000000_4となります。つまり、
- 桁がd(13)のときは、リンク先1 (1899)に飛ぶ
- 桁がf(15)のときは、リンク先2 (8377)に飛ぶ
- それ以外はリンク先0 (-1)に飛ぶ。つまり無所属とする と解釈できます。
これにより、頂点数は変化しませんが、多くのノードが16要素から3要素になる等して4割弱程度のメモリ圧縮が期待できます。
実装は、shrink_by_pattern関数が該当します。
新しい実装の結果
このデータ構造の圧縮により、逆ジオコーディングのためのデータ構造が、SNSの実行インスタンスのメモリに乗せることができるようになります。計測したメモリは56.9MiB程度になりました。 その結果、GCSのリクエストのオーバーヘッド分がなくなり、リクエストの高速化ができました。コールドスタートを避けて10回測った結果では、以下のようになり、0.2~0.3秒の高速化になりました。
最終的なコード
pirika-inc/revgeocode_trie_builder: 日本国内向け逆ジオコーディング用インメモリデータ構造構築ツール
まとめ・雑感
逆ジオコーディングの処理の高速化について紹介しました。最適化している間は楽しくて夢中になっていました。弊社の業務に限らず汎用性の高いタスクだと感じました。
