Uber H3 index 地图索引思考

H3 是 uber 设计的六边形空间索引，go 语言操作包是 h3-go，可以通过经纬度获取所在的 h3 六边形，每个经纬度对应的六边形都是确定的，他们唯一对应一个 h3index。

在业务开发中，我们可以通过 h3index 来对地理空间中的对象做聚合，本质还是将经纬度查询转换成了 h3index 的查询。但这种转换其实有很多种变通，比如说，获取两个经纬度连线上的所有 h3 网格，在继续扩展形式的话，获取某个围栏内的所有 h3 网格。好在官方库已经提供了现成的方法供我们使用。

大家都清楚经维度，经度范围 [-180, 180]，维度范围 [-90, 90]，它表示地理空间的一个固定点，类似二维坐标系中的 xy 坐标。而 h3index 表示一个面，它按不同的维度将地理空间划分成固定的面，已知一个 h3index 编码，它表示的地理空间是确定的。

应用场景

h3 网格可以圈定一个正六边形范围，并唯一对应一个 h3index，我们将 h3index 作为数据存储索引，以此来实现地理空间的范围查询。

比如，我们存储附近商圈的信息，使用商圈的中心点坐标计算 h3index，存储到底层数据库中。当计算用户附近的商圈时，使用用户所在的位置计算 h3index，然后扩圈兼容 h3 边界情况，就可以计算出附近的商圈。

不同维度的 h3index 面积不同，维度越大，网格的面积越小。详细的对应关系可以点击查看 Average area in m2，打趣的讲，我们的住房面积很少超过 12 级网格。

以12级网格为例，网格平均面积为 307 $m^2$，这差不多算是比较小的面积范围。如果是打车场景，我们可以假设该范围内的用户，上车点的位置大概率都是相同的。

更宽泛的说，在这个网格内的大部分用户，在某些场景特性下具有相似性，可以使用范围来做协同推荐。同一个网格内的用户，大概率会逛相同的超市或餐厅。

如上图，在之前的基础上，我向外扩一圈，会多出来了 6 个h3index。通过图中的 7 个 h3index，对虚线圆圈内的商圈做最热最优计算，具体的计算过程可以简化如下：

假设用户在最中心的网格，在虚线框定的范围内有 3 个商圈可以去，该给用户推荐那个商圈呢？可以结合扩圈加权打分的逻辑，和用户同网格的权重可以高一下，向外扩一圈的权重做递减，最终选出策略上的最优点。

网格本质

网格的本质是面，表示一个范围内的所有点。我们当前的位置需要使用经纬度表示，经纬度可以明确标志出地理空间的绝对位置，但点的粒度特别细，我稍微挪动一步，经纬度就变化了。

如果现在对用户位置做缓存，缓存的 key 是经度和维度的组合，那这个缓存的体量可能会非常大，而且是无限的。如果缓存用户当前所在的 h3 网格，缓存就会有明确的上限。而且，绝大多数情况下，选择最大维度的网格，其实也可以满足我们的业务需求。

通过用户所在的位置 lng、lat，计算出用户所在的网格 h3index，还可以计算出网格的中心点坐标。可以基于网格进行扩圈，也可以基于网格的中心点画圆。最终，框定出一个固定的范围。

另外，h3 网格是固定编码的，是不会随现实世界的变化而变化的。现实世界中往往有一套对应的 AOI、POI 编码，对应了一座固定的住宅社区、或者写字楼等，但这些 AOI、POI 可能会因为各种地理规划原因发生变化，但 h3 网格不会。

建立 h3 网格和 AOI、POI以及道路的关系，主要还是方便我们做数据召回。建立用户位置和 h3 网格的关系，不同的业务有不同的想法。

网格的核心便是点到面的转换，主观上，我们对面的感觉应该更清晰，业务上经常需要一个范围或者区域来划定边界。假设给你提供 3 个位置点信息，相互连线可以绘制一个三角形，但如果业务上需要绘制的是一个大于3的多边形，你有什么好的处理思路吗？希望可以看到大家的留言

网格维度

不同维度的 h3index 表示不同大小的地理空间，我们可以查看官方介绍来了解详细信息。h3 总共有 16 个维度，维度越高，表示的面积越小，也越精确。下图的对应关系表示，使用不同维度的 h3index 表示地球所需要的数量。

拿维度 15 来说，已经是 h3 能表示的最小单元了。用这个维度来表示地球的话，需要存储 569,707,381,193,162 个 h3index 索引。第3、4列表示其中六边形（hexagon）和五边形（pentagon）的组成个数。为什么组成中还会有五变形呢？ 我们应该从设计的根源上去探索解释。

使用 h3-go 包，我们通过指定经纬度、以及维度就可以获取到对应的 h3index 索引。还可以通过 h3index 获取所在六边形的组成点。下面代码中 FromGeo 就返回了对应的 h3index 。

geo := GeoCoord{
	Latitude:  37.775938728915946,
	Longitude: -122.41795063018799,
}
resolution := 9
fmt.Printf("%#x\n", FromGeo(geo, resolution))

相比较 h3index 能表示的最小单元格个数而言，我们最关心的还是单元格所能表示的面积大小，不同维度的 h3index 究竟能表示多大的范围。

如何计算网格

以计算 8 级网格为例，相关的函数FromGeo和ToGeoBoundary：

现在在地图上可视化查看 8 级网格的大小，主要就是确定8级网格的各个顶点坐标。对应的 sdk 都提供了对应的获取方法，我们选取一个经纬度，绘制这个经纬度所在的 8 级网格。我先在地图上拾取天台公园的坐标：116.410732,39.881836，根据这个坐标来查看对应的8级 h3 网格。

FromGeo 参数1指定位置信息，参数2设置维度为 8 级网格，计算得出对应的 h3index，然后调用 ToGeoBoundary 计算的出网格的各个顶点坐标（形成六边形围栏）。我们通过地图坐标的连线工具可以查看具体的边界。

geo := h3.GeoCoord{
		Longitude: 116.410732,
		Latitude:  39.881836,
	}
h3Index := h3.FromGeo(geo, 8)
boundary := h3.ToGeoBoundary(h3Index)

这样看 8 级网格还是少点意思，没有距离信息，啥都不是。但要在正六边形上标注出距离，我估计了一下，要花费不少时间，还是乖乖地当一个“调包侠”吧。用现成的API集成工具，查出两点间的距离，上下图对比的来看一下，也很容易看出对应标注的距离。

正六边形边的长度是 491米，两条平行边间的距离是814米，正对着两个点的距离是873米。现在我们去看官方对于 8 级网格平均边长的描述：0.531414010 km，看来这个正六边的边长小于平均值，可以点击 restable 查看详细信息。

日常工作中，最常用的还是网格的边长和面积，它们均能反应出我们想要覆盖的地理范围。13级网格的平均面积只有 43$m^2$，12级网格的平均面积是 307 $m^2$。在北上广深，大部分居住房子的面积都不会超过 12 级网格。

网格扩圈（KRing）

为什么要扩圈呢？假如你基于当前位置计算所在的 h3 8级网格，你虽然只能获取到一个 h3index，但你可能正好在网格的正中心，也可能恰好在网格的边缘。如果你想计算附近400米范围内的商店，如果你正好在网格的正中心，当前网格就能覆盖你的需求，但如果你在当前 h3 网格的边缘，那400米的范围会落到临近的网格内。所以，我们一般会选择扩圈，根据想要覆盖的范围，来决策扩几圈。

函数 KRing 方法用来计算 h3 扩圈，扩1圈就是当前当前h3网格按照同维度向外扩一圈，它会返回扩圈范围内的所有h3index。在上述代码的基础上，我们展示针对 h3Index 扩1圈和扩2圈的结果。

	indexs := h3.KRing(h3Index,1)
	fmt.Println(indexs)

	indexs := h3.KRing(h3Index, 2)
	fmt.Println(indexs)

为了让代码展示的更加直观，我专门处理一下返回的 indexs，能够在地图上勾勒出一个个网格图形。上面的代码 Println 也能展示出来，但展示的效果很不直观。我们还是配合地图多边形绘制工具，将扩圈后的网格展示出来。

func drawLine(indexes []h3.H3Index) string {
	draw := make([]string, 0)

	for i := range indexes {
		boundary := h3.ToGeoBoundary(indexes[i])
		points := make([]string, 0)
		for j := range boundary {
			points = append(points, fmt.Sprintf("[%f,%f]", boundary[j].Longitude, boundary[j].Latitude))
		}

		jsStr := `
		var path_%d = [
			%s
		]

		var polygon_%d = new AMap.Polygon({
		path: path_%d,
		strokeColor: "#FF33FF",
		strokeWeight: 6,
		strokeOpacity: 0.2,
		fillOpacity: 0.4,
		fillColor: '#1791fc',
		zIndex: 50,
		})

		map.add(polygon_%d)`

		draw = append(draw, fmt.Sprintf(jsStr, i, strings.Join(points, ","), i, i, i))
	}
	return  strings.Join(draw, "\n")
}

左图是基于中心点 h3index 扩1圈后的效果，右图是扩 2 圈后的效果。扩1圈返回 7 个六边形，扩 2 圈返回 19 个六边形。如果要解决边界边界问题，只需要向外扩一圈就可以了。

继续延伸的话，如果继续向外扩一圈，会覆盖多少个网格呢？我们通过上面提到的扩圈函数可以计算一下试试。还是以天坛公园为例，来看看扩圈后所覆盖的网格数量：

扩圈数	覆盖的总网格数
0	1
1	7
2	19
3	37
4	61
5	91
6	127
7	169

只要扩圈没有超出该维度网格的上限，理论上得有一个公式来计算出扩圈后能覆盖的网格数。话又说回来，我们一般最多也就扩2圈，有没有这个公式意义不大。

连线经过的网格 （Line）

Line函数计算两个经纬度连线所经过的所有H3索引，包含起点和终点所在的H3索引。生成的 H3 索引网格序列代表了从起点到终点六边形网格中心的近似“直线”路径。这个序列是有序的，按照从起点到终点的顺序排列。每个索引代表路径上连续的一个 H3 六边形网格。

逻辑上使用“线性插值”的技术，直白的将就是两个点连线，然后在这个线段上找点，计算这些点所在的网格，最后就是两个经纬度连线上的所有网格。特别需要注意的是：函数的入参是两个网格索引，而不是经纬度

// Line returns the line of h3 indexes connecting two indexes
func Line(start, end H3Index) []H3Index {
	n := C.h3LineSize(start, end)
	cout := make([]C.H3Index, n)
	C.h3Line(start, end, &cout[0])
	return h3SliceFromC(cout)
}

这个方法的好处还是转换，将线段转换为网格序列。线段在现实世界中可以表示路线，最终就形成了路线到网格的转换。

H3 网格设计

将真实世界的 3 维地理空间映射到 2 维的网格体系，是每个地图索引都需要面对的问题。H3 通过球心投影的方式将地球抽象成一个空间正20 面体，但地球本身也不是极致的圆形，抽象肯定会牺牲一部分准确性。

正20面体的每个面都是一个三角形平面，抽象到地图上的话，每个平面应该都是二维的平面。然后就变成了一个圆形的正20面体，就正好可以接受地图的投影。而且，如果把正 20 面体的每个面进行编码，我们可以通过选择正 20 面体，将合适的面对应到合适的地理区域上。

正 20 面体展开到二维平面，就是地球的投影，下面就是其中的一种展开方式。因为有陆地和海洋的区别，而我们要分析的地理数据基本都是在陆地上的。所以，正 20 面体通过合理的转动，可以将它的顶点都落在海洋的区域。为什么要特别处理这些顶点数据呢？

H3网格索引没有直接利用正 20 面体，而是在每个面的基础上继续做了网格拆分。为什么最终选择使用正六边形作为拆分的单元呢，还有没有别的多边形可以选择呢？其实还可以有三角形和正方形，GeoHash 采用的就是正方面的网格，那么正方形网格和正六边形有什么区别吗？

假设一个 n 面体，多面体的内角和计算公式为 (n-2)*180，假设 m 个正多边形在任意顶点出拼接形成一个 360 度才能形成网格，最终计算的表达式：(n-2)*180/n*m=360。如果 n 表示三角形，就需要有 6 个面拼接；如果 n 是正方形，就需要 4 个面拼接，如果 n 是六边形，就需要 3 个面拼接。

如果按照三角形或者四边形的拆分方式，某一个网格单元和临边的网格单元的中心距离是不相同的，如上图。但采用正六边形，网格单元和临边的其它网格单元是相同的。那么，在网格体系中，我们如何获取某个网格单元周边的网格呢？

H3 网格独立于现实中的区域，真实世界的区域边界往往会发生一些变化，区域的边界形态各异，经常还会发生变动调整，如果保持现实地理的区域跟踪，保证区域边界数据动态更新难度可想而知。试想一下，如果我们根据北京的每个区、每个街道进行数据统计，复杂性是非常大的。

而 H3 是在地球基础上的网格系统，它不用关心现实的街道、区域边界变化，无论现实中的地域如何变化，H3 都不会发生变化，也就是将地图上的某一个区域固定化了。如果要统计某个区域的维度信息，只需要先获取该区域的所有 H3 网格，汇总所有网格的信息就可以了。


在 0 级索引下，每个面包含了 10 个单元，如上图，其中边界线上的一些单元会被多个平面共同包含，这个面总共包含了 5.5 个六边形和 3 个三角形。上文提到的，正20面体所包含的 110 个六边形和 12 个五边形就是这么计算的出来的。正20面体包含有 20 个面、30条边、12 个顶点，其中的 12 个顶点就正好对应了这 12 个五边形。

网格关系

上图形象的描述了，2个等级 H3 网格的大小关系，高级别的 H3 网格大概需要 7 个才接近一个低级别的网格，两个级别之间不是完全的包含关系。

编码模式

h3index 提供了3种编码方式，最常用的属 cell mode 模式，就是我们上面示例所用到的，每个 h3index 表示一个正六边形网格。结合 KRing 的方法，我们可以在这个网格的基础上进行扩圈，返回的 h3index 还是 cell mode 模式

关于 cell model index ，官方文档中对编码做了详细解释。表格中每一行表示 16 bit，每一列表示将对应 16 bit 当做一组。0x 表示 16进制，0x30 正好等于 64。

我们以 85283473fffffff 为例，拆开成 64 bit 来看看。 首先，将 85283473fffffff 这个 16进制的整数转换成2进制的整数，我已经按每8位一组做了分隔：

0001000 01010010 10000011 01000111 00111111 11111111 11111111 11111111

从高位到低位依次看：

• 第一个bit （0）是保留位，被设置为0
• 4 个 bit 位（0001）表示 h3index 单元格的编码模式，cell mode 的编码为1，Directed edge mode 编码为2，Vertex mode 为4
• 3 个 bit （000）被保留，全部为 0。到这里，第一组编码就结束了
• 4 个 bit （0101）表示网格的维度，表示 0-15这16中情况，这个网格是 5 级网格。因为 4 个bit 能表示的最大整数是 15，所以，4个 bit 满打满算刚刚好。
• 7个 bit（0010 100）表示基础网格 0-121。0级别的网格总共有 122 个，现在位于 20 的基础网格
• 接下来还剩余 45 个bit，然后每3个bit为以小组，标识在不同维度下的坐标位置，其中未使用到的用 7 表示。所以，有效的部分为（00011 01000111 00），换算下来就是（0,6,4,3,4），对应的维度分别是（1,2,3,4,5），具体的坐标方位如下图所示

关于后续的模式的编码，我们就不一一深入解释了。从 h3index 的编码中也能给我们一些指导，首先确保编码唯一，其次确保编码有意义。

我们经常使用的打车场景，会给我们推荐一些上车点，其实可以采用这样的模式来对点进行编码。比如，点包含的一些基础属性，经度、纬度、道路、热度等信息，现在给这个点新加一个唯一ID来标识，怎么设计这个点ID？无脑的简单设计就是经度、纬度的简单拼接。但其实设计的可以更有技术些。可以将对应的h3index网格信息加入到ID中，将经度和纬度的浮点数换算成更小的整数，在加入一些点的类型、热度等信息。

第二种编码模式是 Directed edge mode，directed 含义是有方向的，姑且称呼为有向边模式，是对六边形的边进行编码。h3 package 中提供了获取单向边的方法：ToUnidirectionalEdges ，我们通过方法调用，查看一下 85283473fffffff 的单向边。

h3Index := h3.FromString("85283473fffffff")
edges := h3.ToUnidirectionalEdges(h3Index)
for e := range edges {
	fmt.Println(h3.ToString(edges[e]))
}

// output:
115283473fffffff
125283473fffffff
135283473fffffff
145283473fffffff
155283473fffffff
165283473fffffff

程序会输出 6 条到附近 h3index 的单向边，通过单向边，我们可以获取单向边指向的 h3index，我们以 115283473fffffff 为例，查看它所指向的目的 h3index。

edgeIndex := h3.FromString("115283473fffffff")
h3index := h3.DestinationFromUnidirectionalEdge(edgeIndex)
fmt.Println(h3.ToString(h3index))

// output:
85283477fffffff

我们在地图上绘制这两个 h3index 的位置，通过这个边可以查到的相邻的 cell mode 模式下的 h3index。我看对应的 edge 相关的方法还有FromUnidirectionalEdge，它会返回两个h3index，理论上就是下面这两个了。

edgeIndex := h3.FromString("115283473fffffff")
origin,dest := h3.FromUnidirectionalEdge(edgeIndex)
fmt.Println(h3.ToString(origin), "->", h3.ToString(dest))

// output:
85283473fffffff -> 85283477fffffff

第三种属于对六边形的点进行编码，Vertex mode，Vertex 表示顶点，姑且称呼为顶点编码。比较遗憾的是，h3 包没有提供顶点相关的方法。下面是顶点的编码，按道理来说，通过顶点可以获取到 3 个 cell mode。

Mongo

mongo 自身集成了空间索引的能力，说到底，地理空间的数据都需要被存储和被索引，如何存储位置数据以及如何查询位置数据，是每个业务都要面临的难题。我们专门去看看 mongo 提供了哪些空间能力。

GeoJSON 对象

地理空间中经度和维度可以唯一确定一个位置，这属于位置点信息。市面上很多基于用户位置提供的服务，都是经纬度的一些聚合运算。GeoJSON 是 Mongo 中存储空间数据的对象，这种对象类型的存储模式，我一直都觉得应该是一种趋势的。

最基本的对象，也是地理空间里最基本的概念：点、线、面。参考一下它面的结构，一般来说，面是线段收尾拼接组成的，关键是看看 Mongo 中如何表示面这个结构。下面就是多边形的表示，比较有特点的是首尾点是重复的。

{
  type: "Polygon",
  coordinates: [ [ [ 0 , 0 ] , [ 3 , 6 ] , [ 6 , 1 ] , [ 0 , 0  ] ] ]
}

想一想，给我们一堆点，我们如何把这堆点两两连接起来，最后形成一个面呢？关键就是看最后一个点是否要明确的声明出来。从这个面的对象也可以看出，Mongo 它支持不规则的面空间，这个多边形完全由用户自己定义。

Geospatial

Mongo 提供了哪些便捷的空间查询呢，既然有了点、线、面的空间对象定义，对应的操作肯定也少不了。

geoIntersects

这个方法会应用在哪些业务场景下，或者说，我们在哪些场景中才需要用到这个方法，是我们关注的重点。函数名称是空间相交，但其实就是圈定一个地理范围，在这个范围内进行其他位置数据查找。

参考文章：

1. Uber’s Hexagonal Hierarchical Spatial Index
2. h3geo文档
3. mongo geojson