S2算法应用

需求：计算不同区域范围，X公里半径内实体店或场站覆盖率。

实现思路：

---

• 为了便于理解，将地球看成一个基于经纬度线的坐标系。将经度和纬度看成二维坐标系中的两个纬度，横轴表示经度[-180o, 0o)，（0o, 180o]，纵轴表示纬度[-90o, 0o)，（0o, 90o]。
• 以最小纬度和经度对应坐标为第一个六角形中心点，在经度方向循环计算六角形各顶点（顶点开始，顺时针，命名：Point1,Point2,Point3,Point4,Point5,Point6)及中心点(Point0)坐标，直至六角形中心点经度大于等于最大经度。 存储六角形标记为（0,0),(0,1),(0,2) ......., 表示六边形位于第0行第N列
• 第一行计算完成后，开始计算第2行，如下图，第2行（奇数行）, 除第1个和最后一个点作特殊处理外，其它点的 Point3、Point4、Point5是重叠的，注意坐标的处理，否则会出现偏差
• 
• 在维度方向上按上一步循环。完成整个区域范围内六角边分割，注意：为了快速定位，还需计算出每个六角形中心点对应的Geohash，根据半径不一致，可选择不同的Geohash级别
• 根据场站经纬度计算出30级 Cell ID 值 及 Geohash 码
• 根据 Geohash 码找出附近的六角形，通过六顶点坐标 构造 IRegion, 判断场站是否包含在六边形内，如果不包含，再次计算出当前Geohash码周边8个Geohash框，再次计算
• 至此完成此区域内场站命中的六角形。

关键代码

---

• 根据中心点坐标、边长、偏差角度（中心点至顶点开始）计算下一个点坐标

•         public const double Ea = 6378137;     //   赤道半径（米）  
          public const double Eb = 6356725;     //   极半径 （米）
  
          /// <summary>
          /// 
          /// </summary>
          /// <param name="lat"></param>
          /// <param name="lng"></param>
          /// <param name="distance"></param>
          /// <param name="angle"></param>
          /// <returns></returns>
          public static Point GetPoint(double lat, double lng, double distance, double angle)
          {
  
              double dx = distance * 1000 * Math.Sin(angle * Math.PI / 180.0);
              double dy = distance * 1000 * Math.Cos(angle * Math.PI / 180.0);
  
              double ec = Eb + (Ea - Eb) * (90.0 - lat) / 90.0;
              double ed = ec * Math.Cos(lat * Math.PI / 180.0);
  
              double newLon = (dx / ed + lng * Math.PI / 180.0) * 180.0 / Math.PI;
              double newLat = (dy / ec + lat * Math.PI / 180.0) * 180.0 / Math.PI;
  
              return new Point(newLat, newLon);
          }

   

• 计算场站所属六边形

• /// <summary>
          /// 
          /// </summary>
          /// <param name="destRows"></param>
          /// <param name="cell"></param>
          /// <param name="level"></param>
          /// <param name="staid"></param>
          /// <param name="hashMap"></param>
          /// <param name="geohashValue"></param>
          /// <returns></returns>
          private string GetPgID(DataTable dest, S2Cell cell,string geohashValue)
          {
              //先找当前geohash4的值
              string pgID = this.GetPGIDByHash(dest, cell, geohashValue);
              if (string.IsNullOrEmpty(pgID) == false)
                  return pgID;
              
  
              //当前hash未命中时，找相邻8格
              List<string> hashLst = GeoHashService.Default.GetGeoHashExpand(geohashValue);
              foreach (string ghValue in hashLst)
              {
                  pgID = this.GetPGIDByHash(dest, cell, ghValue);
                  if (string.IsNullOrEmpty(pgID) == false)
                      return pgID;
              }
  
              return string.Empty;
          }
          
         /// <summary>
         /// 
         /// </summary>
         /// <param name="dest"></param>
         /// <param name="cell"></param>
         /// <param name="geohashValue"></param>
         /// <returns></returns>
          private string GetPGIDByHash(DataTable dest, S2Cell cell, string geohashValue)
          {
              DataRow[] destRows = dest.Select(string.Format("{0} = '{1}'", M_GEOHASH, geohashValue));  //城市均分的网格
  
              foreach (DataRow dRow in destRows)
              {
                  string pgID = Convert.ToString(dRow["ID"]);
                  IS2Region cells = this.BuildPolygon(dRow);
  
                  if (cells.Contains(cell) == true)
                  {
                      return pgID;
                  }
              }
  
              return string.Empty;
          }
  
          /// <summary>
          /// 构造容器
          /// </summary>
          /// <param name="row"></param>
          /// <returns></returns>
          private IS2Region BuildPolygon(DataRow row)
          {
              List<S2Point> lst = new List<S2Point>();
              lst.Add(this.GetPoint(Convert.ToDouble(row["lat1"]), Convert.ToDouble(row["lng1"])));
              lst.Add(this.GetPoint(Convert.ToDouble(row["lat2"]), Convert.ToDouble(row["lng2"])));
              lst.Add(this.GetPoint(Convert.ToDouble(row["lat3"]), Convert.ToDouble(row["lng3"])));
              lst.Add(this.GetPoint(Convert.ToDouble(row["lat4"]), Convert.ToDouble(row["lng4"])));
              lst.Add(this.GetPoint(Convert.ToDouble(row["lat5"]), Convert.ToDouble(row["lng5"])));
              lst.Add(this.GetPoint(Convert.ToDouble(row["lat6"]), Convert.ToDouble(row["lng6"])));
  
              S2Loop loop = new S2Loop(lst);
              loop.Normalize();
              return loop;
          }

示例效果

---

•  

 

参考资料

---

• S2演示地址：http://s2map.com/#order=latlng&mode=polygon&s2=false&points=