S2其實是來自幾何數學中的一個數學符號 S²，它表示的是單位球。S2 這個庫其實是被設計用來解決球面上各種幾何問題的。值得提的一點是，除去 golang 官方 repo 里面的 geo/s2 完成度目前只有40%，其他語言，Java，C++，Python 的 S2 實現都完成100%了。看看怎么用 S2 來解決多維空間點索引的問題。通常地球上的點我們會用經緯度來表示，將經緯度坐標轉換為希爾伯特曲線 Cell ID包含如下步驟：
- 把經緯度轉換成弧度。
- 經緯弧度轉換成坐標系上的一個點 S(lat,lng) -> f(x,y,z)
- 球面變平面 f(x,y,z) -> g(face,u,v)
- 球面矩形投影修正:g(face,u,v) -> h(face,s,t)
- 點與坐標軸點相互轉換:h(face,s,t) -> H(face,i,j)
- 坐標軸點與希爾伯特曲線 Cell ID 相互轉換:H(face,i,j) -> CellID
1. 球面坐標轉換 S(lat,lng) -> f(x,y,z)
按照之前我們處理多維空間的思路，先考慮如何降維，再考慮如何分形。眾所周知，地球是近似一個球體。球體是一個三維的，如何把三維降成一維呢？球面上的一個點，在直角坐標系中，可以這樣表示：
-
x = r * sin θ * cos φ y = r * sin θ * sin φ z = r * cos θ - 通常地球上的點我們會用經緯度來表示。
- 再進一步，我們可以和球面上的經緯度聯系起來。不過這里需要注意的是，緯度的角度 α 和直角坐標系下的球面坐標 θ 加起來等於 90°。所以三角函數要注意轉換。於是地球上任意的一個經緯度的點，就可以轉換成 f(x,y,z)。
- 在 S2 中，地球半徑被當成單位 1 了。所以半徑不用考慮。x，y，z的值域都被限定在了[-1,1]這個區間之內了。
- 示例：坐標 （36.683, 117.1412）
- 第一步：把經緯度轉換成弧度。由於經緯度是角度，弧度轉角度乘以 π / 180°
- 第二步：經緯弧度轉換成坐標系上的一個點 S(lat,lng) -> f(x,y,z)
2. 球面變平面 f(x,y,z) -> g(face,u,v)
接下來一步 S2 把球面碾成平面。首先在地球外面套了一個外切的正方體，如下圖。
- 從球心向外切正方體6個面分別投影。S2 是把球面上所有的點都投影到外切正方體的6個面上。
- 這里簡單的畫了一個投影圖，上圖左邊的是投影到正方體一個面的示意圖，實際上影響到的球面是右邊那張圖。
- 從側面看，其中一個球面投影到正方體其中一個面上，邊緣與圓心的連線相互之間的夾角為90°，但是和x，y，z軸的角度是45°。我們可以在球的6個方向上，把45°的輔助圓畫出來，見下圖左邊。
-
上圖左邊的圖畫了6個輔助線，藍線是前后一對，紅線是左右一對，綠線是上下一對。分別都是45°的地方和圓心連線與球面相交的點的軌跡。這樣我們就可以把投影到外切正方體6個面上的球面畫出來，見上圖右邊。投影到正方體以后，我們就可以把這個正方體展開了。
- 一個正方體的展開方式有很多種。不管怎么展開，最小單元都是一個正方形。 以上就是 S2 的投影方案。
- 接下來講講六邊形的投影方案。按照六邊形來投影可能是目前最好的方式，不過也可能是最復雜的方式。
- 六邊形的菱角比較少，六個邊也能相互銜接其他的六邊形。看上圖最后邊的圖可以看出來，六邊形足夠多以后，非常近似球體。
-
六邊形展開以后就是上面這樣。當然這里只有12個六邊形。六邊形個數越多越好，粒度越細，就越貼近球體。Uber 在一個公開分享上提到了他們用的是六邊形的網格，把城市划分為很多六邊形。這塊應該是他們自己開發的。也許滴滴也是划分六邊形，也許滴滴有更好的划分方案也說不定。
-
回到 S2 上面來，S2是用的正方形。這樣第一步的球面坐標進一步的被轉換成 f(x,y,z) -> g(face,u,v)，face是正方形的六個面，u，v對應的是六個面中的一個面上的x，y坐標。
3. 球面矩形投影修正:g(face,u,v) -> h(face,s,t)
上一步我們把球面上的球面矩形投影到正方形的某個面上，形成的形狀類似於矩形，但是由於球面上角度的不同，最終會導致即使是投影到同一個面上，每個矩形的面積也不大相同。
- 上圖就表示出了球面上個一個球面矩形投影到正方形一個面上的情況。
-
經過實際計算發現，最大的面積和最小的面積相差5.2倍。見上圖左邊。相同的弧度區間，在不同的緯度上投影到正方形上的面積不同。現在就需要修正各個投影出來形狀的面積。如何選取合適的映射修正函數就成了關鍵。目標是能達到上圖右邊的樣子，讓各個矩形的面積盡量相同。
-
線性變換是最快的變換，但是變換比最小。
-
tan() 變換可以使每個投影以后的矩形的面積更加一致，最大和最小的矩形比例僅僅只差1.414。可以說非常接近了。但是 tan() 函數的調用時間非常長。如果把所有點都按照這種方式計算的話，性能將會降低3倍。
-
最后谷歌選擇的是二次變換，這是一個近似切線的投影曲線。它的計算速度遠遠快於 tan() ，大概是 tan() 計算的3倍速度。生成的投影以后的矩形大小也類似。不過最大的矩形和最小的矩形相比依舊有2.082的比率。
-
上表中，ToPoint 和 FromPoint 分別是把單位向量轉換到 Cell ID 所需要的毫秒數、把 Cell ID 轉換回單位向量所需的毫秒數。
-
（Cell ID 就是投影到正方體六個面，某個面上矩形的 ID，矩形稱為 Cell，它對應的 ID 稱為 Cell ID）。ToPointRaw 是某種目的下，把 Cell ID 轉換為非單位向量所需的毫秒數。在 S2 中默認的轉換是二次轉換。
- 投影之后的修正函數三種變換應該如下：
- 注意上面雖然變換公式只寫了u，不代表只變換u。實際使用過程中，u，v都分別當做入參，都會進行變換。
- 經過修正變換以后，u，v都變換成了s，t。值域也發生了變化。u，v的值域是[-1,1]，變換以后，是s，t的值域是[0,1]。
- S2 可以優化的點有兩處，一是投影的形狀能否換成六邊形？二是修正的變換函數能否找到一個效果和 tan() 類似的函數，但是計算速度遠遠高於 tan()，以致於不會影響計算性能？
4. 點與坐標軸點相互轉換:h(face,s,t) -> H(face,i,j)
在 S2 算法中，默認划分 Cell 的等級是30，也就是說把一個正方形划分為 2^30 * 2^30個小的正方形。那么上一步的s，t映射到這個正方形上面來，對應該如何轉換呢？
-
- s，t的值域是[0,1]，現在值域要擴大到[0,2^30^-1]。
5. 坐標軸點與希爾伯特曲線 Cell ID 相互轉換:H(face,i,j) -> CellID
最后一步，把 i，j 和希爾伯特曲線上的點關聯起來，如下畫一個局部的圖，i，j從0-7變化
- 推理過程略：參見官方PPT: https://docs.google.com/presentation/d/1Hl4KapfAENAOf4gv-pSngKwvS_jwNVHRPZTTDzXXn6Q/view#slide=id.i22
6. S2 Cell ID 數據結構
S2 Cell ID 數據結構，這個數據結構直接關系到不同 Level 對應精度的問題。如下圖：
- 在 S2 中，每個 CellID 是由64位的組成的。可以用一個 uint64 存儲。開頭的3位表示正方體6個面中的一個，取值范圍[0,5]。3位可以表示0-7，但是6，7是無效值。
- 64位的最后一位是1，這一位是特意留出來的。用來快速查找中間有多少位。從末尾最后一位向前查找，找到第一個不為0的位置，即找到第一個1。這一位的前一位到開頭的第4位（因為前3位被占用）都是可用數字。
- 綠色格子有多少個就能表示划分多少格。上圖左圖，綠色的有60個格子，於是可以表示[0,2^30^ -1] [0,2^30^ -1]這么多個格子。上圖右圖中，綠色格子只有36個，那么就只能表示[0,2^18^ -1][0,2^18^ -1]這么多個格子。
- level 0 就是正方體的六個面之一。地球表面積約等於510,100,000 km^2^。level 0 的面積就是地球表面積的六分之一。level 30 能表示的最小的面積0.48cm^2^，最大也就0.93cm^2^ 。
7. C#示例 代碼
-
private static void GetCellDetail(double lat,double lng) { Console.WriteLine(string.Format("-------------坐標： 36.683，117.1412---------------- ")); /*第一步：把經緯度轉換成弧度。由於經緯度是角度，弧度轉角度乘以 π / 180° */ S2LatLng ll = S2LatLng.FromDegrees(36.683, 117.1412); Console.WriteLine(string.Format("-------------第一步：把經緯度轉換成弧度--------------")); Console.WriteLine(string.Format("latr:{0},lngr{1}", ll.LatRadians, ll.LngRadians)); //第二步：經緯弧度轉換成坐標系上的一個點 S(lat,lng) -> f(x,y,z) */ S2Point point = ll.ToPoint(); Console.WriteLine(string.Format("-------------第二步：經緯弧度轉換成坐標系上的一個點 S(lat,lng) -> f(x,y,z) --------------")); Console.WriteLine(string.Format("f(x,y,z): x:{0},y:{1},z:{2}", point.X, point.Y, point.Z)); //第三步：進行投影 f(x,y,z) -> g(face,u,v) int face = S2Projections.XyzToFace(point); R2Vector vector = S2Projections.ValidFaceXyzToUv(face, point); Console.WriteLine(string.Format("-------------進行投影 f(x,y,z) -> g(face,u,v) --------------")); Console.WriteLine(string.Format("g(face,u,v): {0},{1},{2}", face, vector.X, vector.Y)); //第四步： 投影面積修正 g(face,u,v) -> h(face,s,t) double s = S2Projections.UvToSt(vector.X); double t = S2Projections.UvToSt(vector.Y); Console.WriteLine(string.Format("-------------第四步： 投影面積修正 g(face,u,v) -> h(face,s,t) --------------")); Console.WriteLine(string.Format("h(face,s,t): {0},{1},{2}", face, s, t)); //第五步：點與坐標軸點相互轉換 h(face,s,t) -> H(face,i,j) int i = S2CellId.StToIj(s); int j = S2CellId.StToIj(t); Console.WriteLine(string.Format("-------------第五步：點與坐標軸點相互轉換 h(face,s,t) -> H(face,i,j) --------------")); Console.WriteLine(string.Format("H(face, i, j): {0},{1},{2}", face, i, j)); //第六步：坐標軸點與希爾伯特曲線 Cell ID 相互轉換 S2CellId cellid = S2CellId.FromFaceIj(face, i, j); Console.WriteLine(string.Format("-------------第六步：坐標軸點與希爾伯特曲線 Cell ID 相互轉換 --------------")); Console.WriteLine(string.Format("CellID.id: {0},level:{1}", cellid.Id, cellid.Level)); //驗證轉坐標 S2LatLng lan = cellid.ToLatLng(); Console.WriteLine(string.Format("lat: {0},lng:{1}", lan.LatDegrees, lan.LngDegrees)); }
8、S2 與 Geohash 對比
- Geohash 有12級，從5000km 到 3.7cm。中間每一級的變化比較大。有時候可能選擇上一級會大很多，選擇下一級又會小一些。比如選擇字符串長度為4，它對應的 cell 寬度是39.1km，需求可能是50km，那么選擇字符串長度為5，對應的 cell 寬度就變成了156km，瞬間又大了3倍了。這種情況選擇多長的 Geohash 字符串就比較難選。選擇不好，每次判斷可能就還需要取出周圍的8個格子再次進行判斷。Geohash 需要 12 bytes 存儲。
- S2 有30級，從 0.7cm² 到 85,000,000km² 。中間每一級的變化都比較平緩，接近於4次方的曲線。所以選擇精度不會出現 Geohash 選擇困難的問題。S2 的存儲只需要一個 uint64 即可存下。
- S2 庫里面不僅僅有地理編碼，還有其他很多幾何計算相關的庫。地理編碼只是其中的一小部分。本文沒有介紹到的 S2 的實現還有很多很多，各種向量計算，面積計算，多邊形覆蓋，距離問題，球面球體上的問題，它都有實現。
參考資料