効率的な多次元空間ポイントインデックスアルゴリズム - Geohash と Google S2

[[201793]]

導入

毎晩残業して家に帰るときは、Didiやシェア自転車を使うこともあります。アプリを開くと、次のインターフェースが表示されます。

アプリのインターフェースには、近くの特定の範囲内で利用可能なタクシーまたはシェア自転車が表示されます。地図には、自分を中心に 5 キロメートルの範囲内の車が表示されるとします。どのように実現するのでしょうか? 最も直感的なアイデアは、データベース内のテーブルを検索し、ユーザーからの距離が 5 キロメートル以下の車を計算してクエリし、それらをフィルタリングして、データをクライアントに返すことです。

このアプローチはかなり愚かであり、一般的には行われません。なぜでしょうか? このアプローチでは、テーブル全体のすべての項目の相対距離を計算する必要があるためです。時間がかかりすぎます。データ量が多すぎるため、分割して処理する必要があります。次に、マップをブロックに分割することを考えます。この方法では、各ブロック内の各データに対して相対距離を 1 回計算する場合でも、テーブル全体を 1 回計算するよりもはるかに高速になります。

現在最も一般的に使用されているデータベースである MySQL と PostgreSQL はどちらもネイティブで B+ ツリーをサポートしていることは誰もが知っています。このデータ構造は効率的にクエリできます。マップのセグメンテーションのプロセスは、実際にはインデックスを追加するプロセスです。マップ上のポイントに適切なインデックスを追加して並べ替える方法を考えれば、バイナリ検索に似た方法を使用して高速クエリを実行できます。

問題は、地図上の点が経度と緯度の 2 次元であることです。どのようにインデックスを作成するのでしょうか。経度または緯度のいずれかの次元のみを検索すると、一度検索した後で再度検索する必要があります。もっと高次元はどうでしょうか? 3次元です。ジョイントキーの連結など、次元の優先順位を設定できるという人もいるかもしれません。では、3次元空間では、x、y、zのどの次元の優先順位が高いのでしょうか。優先順位を設定することは、あまり合理的ではないようです。

この記事では、さらに 2 つの一般的な空間ポイント インデックス アルゴリズムを紹介します。

1. ジオハッシュアルゴリズム

1. Genhashアルゴリズムの紹介

Genhash は Gustavo Niemeyer によって発明されたジオコーディング技術です。空間をグリッドに分割する階層型データ構造です。 Genhash は、空間充填曲線における Z 次曲線の実用的応用です。

Z オーダー曲線とは何ですか?

上の図はZオーダー曲線です。この曲線は比較的シンプルで、生成も簡単です。各 Z の端と端を接続するだけです。

Z オーダー曲線は 3 次元空間に拡張することもできます。 Z 形状が十分に小さく、密度が高ければ、3 次元空間全体を埋めることができます。

この時点で、読者はまだ混乱していて、ジオハッシュと Z 曲線の関係がわからないかもしれません。実際、ジオハッシュ アルゴリズムの理論的根拠は、Z 曲線の生成原理に基づいています。 Geohashに戻りましょう。

Geohash は任意の精度のセグメンテーション レベルを提供できます。一般的な評価は 1 から 12 までです。

引用文で言及されている問題を覚えていますか? ここでは、Geohash を使用してこの問題を解決できます。

Geohash 文字列の長さを使用して、分割する領域のサイズを決定できます。この対応は、上の表のセルの幅と高さを参照できます。セルの幅と高さを選択すると、Geohash 文字列の長さが決定されます。このようにして、地図を長方形の領域に分割します。

地図上でエリアは分割されていますが、ポイントとポイントの近くのエリアをすばやく見つける方法という、まだ解決されていない問題があります。

Geohash には、Z オーダー曲線に関連するプロパティがあります。つまり、あるポイントに近いハッシュ文字列 (絶対的ではありませんが) には常に共通のプレフィックスがあり、共通のプレフィックスの長さが長いほど、2 つのポイントは近くなります。

この特性により、Geohash は一意の識別子としてよく使用されます。 Geohash はデータベース内のポイントを表すために使用できます。 Geohash の共通プレフィックス機能を使用すると、近隣のポイントをすばやく検索できます。ポイントが近いほど、ターゲット ポイントとのジオハッシュ文字列の共通プレフィックスが長くなります (ただし、必ずしもそうとは限りません。特別なケースもあります。次の例で説明します)

Geohash にもいくつかのエンコード形式があり、そのうちの 2 つは base 32 と base 36 です。

ベース 36 バージョンでは大文字と小文字が区別され、36 文字 (「23456789bBCdDFgGhHjJKlLMnNPqQrRtTVWX」) が使用されます。

2. Geohashアプリケーションの例

次の例では、base-32 を例として使用します。例えば。

上の写真は、中央にメトロ シティがある地図です。メトロ シティに最も近いレストランを探す必要がある場合、どうすれば見つけられるでしょうか?

最初のステップは、ジオハッシュを使用してマップをグリッド化することです。表を参照して、文字列の長さが 6 の長方形を選択し、マップをグリッド化します。

クエリの結果、Metro Cityの緯度と経度は[31.1932993, 121.43960190000007]です。

まず緯度について考えてみましょう。地球の緯度範囲は[-90,90]です。この区間を[-90,0)と[0,90]の2つの部分に分割します。 31.1932993 は区間 (0,90)、つまり右の区間に位置し、1 とマークされています。次に、区間 (0,90) を [0,45]、[45,90] と 2 つに分割し続けます。31.1932993 は区間 [0,45)、つまり右の左の区間に位置し、0 とマークされています。分割を続けます。

左区間の中央値 右区間のバイナリ結果

次に、同じ方法で経度を処理します。地球の経度範囲は[-180,180]です

左区間の中央値 右区間のバイナリ結果

緯度によって生成されるバイナリ文字列は 101011000101110 で、経度によって生成されるバイナリ文字列は 110101100101101 です。「経度は偶数位置、緯度は奇数位置」という規則に従って、経度と緯度のバイナリ文字列が再結合され、新しい文字列 111001100111100000110011110110 が生成されます。最後のステップは、この最終的な文字列を文字に変換することです。これには、base-32 テーブルを参照する必要があります。 11100 11001 11100 00011 00111 10110 を 10 進数に変換すると 28 25 28 3 7 22 になります。表を参照して最終結果 wtw37q を取得します。

周囲の 8 つのグリッドを個別に計算することもできます。

地図からわかるように、これら 9 つの隣接するグリッドのプレフィックスはまったく同じです。それらはすべてwtw37です。

文字列にもう 1 桁追加するとどうなるでしょうか? Geohash は 7 桁に増えます。

ジオハッシュが 7 桁に増えるとグリッドが小さくなり、Metro City のジオハッシュは wtw37qt になります。

これを読んだ後、読者は Geohash のアルゴリズム原理を明確に理解できるはずです。 6桁と7桁を1つの絵に組み合わせてみましょう。

中央の大きなグリッドのジオハッシュ値は wtw37q なので、その中のすべての小さなグリッドのプレフィックスも wtw37q であることがわかります。 Geohash 文字列の長さが 5 の場合、Geohash は wtw37 になるはずです。

次に、先ほど述べたジオハッシュとZオーダー曲線の関係について説明します。最後のステップで経度と緯度の文字列を結合するルールを思い出してください。「偶数桁は経度、奇数桁は緯度です。」読者は、このルールがどこから来たのか、少し興味があるはずです。ただ空想で作られたのでしょうか? 実はそうではありません。このルールは、Z オーダー曲線です。下の写真をご覧ください:

x 軸は緯度、y 軸は経度です。このように、経度は偶数で配置され、緯度は奇数で配置されます。

最後に、正確性の問題があります。以下の表のデータの一部は Wikipedia から引用したものです。

3. ジオハッシュの実装

この時点で、読者は Geohash アルゴリズムを明確に理解しているはずです。次に、Go を使用して Geohash アルゴリズムを実装します。

パッケージ geohash 
 
輸入 （
 「バイト」  
 ） 
 
定数(
    BASE32 = "0123456789bcdefghjkmnpqrstuvwxyz"  
    MAX_LATITUDE float64 = 90
    MIN_LATITUDE フロート64 = -90
    MAX_LONGITUDE float64 = 180
    MIN_LONGITUDE float64 = -180
 ） 
 
 var (
    ビット = [] int {16, 8, 4, 2, 1}
    base32 = []バイト(BASE32)
 ） 
 
ボックス構造体型{
    MinLat、MaxLat float64 // 緯度
    MinLng、MaxLng float64 // 経度
} 
 
 func (this *Box) 幅() float64 {
 this.MaxLng - this.MinLngを返す
} 
 
 func (this *Box) Height() float64 {
 this.MaxLat - this.MinLatを返す
} 
 
 // 入力値: 緯度、経度、精度 (ジオハッシュの長さ)
 // ポイントが位置するジオハッシュとエリアを返します
func Encode(緯度, 経度 float64,精度  int ) (文字列、*ボックス) {
    var geohash バイト.バッファ
    var minLat、maxLat float64 = MIN_LATITUDE、MAX_LATITUDE
    var minLng、maxLng float64 = MIN_LONGITUDE、MAX_LONGITUDE
    var mid float64 = 0 
 
ビット、ch、長さ、isEven := 0、0、0、 true  
長さ<精度{
        偶数の場合{
            中間値 = (最小長 + 最大長) / 2 の場合、中間値 < 経度 {
                ch |= ビット[ビット]
                最小長さ = 中
            }それ以外{
                最大長 = 中
            }
        }それ以外{
            中間 = (最小緯度 + 最大緯度) / 2 の場合、中間 < 緯度 {
                ch |= ビット[ビット]
                最小緯度 = 中央
            }それ以外{
                最大緯度 = 中央
            }
        } 
 
        偶数 = !偶数
        ビット<4の場合{
ビット++
        }それ以外{
            ジオハッシュ.WriteByte(base32[ch])
            長さ、ビット、ch = 長さ+1、0、0
        }
    } 
 
    b := &ボックス{
        最小緯度: 最小緯度、
        最大緯度: 最大緯度、
        最小長さ: 最小長さ、
        最大長: 最大長、
    } 
 
 geohash.String()を返す、b
 }

4. ジオハッシュの長所と短所

Geohash の利点は明らかで、エンコードに Z オーダー曲線を使用します。 Z オーダー曲線は、2 次元または多次元空間内のすべての点を 1 次元曲線に変換できます。数学的にはフラクタル次元と呼ばれます。また、Z オーダー曲線には局所的な順序保存もあります。

Z オーダー曲線は、点の座標値のバイナリ表現をインターリーブすることによって、複数の次元で点の Z 値を単純に計算します。データがソートに追加されると、バイナリ検索ツリー、B ツリー、スキップ リスト、または (下位ビットが切り捨てられた) ハッシュ テーブルなどの任意の 1 次元データ構造を使用してデータを処理できます。 Z オーダー曲線によって得られる順序は、四分木の深さ優先トラバーサルから得られる順序と同等に記述できます。

これは Geohash のもう 1 つの利点であり、これにより近隣のポイントの検索が高速化されます。

Geohash の欠点の 1 つは、Z オーダー曲線からも生じます。

Z オーダー曲線にはさらに深刻な問題があります。局所的には順序が維持されますが、突然変異の影響を受けやすいのです。 Z の文字が変わるたびに、順序が変化する可能性があります。

上の写真にマークされている青い点を見てください。どの 2 つの点も隣接していますが、離れています。右下隅の図を見ると、隣接する値を持つ 2 つの赤い点の間の距離は、正方形全体の辺の長さとほぼ同じです。緑の点に隣接する 2 つの値も、正方形の長さの半分に達します。

Geohash のもう 1 つの欠点は、グリッド サイズが適切に選択されていない場合、隣接するポイントを特定するのが困難になる可能性があることです。

上の図を見ると、ジオハッシュ文字列が 6 として選択されている場合、大きな青いグリッドが表示されます。赤い星はメロ市、紫色の点は探索中に見つかった目標地点です。 Geohash アルゴリズムを使用してクエリを実行すると、近いものは wtw37p、wtw37r、wtw37w、wtw37m になる可能性があります。しかし、実際には最も近いポイントは wtw37q です。このような大きなグリッドを選択した場合は、周囲の 8 つのグリッドを検索する必要があります。

Geohash 文字列が 7 の場合、小さな黄色のボックスになります。このように、赤い星に最も近い点は 1 つだけです。 wtw37qwです。

グリッドのサイズと精度が適切に選択されていない場合、最も近いポイントを照会するには、周囲の 8 つのポイントを再度照会する必要があります。

2. 空間充填曲線とフラクタル

2 番目の多次元空間ポイント インデックス アルゴリズムを紹介する前に、空間充填曲線とフラクタルについて説明する必要があります。

多次元空間の点のインデックス付けの問題を解決するには、2 つの問題を解決する必要があります。まず、多次元空間を低次元または 1 次元空間に縮小する方法は? 次に、1 次元曲線をフラクタル化する方法はどうでしょうか?

1. 空間充填曲線

数学的分析には、「無限に長い線は、任意の次元空間のすべての点を通過することができるか？」という難しい問題があります。

1890 年、ジュゼッペ ペアノは、単位正方形上のすべての点を通過する連続曲線 (現在ではペアノ曲線と呼ばれています) を発見しました。彼の目標は、単位間隔から単位正方形への連続的なマッピングを構築することでした。ペアノは、単位区間内の無限数の点の濃度は、任意の有限次元多様体内の無限数の点の濃度と同じであるという、ゲオルク・カントールの以前の直観に反する結果に触発されました。ペアノが解決する問題の本質は、そのような連続的なマッピング、つまり平面全体を埋めることができる曲線が存在するかどうかです。上の写真は彼が見つけた曲線です。

一般的に言えば、1 次元のオブジェクトが 2 次元のグリッドを埋めることは不可能です。しかし、ペアノ曲線は反例を与えます。ペアノ曲線は連続だがどこでも微分不可能な曲線です。

ペアノ曲線の構築方法は次のとおりです。正方形を 9 つの等しい小さな正方形に分割し、左下隅の正方形から右上隅の正方形までの線分を使用して、小さな正方形の中心を接続します。次のステップでは、各小さな正方形を 9 つの等しい正方形に分割し、上記の方法でそれらの中心を接続します。この操作手順を無限に繰り返し、最終的に得られる極限ケースの曲線をペアノ曲線と呼びます。

ペアノは、[0, 1]区間上の点と正方形上の点の間のマッピングについて詳細な数学的記述を与えた。実際、これらの正方形の点は

、 x と y が単位正方形に属するすべての値を取るような 2 つの連続関数 x = f(t) と y = g(t) を見つけることができます。

1年後の1891年、ヒルベルトはヒルベルト曲線と呼ばれるこの曲線を作成しました。

上の図は1-6次のヒルベルト曲線です。具体的な構築方法については次の章で説明します。

上の図は、3D 空間を埋め尽くすヒルベルト曲線を示しています。

空間充填曲線には、ドラゴン曲線、ゴスパー曲線、コッホ曲線、ムーア曲線、シェルピンスキー曲線、オスグッド曲線など、さまざまなバリエーションがあります。これらの曲線はこの記事とは関係がないので、詳しくは紹介しません。

数学的解析において、空間充填曲線は、単位間隔を単位正方形、立方体、またはより一般的には n 次元超立方体内の連続曲線にマッピングし、パラメータが増加するにつれて単位立方体内の特定の点に任意に近づくことができるパラメータ化された注入関数です。空間充填曲線は、数学的な重要性に加えて、次元削減、数理計画法、スパース多次元データベースのインデックス作成、電子工学、生物学でも使用されます。空間充填曲線は現在、インターネット マップで使用されています。

2. フラクタル

ペアノ曲線の出現は、人々の次元に対する理解に誤りがあることを示しており、次元の定義を再検討する必要があることを示しています。これがフラクタル幾何学が考えることです。フラクタル幾何学では、次元は分数になる場合があり、これをフラクタル次元と呼びます。

多次元空間で次元削減した後にフラクタル化する方法も問題です。フラクタルを作成する方法はたくさんあります。ここでは、フラクタルを作成する方法と、各フラクタルのフラクタル次元、つまりハウスドルフ フラクタル次元とトポロジカル次元のリストを示します。ここではフラクタルについて詳しく説明しません。興味のある方はリンク先の内容を注意深く読んでみてください。

次に、多次元空間ポイントインデックスアルゴリズムについて説明します。次のアルゴリズムの理論的根拠はヒルベルト曲線にあります。まずヒルベルト曲線について詳しく説明します。

3. ヒルベルト曲線

1. ヒルベルト曲線の定義

ヒルベルト曲線は、1891 年にデイヴィッド・ヒルベルトによって提案された、平面の正方形を埋めるフラクタル曲線 (空間充填曲線) です。

平面を埋め尽くすので、ハウスドルフ次元は 2 です。埋める正方形の辺の長さを 1 とすると、n 番目のステップでのヒルベルト曲線の長さは 2^n - 2^(-n) です。

2. ヒルベルト曲線の構築法

一次ヒルベルト曲線は、正方形を 4 つの等しい部分に分割し、その 1 つのサブ正方形の中心から始めて、残りの 3 つの正方形の中心を通る線を引くことによって生成されます。

2 次ヒルベルト曲線は、各サブ正方形を 4 つの等しい部分に分割し、4 つの小さな正方形ごとに 1 次ヒルベルト曲線を生成することによって生成されます。次に、4 つの 1 次ヒルベルト曲線を端から端まで接続します。

3 次ヒルベルト曲線の生成方法は 2 次ヒルベルト曲線の生成方法と似ています。まず、2 次ヒルベルト曲線を生成します。次に、4 つの 2 次ヒルベルト曲線を端から端まで接続します。

n 次ヒルベルト曲線を生成する方法も再帰的です。まず、n-1 次ヒルベルト曲線を生成し、次に 4 つの n-1 次ヒルベルト曲線を端から端まで接続します。

3. ヒルベルト曲線を選択する理由は何ですか?

これを読んで、次のような疑問を持つ読者もいるかもしれません。空間充填曲線は数多くありますが、なぜヒルベルト曲線が選ばれたのでしょうか。

ヒルベルト曲線は非常に優れた特性を持っているからです。

（１）次元削減

まず、空間充填曲線として、ヒルベルト曲線は多次元空間の次元を効果的に削減できます。

上の図は、平面を塗りつぶし、平面上のすべての点を 1 次元の線に拡張した後のヒルベルト曲線を示しています。

上の図のヒルベルト曲線は 16 点しか通らないのに、どうやって平面を表現できるのかと疑問に思う人もいるかもしれません。

もちろん、n が無限大に近づくと、n 次ヒルベルト曲線は平面全体をほぼ埋めることができます。

（２）安定性

n 次のヒルベルト曲線において、n が無限大に近づくと、曲線上の点の位置は基本的に安定する傾向があります。例えば：

上の図の左側がヒルベルト曲線、右側が蛇行曲線です。 n が無限大に近づくと、理論的には両方とも平面を埋めることができます。しかし、なぜヒルベルト曲線の方が優れているのでしょうか?

蛇行曲線上の点を考えると、n が無限大に近づくにつれて、蛇行曲線上のこの点の位置は常に変化します。

これにより、ポイントの相対的な位置は常に不確実になります。

ヒルベルト曲線をもう一度見てみましょう。これも n が無限大に近づく点です。

上の図からわかるように、点の位置はほとんど変わっていません。したがって、ヒルベルト曲線の方が優れています。

（３）連続

ヒルベルト曲線は連続しているので、空間を埋めることが保証されます。連続性には数学的な証明が必要です。具体的な証明方法についてはここでは詳しく述べません。興味のある方は、記事の最後にあるヒルベルト曲線に関する論文をクリックしてください。そこに連続性の証明が含まれています。

次に紹介する Google S2 アルゴリズムはヒルベルト曲線に基づいています。これで読者はヒルベルト曲線がなぜ選ばれたのか理解できるはずです。

4. S²アルゴリズム

• GoogleのS2ライブラリは、空間インデックス機能だけでなく、4年以上前にリリースされたライブラリであり、十分な注目を集めなかったという点でも、真の宝物です。

上記の文章は、2015 年に Google のエンジニアが投稿したブログ記事からの抜粋です。彼は、S2 アルゴリズムがリリースから 4 年経っても、十分な評価を受けていないことを心から嘆きました。しかし、現在ではS2は大手企業でも利用されています。

このヘビーウェイトアルゴリズムを紹介する前に、まずこのアルゴリズムの名前の由来を説明したいと思います。 S2 は、実際には単位球を表す幾何数学の数学記号 S² に由来しています。 S2 ライブラリは、実際には球面上のさまざまな幾何学的問題を解決するために設計されています。 40% しか完成していない golang 公式リポジトリの geo/s2 を除き、Java、C++、Python を含む他の言語での S2 実装はすべて 100% 完成していることは注目に値します。この記事では、このバージョンの Go に焦点を当てます。

次に、S2 を使用して多次元空間の点のインデックス付けの問題を解決する方法を見てみましょう。

1. 球面座標変換

多次元空間を扱うためのこれまでの考え方によれば、まず次元を削減する方法を検討し、次にフラクタル化する方法を検討します。

ご存知のとおり、地球はほぼ球体です。球体は3次元の物体です。この3次元の物体を1次元の物体にするにはどうすればよいでしょうか。

球面上の点は、直交座標系では次のように表すことができます。

 x = r * sinθ * cosφ
 y = r * sinθ * sinφ
 z = r * cosθ

通常、地球上の地点を表すには経度と緯度を使います。

さらに一歩進んで、これを球面上の経度と緯度に関連付けることができます。ただし、ここで注目すべきは、緯度の角度 α と直交座標系の球面座標 θ を合計すると 90° になるということです。したがって、三角関数の変換に注意してください。

したがって、地球上の任意の経度と緯度の点は f(x, y, z) に変換できます。

S2では、地球の半径を単位1とします。したがって、半径を考慮する必要はありません。 x、y、zの値の範囲は[-1,1] x [-1,1] x [-1,1]の間隔に制限されます。

2. 球が平面になる

次のステップS2は、球面を平面に研磨することです。それはどうやってやるのですか?

まず、下の図のように、外接立方体を地球の外側に配置します。

接線立方体の 6 つの面を球の中心から投影します。 S2 は球面上のすべての点を外接立方体の 6 つの面に投影します。

ここでは簡単な投影図を描いています。上の図の左側は、立方体の 1 つの面への投影の概略図です。実際に影響を受ける球面は右側のものです。

側面から見ると、球面の 1 つが立方体の面の 1 つに投影され、円の端と中心を結ぶ線は互いに 90° の角度をなしますが、x 軸、y 軸、z 軸に対しては 45° の角度をなします。下の図の左側に示すように、球体の 6 つの方向に 45° の補助円を描くことができます。

上の写真の左側には6本の補助線が引かれています。青い線は前後のペア、赤い線は左右のペア、緑の線は上下のペアです。これらは、45°の角度と円の中心を結ぶ線が球と交差する点の軌跡です。このようにして、上の図の右側に示すように、外接立方体の 6 つの面に投影された球面を描くことができます。

立方体に投影した後、立方体を展開することができます。

立方体を展開する方法はたくさんあります。どのように拡大しても、最小単位は正方形です。

上記はS2の射影解です。次に、他の投影ソリューションについて説明します。

まず、三角形と四角形を組み合わせた次のような方法があります。

この方法の展開図を下図に示します。

この方法は実際には非常に複雑で、サブグラフは 2 種類のグラフィックで構成されます。座標変換は少し複雑です。

もう一つの方法は、すべて三角形を使うことです。三角形の数が多いほど、球体に近づきます。

上の画像の一番左の絵は、20 個の三角形で構成されています。角が多く、球体とはかなり異なっていることがわかります。三角形の数が増えるにつれて、球体にどんどん近づいていきます。

20 個の三角形を展開すると、このように見えるようになります。

最後の方法は、現時点ではおそらく最良の方法ですが、最も複雑な方法でもあるかもしれません。六角形に合わせて投影します。

六角形は角の数が比較的少なく、6 つの辺が他の六角形と接続できます。上記の最後の図から、六角形の数が十分であれば、球体に非常に近くなることがわかります。

六角形を展開するとこのようになります。もちろん、ここには六角形が 12 個しかありません。六角形の数が多いほど良く、粒度が細かくなるほど球体に近くなります。

Uber は公開情報の中で、都市を多数の六角形に分割し、六角形のグリッドを使用していると述べました。彼らはこれを自分たちで開発すべきだった。おそらく Didi もエリアを六角形に分割しているのでしょう、あるいは Didi の方がより優れた分割プランを持っているのかもしれません。

S2 に戻ると、S2 は正方形を使用します。このようにして、最初のステップの球面座標はさらに f(x,y,z) -> g(face,u,v) に変換されます。ここで、face は正方形の 6 つの面であり、u と v は 6 つの面のうちの 1 つの面の x 座標と y 座標に対応します。

3. 球面長方形投影補正

前のステップでは、球面上の球面長方形を正方形の特定の面に投影しました。結果の形状は長方形に似ていますが、球面上の角度が異なるため、同じ面に投影された場合でも、各長方形の面積は異なります。

上の図は、球面上の正方形の面に球状の長方形を投影した状況を示しています。

実際に計算してみると、最大面積と最小面積の差は5.2倍になることがわかりました。上の写真の左側をご覧ください。同じ円弧間隔でも、異なる緯度の正方形に投影すると面積が異なります。

ここで、投影された各形状の面積を修正する必要があります。適切なマッピング補正関数をどのように選択するかが鍵となります。目標は、各長方形の面積をできるだけ同じにして、上の図の右側のような外観を実現することです。

この変換コードは C++ バージョンでは詳細に説明されていますが、Go バージョンでは簡単にしか説明されていません。長い間、私は混乱していました。

線形変換は最も高速な変換ですが、変換比は最も小さくなります。 tan() 変換により、投影された各長方形の面積の一貫性が高まり、最大長方形と最小長方形の比率はわずか 0.414 しか違わなくなります。非常に近いと言えるでしょう。しかし、tan() 関数の呼び出し時間は非常に長くなります。すべてのポイントをこのように計算すると、パフォーマンスは 3 倍低下します。

最終的に、Google は接線を近似する投影曲線である二次変換を選択しました。 tan() よりもはるかに高速に計算され、約 3 倍高速です。結果として得られる投影長方形のサイズも同様になります。ただし、最大の長方形と最小の長方形の比率は依然として 2.082 です。

上記の表で、ToPoint と FromPoint は、それぞれ単位ベクトルをセル ID に変換するのに必要なミリ秒数と、セル ID を単位ベクトルに戻すのに必要なミリ秒数です (セル ID は、立方体の 6 つの面の 1 つに投影された長方形の ID です。長方形はセルと呼ばれ、それに対応する ID はセル ID と呼ばれます)。 ToPointRaw は、何らかの目的のためにセル ID を非単位ベクトルに変換するのに必要なミリ秒数です。

S2 のデフォルトの変換は、セカンダリ変換です。

 #S2_PROJECTION と S2_QUADRATIC_PROJECTION を定義します

これら 3 つの変換がどのように実行されるかを詳しく見てみましょう。

 #S2_PROJECTION == S2_LINEAR_PROJECTIONの場合
 
インラインダブルS2::STtoUV(ダブルs) {
 2 * s - 1 を返します。
 } 
 
インラインダブルS2::UVtoST(ダブルu) {
 0.5 * (u + 1)を返します。
 } 
 
 #elif S2_PROJECTION == S2_TAN_PROJECTION 
 
インラインダブルS2::STtoUV(ダブルs) {
  // 残念ながら、tan(M_PI_4)は1.0よりわずかに小さいです。これは 
  // tan()の実装に欠陥があり、  
  // x=pi/4でのtan(x)は2であり、 2つの隣接する浮動小数点
  // π/4の無限精度値の両側の点数は
  // 1.0 よりわずかに下とわずかに上になる接線 
  // 最も近い倍精度の結果。 
 
  s = tan(M_PI_2 * s - M_PI_4);
 s + (1.0 / (GG_LONGLONG(1) << 53)) * sを返します。
 } 
 
インラインダブルS2::UVtoST(ダブルu) {
  揮発性ダブルa = atan(u);
 (2 * M_1_PI) * (a + M_PI_4)を返します。
 } 
 
 #elif S2_PROJECTION == S2_QUADRATIC_PROJECTION 
 
インラインダブルS2::STtoUV(ダブルs) {
  s >= 0.5 の場合、(1/3.) * (4*s*s - 1)を返します。
それ以外           (1/3.) * (1 - 4*(1-s)*(1-s))を返します。
 } 
 
インラインダブルS2::UVtoST(ダブルu) {
  u >= 0 の場合、0.5 * sqrt(1 + 3*u)を返します。
それ以外         1 - 0.5 * sqrt(1 - 3*u) を返します。
 } 
 
 ＃それ以外  
 
 #error S2_PROJECTIONの値が不明です
 
 #終了

上記ではtan(M_PI_4)の処理が行われていますが、精度上の理由から1.0よりわずかに小さくなっています。

したがって、投影後の補正関数の 3 つの変換は次のようになります。

 // 線形変換
u = 0.5 * ( u + 1 ) です。 
 
 // tan() 変換
u = 2 / π * (atan(u) + π / 4) = 2 * atan(u) / π + 0.5 
 
 // 二次変換
u >= 0、u = 0.5 * sqrt(1 + 3*u)
 u < 0、u = 1 - 0.5 * sqrt(1 - 3*u)

上記の変換式では u のみが書かれていますが、u のみが変換されるわけではないことに注意してください。実際の使用では、u と v は入力パラメータとして扱われ、変換されます。

この修正関数は、Go バージョンで二次変換のみを実装します。他の 2 つの変換方法については、ライブラリ全体ではまったく言及されていません。

 // stToUV は sまたはt 値を対応する uまたはv 値に変換します。
 // これは[-1,1]から[-1,1]への非線形変換であり、
 //セルのサイズをより均一にしようとします。
 // これは、C++ バージョンで「二次変換」と呼ばれるものを使用します。
関数stToUV(s float64) float64 {
    s >= 0.5の場合{
 （1 / 3.）*（4*s*s - 1）を返します
    }
 （1 / 3.）*（1 - 4*（1-s）*（1-s））を返す
} 
 
 // uvToSTはstToUV変換の逆変換です。
 //は  uvToST(stToUV(x)) == x が常に真であるとは限りません。
 // エラー。
関数 uvToST(u float64) float64 {
    u >= 0 の場合 {
 0.5 * math.Sqrt(1+3*u)を返します。
    }
 1 - 0.5*math.Sqrt(1-3*u)を返します
}

補正変換後、u と v は s と t に変換されます。値の範囲も変更されました。 uとvの値の範囲は[-1,1]です。変換後、sとtの値の範囲は[0,1]です。

ここまでをまとめると、球面上の点は S(lat,lng) -> f(x,y,z) -> g(face,u,v) -> h(face,s,t) となります。現在、変換には合計 4 つのステップがあります。球面の緯度と経度の座標が球面の xyz 座標に変換され、次に外接立方体の投影面上の座標に変換され、最後に修正された座標に変換されます。

これまでのところ、S2 が最適化できるポイントは 2 つあります。まず、投影の形状を六角形に変更できるかどうか。次に、変更された変換関数は、計算パフォーマンスに影響を与えないように、実際には tan() に似ているが、計算速度が tan() よりもはるかに速い関数を見つけることができますか。

4. ポイントと軸ポイント間の変換

S2 アルゴリズムでは、セル分割のデフォルトのレベルは 30 です。つまり、正方形は 2^30 * 2^30 の小さな正方形に分割されます。

では、前のステップの s と t をこの正方形にどのように変換すればよいでしょうか?

sとtの値の範囲は[0,1]ですが、値の範囲を[0,230-1]に拡張する必要があります。

 // stToIJ はST 座標の値を IJ 座標の値に変換します。
 stToIJ関数(s float64) int {
戻りクランプ( int (math.Floor(maxSize*s)), 0, maxSize-1)
 }

C++実装でも同様である。

インラインint S2CellId::STtoIJ( double s) {
  // static_castを介してフローティングポイントから整数への変換は非常に遅い
  //丸めモードを変更する必要があるため、Intelプロセッサで。
  // FastIntround（）を使用して最寄りの整数への丸めははるかに高速です。
  //ここで0.5を減算して丸めます
戻る  max （0、 min （kmaxsize -1、mathutil :: fastintround（kmaxsize * s -0.5）））;
 }

このステップでは、h（face、s、t） - > h（face、i、j）です。

5。座標軸ポイントとヒルベルト曲線セルID間の変換

最後のステップは、J、JをHilbert Curveのポイントと関連付ける方法です。

 const（
    ルックアップ= 4
    swapmask = 0x01
    InvertMask = 0x02
 ） 
 
 var（
    ijtopos = [4] [4] int {
        {0、1、3、2}、//正規順序 
        {0、3、1、2}、//軸が交換されます
        {2、3、1、0}、// bits反転
        {2、1、3、0}、//交換および反転
    }
    postoij = [4] [4] int {
        {0、1、3、2}、//正規順序：（0,0）、（0,1）、（1,1）、（1,0）
        {0、2、3、1}、// axes交換：（0,0）、（1,0）、（1,1）、（0,1）
        {3、2、0、1}、// bits反転：（1,1）、（1,0）、（0,0）、（0,1）
        {3、1、0、2}、//交換および逆転：（1,1）、（0,1）、（0,0）、（1,0）
    }
    Postoorientation = [4] int {swapmask、0、0、invertmask |
    lookupij [1 <<（2*lookupbits + 2）] int  
    lookuppos [1 <<（2*lookupbits + 2）] int  
 ）

変換の前に、定義された変数のいくつかを説明しましょう。

Postoijは、いくつかのユニットヒルベルト曲線の位置情報を記録するマトリックスを表します。

以下に示すように、Postoijアレイの情報はグラフで表されます。

同様に、IJTOPOSアレイの情報は、以下に示すようにグラフで表されます。

Postoorientationアレイには、1、0、0、および3の4つの数値が含まれています。

LookupijとLookupposは、それぞれ1024の容量を持つ2つの配列です。これらは、Hilbert曲線IDを座標軸IJに変換する変換表と、座標軸IJをHilbert曲線IDに変換する変換テーブルに対応します。

 func init（）{
    initlookupcell（0、0、0、0、0、0）
    initlookupcell（0、0、0、swapmask、0、swapmask）
    initlookupcell（0、0、0、invertmask、0、invertmask）
    initlookupcell（0、0、0、swapmask | invertmask、0、swapmask | invertmask）
 }

ここに初期化の再帰関数があります。ヒルベルト曲線の標準順序では、4つのグリッドがあり、グリッドが順調であることがわかります。 4番目のパラメーターは0〜3です。

 // initlookupcell lookupijテーブルを初期化します  init時に。​
 func initlookupcell（ level 、i、j、orimorientation、pos、or orientation int ）{ 
 
    level == lookupbitsの場合{
        ij：=（i << lookupbits） + j
        lookuppos [（ij << 2） + orierientation] =（pos << 2） +方向
        lookupij [（pos << 2） + orierientation] =（ij << 2） +方向
     
戻る 
    } 
 
レベル++
    i << = 1
    J << = 1
    pos << = 2 
     
    r：= posoij [オリエンテーション] 
     
    initlookupcell（ level 、i+（r [0] >> 1）、j+（r [0]＆1）、orierientation、pos、or or orientation [0]）
    initlookupcell（ level 、i+（r [1] >> 1）、j+（r [1]＆1）、orierientation、pos+1、方向^郵便通し[1]）
    initlookupcell（ level 、i+（r [2] >> 1）、j+（r [2]＆1）、orierientation、pos+2、orientation^postoorientation [2]）
    initlookupcell（ level 、i+（r [3] >> 1）、j+（r [3]＆1）、orierientation、pos+3、orientation^postoorientation [3]）
 }

上記の関数は、ヒルベルト曲線を生成します。 POS << 2に操作があることがわかります。これにより、位置が最初の4つの小さなグリッドに変換されるため、位置に4を掛けます。

Lookupbits = 4の最初の設定は、[0,15]から、合計16 = 256です。これは、まさにlookupijとlookupposの総容量です。

ローカルグラフを描画します。ここで、J、Jは0〜7で変化します。

上の図は、4次ヒルベルト曲線です。初期化の実際のプロセスは、4次ヒルベルトの1024ポイントの座標と座標軸上のXおよびY軸の間の対応テーブルを初期化することです。

たとえば、次の表は、再帰プロセスにおけるIとJの中間プロセスを示しています。次の表は次のとおりです

行を取り出して、計算プロセスを詳細に分析します。

電流（i、j）=（0,2）、IJの計算プロセスは、Iを4ビットで左にシフトしてJを追加し、全体的な結果を2ビットだけシフトすることを仮定します。目的は、2つの方向位置を離れることです。 IJの最初の4ビットはIです。次の4ビットはJで、最後の2ビットは方向です。この方法で計算されたIJの値は8です。

次に、Lookuppos [IJ]の値を計算します。上記の写真から、（0,2）で表されるセルの4つの数字が16、17、18、および19であることがわかります。計算のこのステップでは、POSの値は4です（POSは生成されたグリッドの数を記録するために使用され、POSの合計値は0から255にサイクリングします）。 POは、現在のグリッドで構成される現在のグリッドを表します。現在のグリッドは4番目のグリッドです。したがって、4*4は、現在のグリッドの最初の数（16）にオフセットできます。現在のグリッドの形状は、Postoijアレイに記録されています。ここからオリエンテーションを取り出します。

上の図を見ると、16、17、18、および19は、Postoij配列軸が回転する場合に対応するため、17は軸回転図の数字1で表されるグリッドにあります。この時点で、オリエンテーション= 1。

このように、Lookuppos [IJ]で表される数は、4*4+1 = 17で計算されます。ここでは、I、J、およびHilbert曲線の数値との対応が完了しました。

では、ヒルベルト曲線の数値を実際の座標にどのように対応するのでしょうか？

lookupijアレイは逆情報を記録します。 Lookupijアレイに保存されている情報とLookupposアレイはまったく逆です。 Lookupijアレイの次の表に保存されている値は、Lookupposアレイの次の表です。 Lookupijアレイを調べてみましょう。Lookupij[17]の値は8で、それに応じて計算されます（i、j）=（0,2）。この時点で、私とJはまだ大きなグリッドです。 Postoijアレイで説明されている形状情報を使用する必要があります。現在の形状は軸回転であり、各座標に4つの小さなグリッドがあることも知っていました。

この時点で、球状座標全体の座標マッピングが完了しました。

ポイントS（lat、lng） - > f（x、y、z） - > g（face、u、v） - > h（face、s、t） - > h（face、i、j） - > cellid。現在、球状の緯度と経度の座標は球状のXYZ座標に変換され、最終的に補正された座標に変換され、座標系が[0,230-1]の間隔をマッピングするために変換されます。

6。S2セルIDデータ構造

最後に、異なるレベルの対応する精度に直接関連するS2セルIDデータ構造について話す必要があります。

S2では、各Cellidは64ビットで構成されています。 UINT64で保存できます。最初の3桁は、キューブの6つの面の1つであり、値範囲[0,5]を表しています。 3ビットは0-7を表すことができますが、6と7は無効な値です。

最後の64桁目は1で、これは特に予約されています。真ん中にあるビットの数をすばやく見つけるために使用されます。最後の最後のビットを楽しみにして、0ではない最初の位置を見つけてください。つまり、最初の1を見つけます。この数字の最初の数字（最初の3桁が占有されているため）の最初の数字は、利用可能な数字です。

緑色のグリッドの数は、分割されたグリッドの数を表すことができます。上の左の写真には、緑に60のグリッドがあるため、[0,230 -1] * [0,230 -1]非常に多くのグリッドを表すことができます。上の右の写真には、緑色のグリッドが36個しかないため、非常に多くのグリッドしか表しない[0,218 -1]*[0,218 -1]。

では、さまざまなレベルで表現できるグリッドの領域の大きさはどれくらいですか？

前の章から、投影のために、投影後の領域にはまだサイズの違いがあることがわかっています。

ここで計算された式は非常に複雑であるため、詳細については、ドキュメントを読むことができます。

 minareametric = metric {2、8 * math.sqrt2 / 9}
 avgareametric = metric {2、4 * math.pi / 6}
 maxareametric = metric {2、2.63579256963161491}

これは、最大面積と最小面積と平均面積との複数の関係です。

レベル0は、キューブの6つの側面の1つです。地球の表面積は510,100,000 km2に等しい。レベル0の面積は、地球の表面積の6分の1です。レベル30が表すことができる最小面積は0.48cm2、最大値は0.93cm2です。

7。S2とGeohashの比較

Geohashには、5000kmから3.7cmの範囲の12レベルがあります。中央の各レベルの変化は比較的大きいです。前のレベルを選択するときにはるかに大きくなる場合があり、次のレベルを選択するときは小さくなる場合があります。たとえば、文字列の長さを4に選択すると、対応するセルの幅は39.1kmで、需要は50kmになる可能性があります。この場合、選択するgeohash文字列を選択することはより困難です。選択が良くない場合は、判断を下すたびに、周囲の8グリッドを取り出して別の判断を下す必要がある場合があります。 Geohashには12バイトのストレージが必要です。

S2には、0.7cm²から85,000,000km²の30レベルがあります。中央の各レベルの変化は比較的滑らかで、4のパワーに曲線に近い。したがって、選択の精度は、Geohashの選択の難易度の問題を引き起こしません。 S2ストレージでは、保存するには1つのUINT64のみが必要です。

S2ライブラリには、ジオコーディングだけでなく、他の多くの幾何学計算関連ライブラリも含まれています。 Geocodeはほんの少しの部分です。この記事では、さまざまなベクター計算、面積計算、ポリゴンカバレッジ、距離の問題、球状球の問題など、この記事では紹介されていないS2の多くの実装があります。

S2は、貪欲なアルゴリズムを使用して、ローカルな最適なソリューションを見つけることもできます。たとえば、都市を考えると、最適なソリューションを見つけ、ポリゴンは都市を覆うだけです。

上の図に示すように、結果のポリゴンは下の青い領域を覆うだけです。ここで生成されたポリゴンは、大小を帯びています。何があっても、最終結果はターゲットをカバーするだけです。

同じセルを使用することで同じ目的を達成できます。

これらはGeohashができないことです。

8。S2セルが例を示します

まず、緯度と経度とCellidの変換、および長方形領域の計算を見てみましょう。

 latlng：= s2.latlngfromdegrees（31.232135、121.41321700000003）
    Cellid：= s2.cellidfromlatlng（latlng）
    セル：= s2.cellfromcellid（cellid）// 9279882742634381312 
 
    //セル。
    fmt.println（ "latlng =" 、latlng）
    fmt.println（ "cell level =" 、cellid。level （ ））
    fmt.printf（ "cell =％d \ n" 、cellid）
    smallcell：= s2.cellfromcellid（cellid.parent（10））
    fmt.printf （ "smallcell level =％d \ n" 、smallcell。level （））
    fmt.printf（ "smallcell id =％b \ n" 、smallcell.id（））
    fmt.printf（ "smallcell ampxarea =％v \ n" 、smallcell.approxarea（）））
    fmt.printf（ "smallcell averagearea =％v \ n" 、smallcell.averagearea（））
    fmt.printf（ "smallcell exactarea =％v \ n" 、smallcell.exactarea（））

ここでは、親メソッドパラメーターは、変更点を返す対応するレベルのCellidを直接指定できます。

上記の方法では、結果を次のように印刷します。

 latlng = [31.2321350、121.4132170]
細胞レベル= 30
 Cell = 3869277663051577529 
 
 ****Parent **** 10000000000000000000000000000000000000000
 Smallcellレベル= 10
 smallCell id = 11010110110010011011110000000000000000000000000000000000000000
 Smallcell Ampxarea = 1.9611002454714756E-06
 Smallcell Averagearea = 1.997370817559429E-06
 Smallcell Exactarea = 1.9611009480261058E-06

ポリゴンをカバーする別の例を見てみましょう。最初にゾーンを作成しましょう。

 rect = s2.rectfromlatlng（s2.latlngfromdegrees（48.99、1.852）））
    rect = rect.addpoint（s2.latlngfromdegrees（48.68、2.75）） 
 
    rc：=＆s2.regioncoverer {maxlevel：20、maxcells：10、minlevel：2}
    r：= s2.region（rect.capbound（））
    カバー：= rc.covering（r）

オーバーレイパラメーターはレベル2〜20に設定され、セルの最大数は10です。

次に、セルを最大20に変更します。

最後に、30に変更されました。

同じレベルの範囲を見ることができ、セルが多いほど、ターゲット範囲がより正確になります。

一致する長方形の領域は次のとおりです。また、円形の領域の一致にも同じことが言えます。

コードは投稿されず、長方形に似ています。 Geohashはこの機能を行うことができず、自分で手動で実装する必要があります。

9。S2の適用

[[201839]]

現在、S2はより頻繁に使用されており、マップ関連のビジネスで使用されています。 Google MapはS2を大量に直接使用するため、読者は自分でそれを体験できます。 UberはS2アルゴリズムを使用して、最寄りのタクシーを計算します。シーンの例は、この記事の引用で言及されているシーンです。 Didiも関連するアプリケーションを持っている必要があり、おそらくより良い解決策があるかもしれません。これらの空間インデックスアルゴリズムは、現在非常に人気のある共有自転車でも使用されます。

最後に、持ち帰り業界は地図と密接に関連しています。 MeituanとEle.meは、この分野には多くのアプリケーションがあると考えています。

5つ

[[201840]]

この記事では、GoogleのS2アルゴリズムの基本的な実装に焦点を当てています。 Geohashは空間点インデックスアルゴリズムでもありますが、そのパフォーマンスはGoogleのS2よりもわずかに劣っています。さらに、大企業のデータベースは基本的にGoogleのS2アルゴリズムを使用してインデックスを作成し始めています。

実際には、空間検索、多次元の空間表面、および多次元の空間ポリゴンを検​​索する方法があります。街路、背の高い建物、鉄道、川などの実用的な例。これらのものを検索するには、データベーステーブルを設計するにはどうすればよいですか？

答えはもちろん、高効率検索も達成できるため、Rツリー、またはRツリーとB+ツリーが必要です。

この部分は、この記事の範囲内ではありません。次回は「多次元空間ポリゴンインデックスアルゴリズム」を共有できます。

最後に、アドバイスをください。

[この記事は、51CTOコラムニスト「Halfrost」の元の原稿です。

この著者の他の記事を読むにはここをクリックしてください