谷歌的S2,球体上的几何体,细胞和希尔伯特曲线

更新日期:2017年12月5日:谷歌刚刚宣布它将致力于S2库新版本的开发,好消息,可以找到存储库在这里.

谷歌的s2库是真正的财富,这不仅是因为它的空间索引功能,还因为它是一个4年前发布的库,没有得到应有的关注。s2库由谷歌自己在谷歌地图上使用,MongoDB引擎和Foursquare,但是除了Foursquare上的一篇论文外,你在任何地方都找不到关于图书馆的任何文档或文章,一个谷歌演示以及源代码注释。您还将很难找到库的绑定,官方存储库缺少Python库的Swig文件,这要感谢一些叉子我们可以对python语言进行部分绑定(我将在这篇文章中使用它)。我听说Google现在正在积极地开发这个库,我们很可能很快会在他们发布这篇文章时获得更多关于它的细节。但是我决定分享一些关于这个库的例子,以及我认为这个库很酷的原因。

通往牢房的路

您将在s2代码中看到这个“单元”概念。细胞是球体(地球在我们的情况下,但你并不局限于它)到区域或点的紧凑表示。区域也可以用这些相同的单元近似,有一些很好的特点:

  • 它们是紧凑的(由64位整数表示)
  • 他们对地理特征有分辨力
  • 他们是有等级的(他们有等级,相似的水平面也有相似的区域)
  • 任意区域的包含查询非常快

S2库首先将球面的点/区域投影到一个立方体中,立方体的每个面都有一个四叉树球面投影到这里。之后,发生了一些转换(有关原因的详细信息,见谷歌演示)空间是离散的,之后,在希尔伯特曲线,这就是为什么这个图书馆这么好,希尔伯特曲线是一种空间填充曲线,它将多个维度转换为一个具有特殊空间特征的维度:它保留了位置。

希尔伯特曲线

希尔伯特曲线

希尔伯特曲线是空间填充曲线,这意味着它的范围覆盖了整个n维空间。为了理解这是怎么回事,你可以想象一根长的绳子,它以一种特殊的方式排列在空间上,这样绳子就可以穿过空间的每一个正方形。这样就填满了整个空间。要沿希尔伯特曲线转换二维点,您只需要选择字符串中该点所在的点。理解它的一个简单方法是使用这个迭代的例子在这里你可以点击曲线上的任何点,它会显示点在字符串中的位置,反之亦然。

在下图中,希尔伯特曲线最开始处的点(字符串)bepaly亚洲也位于沿着曲线的最开始处(曲线由图像底部的长字符串表示):

希尔伯特曲线
希尔伯特曲线

在下图中我们有更多的点,很容易看出希尔伯特曲线是如何保持空间局部性的。你可以注意到曲线上的点彼此更靠近(在一维表示法中,在二维空间(x,y平面)中,底部的这条线也更近。然而,请注意,相反的情况并不完全正确,因为在X,Y平面上可以有彼此靠近的二维点,而这些点在希尔伯特曲线上不是很近。

希尔伯特曲线
希尔伯特曲线

因为s2使用hilbert曲线来枚举单元格,这意味着靠近值的单元值在空间上也彼此接近。当这种思想与层次分解相结合时,对于索引和查询操作,您bepaly亚洲有一个非常快速的框架。在我们开始实际例子之前,让我们看看这些单元格是如何用64位整数表示的。

如果你对希尔伯特曲线感兴趣,我真的很推荐这篇文章,它非常直观,bepaly亚洲显示了曲线的一些特性。

细胞表达

正如我之前提到的,这些细胞有不同的水平和它们可以覆盖的不同区域。在S2库中,您将发现30个层次分解级别。我在下面复制的幻灯片中的谷歌演示文稿中显示了它们可以覆盖的单元格级别和区域范围:

单元格区域
单元格区域

如你所见,s2bepaly亚洲几何的一个非常酷的结果是,地球上的每厘米平方都可以用64位整数来表示。

使用以下模式表示单元格:

单元表示模式
单元格表示模式(来自原始Google演示文稿的图像)

第一个代表叶细胞,具有最小面积的单元,通常用来表示点。如你所见,这3个初始位保留用于存储投影球面的立方体的面,接着是希尔伯特曲线中单元格的位置,后面总是跟着一个“1”位,这是一个标识单元格级别的标记。

所以,为了检查细胞的水平,所需要做的就是检查最后一个“1”位在单元格表示中的位置。检查安全壳,为了验证一个单元是否包含在另一个单元中,你只需要做一个前缀比较。这些操作真的很快,它们是可能的,因为希尔伯特曲线枚举和层次分解方法使用。

覆盖区域

所以,如果你想生成覆盖一个区域的单元格,可以使用库中指定最大单元格数的方法,将使用的最大单元级别和最小单元级别以及算法将使用指定参数近似该区域。在下面的例子中,我用S2库来提取一些机器学习二进制特征使用15级细胞:

15级单元—机器学习的二进制特性
15级单元—机器学习的二进制特性

单元格区域在上图中使用透明多边形表示,覆盖了我所在城市的整个区域。因为我在最小和最大级别都使用了level 15,这些细胞都覆盖着相似的区域。如果我把最小值改为8(这样就可以使用更大的单元格)。该算法将近似该区域的方式,它将提供更少的单元格,并试图保持近似精度,如下面的例子:

覆盖范围由8至15(级别)
覆盖范围由8至15(级别)

如你所见,我们现在有一个覆盖物,在中心使用较大的单元格,为了处理边界,我们有一个近似值,使用较小的单元格(还要注意四叉树)。

例子

*在本教程中,我使用了在库.编译和安装它的说明出现在存储库的自述文件中,所以我不会在这里重复它。

将纬度/经度点转换为单元表示的第一步如下所示:

>>>导入s2>>latlng=s2.s2 latlng.fromdegrees(-30.043800,>>> cell = s2.S2CellId.FromLatLng(latlng)>>> cell.level()30>>> cell.id()10743750136202470315> >0 >1 cell. totoken ()951977d377e723ab

如你所见,我们首先创建一个类的对象S2LatLng表示lat/lng点,然后将其输入S2CellId类来构建单元格表示。之后,我们可以得到该类的级别和ID。还有一个方法叫做托托肯它将整数表示转换为紧凑的字母数字表示,您可以稍后使用该表示进行解析来自令牌方法。

您还可以获得该单元格的父单元格(比它高一级),并使用包含方法检查某个单元格是否由另一个单元格包含:

>>>parent=cell.parent()>>>打印parent.level()29>>parent.id()1074375013622047316>>parent.totooken()951977d377e723ac>>cell.contains(parent)false>>parent.contains(cell)true

如你所见,父代的级别是子代单元之上的级别(在我们的例子中,一片树叶细胞)。ID也非常相似,除了单元的级别和包bepaly亚洲含检查非常快(它只检查父单元的子单元的范围)。

这些单元格可以存储在数据库中,并且在BTree索引中执行得非常好。为了创建覆盖一个区域的单元格集合,您可以使用S2RegionCoverer类的例子如下:

>>>区域矩形=s2latlng rect(s2latlng.fromdegrees(-51.264871,-30.241701),S2LatLng.FromDegrees(-51.04618,-30.000003))>>>coverer=s2regioncoverer()>>>coverer.set_min_level(8)>>>coverer.set_max_level(15)>>>coverer.set_max_cells(500)>>>covering=coverer.getcovering(region_rect)

首先,我们定义了S2LatLngRect它是一个矩形,定义了我们要覆盖的区域。还可以使用其他类(例如,构建多边形)。的S2RegionCoverer与使用S2区类作为基类。定义矩形后,我们实例化S2RegionCoverer然后设置前面提到的最小/最大电平以及我们想要近似生成的最大单元格数。

如果你想绘制封面,你可以用卡通画,美观和matplotlib,如下面的示例所示:

将matplotlib.pyplot导入为plt from s2 import*from shapely.geometry import polygonimport cartopy.crs导入为ccrsimport cartopy.io.img_tiles as cimgtproj=cimgt.mapquestosm()plt.figure(FigSize=(20,20),dpi = 200) ax = plt.axes(投影= proj.crs) ax.add_image(项目,12)最大设定范围([-51.411886,-50.922470,-30.301314,-29.94364])区域矩形=s2latlng rect(s2latlng.fromdegrees(-51.264871,-30.241701),S2LatLng.FromDegrees(-51.04618,-30.000003))coverer=s2regioncoverer()coverer.set_min_level(8)coverer.set_max_level(15)coverer.set_max_cells(500)covering=coverer.getcovering(region rect)geoms=[]for cellind in covering:new_cell=s2cell(cellind)vertices=[]for i in xrange(0,4):vertex=new_cell.getVertex(i)latlng=s2latlng(vertex)vertices.append((latlng.lat().degrees(),geo =多边形(顶点)geoms.append(geo)print "Total Geometries: {}".format(len(geoms)) ax.add_geometries(geoms,platecarree()固定电流调节器,facecolor='coral',EdgeColor='Black',alpha=0.4)plt.show()。

结果如下:

覆盖细胞。
覆盖细胞。

S2 API中有很多东西,我真的建议你去探索和阅读源代码,这真的很有帮助。S2单元格可用于索引和键值数据库,它可以用在B树上,效率很高,甚至可以用于机器学习(我的例子)不管怎样,它是一个非常有用bepaly亚洲的工具,你应该把它放在工具箱里。希望你喜欢这个小教程!

——基督教。佩隆

克里斯蒂安S佩隆

58岁的评论

  1. 它是四叉树投射到单位球面上,加上使用希尔伯特曲线快速映射整数索引和单元格。

  2. 感谢这篇文章。帮了我不少忙。你在StackExchange上发布了我问题的链接。

    我已经运行了用Java获取CELLID的例子。它运作良好。然而,我有一些问题,你也许可以回答或给同样的方向:

    1) cell ID与geo-hash有什么不同?还是一样?

    2)通过S2可以进行路由吗?

    3)我在用已有的车辆/自行车数据进行ETA(到达时间)分析时,偶然发现了S2库。我读过谷歌,foursquare,Uber使用这个图书馆。我贴出了同样的问题:gis。课件。com/questions/152719/postgis用于gpx分析。几天过去了。你能想出解决这个问题的方法吗?

    1. 你好,我会说s2是一种具有不同属性的geohash。无法使用S2进行路由,S2不是一个路由算法而是球面上的几何,它可以用于生成路由的覆盖,但不能用于路由本身,为此,我建议您研究传统的路由算法,比如A-Star,等。
      你的问题是一个很开放的问题,你可以使用机器学习的方法,利用历史记录来估算每天特定时间到达某个目的地所需要的时间,你可以使用很多方法,但我真的不知道不同模型之间的(结果绩效)基准,试着找相关的文件,你当然应该找一些与你的问题有关的东西。

      1. 你说得太好了——“它可以用来为路线生成覆盖物,但不能用于路由本身,”
        这意味着,如果我在一条路线上有一组点,我可以用S2单元格创建一个覆盖层。
        我必须根据历史数据预先计算,在一天中的某个特定时间移动单元格所需的时间。我可以把它们加起来得到。

        但是,希尔伯特曲线的最佳用途是什么?我仍然对同样的W.R.T用例感到困惑。s2细胞是否能获得更多的信息?就像我的时间问题?希尔伯特曲线也许能让我更快地查找这些信息?

        1. 不确定问题是否还在你的脑海中,但如果是:
          geohash和s2单元的制作目的相同,即在保留位置的同时,将球形地理空间平面划分为分区/单元(意味着两个地理空间相邻的单元应具有“相似”的哈希代码/单元ID)。

          GeoHash将球面分成两部分,一个代码为1,另一个代码为0(lvl1)。采用其中一个部分并继续相同的过程将导致新的LVL具有更细的粒度和更长的geohash代码。
          如您所见,这种方法在嵌入过程中是线性的(lon,lat->geohash代码),这使得它非常快速,但缺乏准确性(主要是在赤道线附近),因此在单元大小上有很大的差异(两个相邻单元在单元区域中可以有很大的差异)。

          S2cell嵌入更复杂,用立方体包围球面,用四叉树将每个立方体面分成4个部分(每个部分分成更小的部分,由此产生了s2单元级别的概念,即树中的级别),最后使用hilbert空间填充曲线将树叶枚举为s2单元ID(长代码)。
          因此,嵌入是二次(四叉树),比geohash慢一点,但考虑到从球体到立方体的投影,该立方体相对于单元区域的方差较小(相邻单元的大小更接近,因此,当将每个单元格中的外观聚合起来并相互比较时,错误会减少——这是一个“热点”吗?bepaly亚洲

          希望这可以帮助!
          s2电池棒极了,这是一个伟大的作品!
          (这里提供的所有数据都是“公共”的,这意味着您可以通过上面的演示文稿和源代码进行一些推断:))

          另一个可能有帮助的资源是这个演示:
          http://www.slideshare.net/danielmarcous/geo-data-analytics

          当做,
          丹尼尔·马库斯。

  3. 然而,注意,相反的情况并不完全正确因为在xy平面上可以有两个相互接近的点而在希尔伯特曲线上却不接近"

    我注意到了这一点,想知道这意味着什么。这似乎会使一些行动“不对称”。这在实践中是一个问题吗?

  4. “除了一张四方形的纸,你在任何地方都找不到关于图书馆的任何文档或文章,谷歌演示文稿和源代码注释。

    这就跟你问声好!我在哪里可以找到四方形的纸?我在java版本的S2CellId类中搞混了,它如何递归地生成查找表,位操作太多,代码很难读。

  5. 你好,Christian -我想micolous库的库把坐标映射到了错误的立方体面,导致更像钻石的细胞形状。我用华盛顿特区的坐标做了个掩护。Micolous图书馆把它们放在101号脸上,而s2map把它们放在100号脸上。我在这里贴了一张对比幻灯片:http://imgur.com/ymvw3o2

  6. 基督徒,这是一篇很棒的文章。我有一个问题-s2单元是否位于理想化球体上?或者它反映了实际表面的轮廓?因为它是球面,是否因为地球的实际形状是一个扁椭球体而存在问题?

  7. 这是奇妙的。谢谢克里斯蒂安,我很抱歉找不到一个比“神奇”更能称颂你的词。

  8. s2可以用来进行geo-radial查询吗?例如,说,我有一套lat/long,我想知道在给定的距离内,它们中的哪些在给定的位置。基本上是geoadd和georadis redis命令的组合。

  9. 伟大的职位。

    如果数据库中有一个lat/long数据负载,那么如何创建s2函数,使其从lat/long转到给定级别的单元格_ID。

    我在谷歌bigquery中有lat/long数据,我想知道我是否可以使用UDF来获得BQ中的s2功能,从而避免将数据发回python或其他地方进行处理。

    干杯
    安迪

  10. 好的文章,Python 3.5和s2sphere (Python)的更新

    将s2sphere导入为s2

    将matplotlib.pyplot导入为plt
    从shapely.geometry导入多边形
    将cartopy.crs导入为ccr
    将cartopy.io.img_图块导入为cimgt

    项目= cimgt.OSM ()
    plt.图(图尺寸=(20,20),DPI=200)
    ax=plt.axs(投影=proj.crs)
    ax.add_image(项目,12)

    ax.set_extent([-51.411886,-50.922470,
    -30.301314,-29.94364])

    区域矩形=s2.latlngrect(
    s2.LatLng.from_degrees(-51.264871,-30.241701),
    s2.LatLng.from_degrees(-51.04618,-30.000003))
    包装工人= s2.RegionCoverer ()
    coverer.min_级别=8
    coverer.max_level = 15
    coverer.max_cells = 500
    覆盖= coverer.get_covering(region_rect)
    GeOMS=

    对于cellid的覆盖:
    new_cell = s2.Cell (cellid)
    顶点=[]
    对于区间(0)内的i,4)中:
    顶点=新的单元。获取顶点(i)
    latlng=s2.latlng.from_点(顶点)
    vertices.append (latlng.lat () .degrees latlng.lng () .degrees))
    geo=多边形(顶点)
    geoms.append(geo)
    打印(“总几何图形:”.格式(len(geoms)))

    ax.add_geometries(几何学,platecarree()固定电流调节器,facecolor='coral',EdgeColor='Black',α=0.4)
    plt.show()

  11. 嗨,克里斯蒂安,

    非常感谢你的博客!

    由于没有关于s2polygon对象的文档,我想知道需要提供什么样的论据。
    你能给我一个指针吗?
    谢谢你!

  12. 所bepaly亚洲有发布的东西都很有意义。然而,考虑这个问题,假设你写了一个很棒的标题?
    我不是说你的信息不好,但是如果你添加了一个吸引人们注意的标题呢?我是说谷歌
    S2,球体上的几何体,Cells和Hilbert曲线Terra Incogbepaly亚洲nita有点香草味。

    你可以查看雅虎的主页并注意他们是如何创建新闻的。
    标题让人们点击。您可以添加视频或相关图片,或者
    二是让人们对你要说的话感兴趣。

    只是我的意见,这会让你的网站更有活力。

  13. 克里斯……我不知道为什么我花了这么长时间才找到这篇精彩的文章,但这对我来说是一个巨大的好处!!!

    非常感谢!bepaly亚洲

  14. 谢谢你的吉祥的评论。事实上,它曾经是一种享受。

    你那复杂的眼神更使我感到愉快!由
    的方式,我们怎么联系?

  15. 非常感谢这个神奇的教程,基督徒!

    不过,我有很多疑问。如果你能帮忙就太好了。.

    在这个例子中,您已经从两个点创建了一个矩形区域。这是怎么回事?
    我想从一千个点(拉特隆)创建一个区域。然后在细胞中分裂它的区域。我该怎么做?

    我可以看到在C++库中有一个用于ReCt的方法加点,但在python中没有。
    你能帮忙吗?谢谢!

  16. 嗨,克里斯蒂安,
    非常感谢这篇文章,非常全面!我唯一的问题是如何选择合适的面(立方体的6个面)来投射给定的(lat,lng)区域。
    最好的,
    –ZP公司

  17. 篇不错的基督徒。
    如何将此库与Android Studio一起使用?

  18. 与教科书相bepaly亚洲比,在网上找到任何东西都是很容易的,当我在这个网页上找到这篇文章时。

  19. 嗨,克里斯蒂安,
    非常感谢你的精彩博客。真是太棒了。

    我有一个用例,我有用户的地理位置和城市中的一个点,我需要找到这个点半径5公里以内的用户。我可以用s2单元ID表示用户和点的位置。但我该如何查询呢?
    是否需要将内存中所有用户的数据加载到应用程序中,然后使用库查找目标用户。或者我可以将所有用户的数据存储在数据库中,并对S2单元id进行查询,因为S2单元id使用Hilbert曲线,所以我们知道,S2单元id越近的地理点彼此之间的距离越近。
    或者还有我没有想到的其他方法吗
    提前感谢您的帮助

  20. 嗨,克里斯蒂安,
    非常感谢你的精彩博客。真是太棒了。

    我有一个用例,我有用户的地理位置和城市中的一个点,我需要找到这个点半径5公里以内的用户。我可以用s2单元ID表示用户和点的位置。但我该如何查询呢?
    是否需要将内存中所有用户的数据加载到应用程序中,然后使用库查找目标用户。或者我可以将所有用户的数据存储在数据库中,并对S2单元id进行查询,因为S2单元id使用Hilbert曲线,所以我们知道,S2单元id越近的地理点彼此之间的距离越近。
    或者还有我没有想到的其他方法吗
    提前感谢您的帮助

  21. 嗨,克里斯蒂安,感谢这篇文章。我只是想说,我必须使用“OSM”而不是“MapQuestOSM”,并在添加“geoms”时切换lat和lng的顺序

    所以不是:
    “`
    .degrees vertices.append(latlng.lat()(),
    .degrees latlng.lng () ()))
    “`
    我做了:
    “`
    顶点.append((latlng.lng().degrees(),
    latlng.lat().degrees()))
    “`

  22. 每bepaly亚洲件事都很有逻辑性。然而,这个怎么样?
    如果你补充一点信息呢?我不是说你的信息不好,但
    假设你添加了一个吸引人们注意的帖子标题?
    我是说谷歌的S2,球体上的几何学,Cells和Hilbert Cuurve Terra bepaly亚洲Incognita有点平原。你可以浏览雅虎的首页
    看看他们是如何创造头条新闻的
    人感兴趣。您不能添加视频或相关
    一两张图片让读者对你的作品感兴趣
    说。在我看来,这会让你的博客更有趣。

  23. 你好。我发现你不经常更新你的网站。我知道写文章既费时又无聊。
    但是你知道有一个工具可以让你创造
    新文章使用现有的内容(从文章目录或其他网站从您的利基)?
    它做得很好。bepaly亚洲新文章是独一无二的,并且通过了copyscape测试。
    你应该试试米夫托罗的工具

  24. 有人有javascript版本吗?我只是一个有reg映射的reg家伙我想实现这个。

留下你的评论

您的电子邮件地址将不会公布。

此网站使用Akismet来减少垃圾邮件。了解如何处理评论数据.