C#实现Google S2算法

OStack程序员社区-中国程序员成长平台 › 门户 › 编程› c#›C#教程

原作者: [db:作者] 来自: [db:来源] 收藏邀请

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 的实现还有很多很多，各种向量计算，面积计算，多边形覆盖，距离问题，球面球体上的问题，它都有实现。