数据, 术→技巧

空间索引之Uber H3

钱魏Way · · 40 次浏览

网格系统(Grid System)对于分析海量空间数据集,将地球空间划分为可识别的网格单元(cell)至关重要。H3是由Uber开源的一个六边形分层索引网格系统,也是最近几年实现数据聚合的主要趋势,在H3出现之前大部分情况采用的是Geohash算法墨卡托投影,还有一些其他投影技术,比如Google S2地理索引。

H3网格系统在Uber的主要应用为:

  • 优化乘车价格和调度(动态定价)
  • 地图空间数据的可视化和挖掘
  • 用于整个市场的分析和优化

GEOHASH存在的问题

GeoHash的原理就是按照区域,把经纬度进行编码,按照不同的网格精度,变成不同位数的编码,在同一区域中的编码相同。

问题1:不同精度下网格的形状不一(长方形和矩形)且精度的变化幅度时小时大

问题2:在不同纬度的地区会出现地理单元单位面积差异较大的情况

在国际化业务中会表现特别明显,比如北京和新加坡的 geohash 对应面积有将近 30% 的差异。这导致业务指标和模型输入的特征存在一定的分布倾斜和偏差。

问题3:存在8邻域到中心网格的距离不相等的问题

Uber H3的设计思路

六边形网格索引

H3是一种基于网格的空间索引,但跟普通的矩形网格索引不同的是,他的每一个网格都是正六边形。为啥要选正六边形呢,因为在基于网格的空间索引中,使用的多边形的边数越多,则一个网格越近似圆形,做缓冲区查询、kNN查询什么的也就越方便。而做网格索引又要求空间能够被网格铺满,不能有缝隙。根据初中数学知识,我们知道一个多边形的内角和公式为:

$$\theta = (x – 2) * 180^\circ$$

其中,x为多边形的边数,$\theta$为多边形的内角和。则一个正多边形的每个角的角度为 $\frac{\theta}{x}=\frac{(x – 2) * 180^\circ}{x}$,而如果需要多边形能够铺满空间,则在多边形的顶点相交处,设每个顶点有y个多边形相交,需要满足以下等式:

$$\frac{360^\circ}{y} = \frac{(x – 2) * 180^\circ}{x}$$

以上等式的求解过程我不再赘述,这个等式只有三组整数解:x=3,y=6;x=4,y=4;x=6,y=3。

因此,能够做网格空间索引的形状只有三角形、矩形和六边形,而六边形因为边数最多,最接近圆,所以理论上来说在某些场景下是最优的选择。

Uber的技术博客上也给出了使用六边形做空间索引的优越性:

上图展示了三种网格到其相邻网格的距离,六边形网格与周围网格的距离有且仅有一个,而四边形存在两类距离,三角形有三类距离。所以,H3就使用六边形作为网格索引的基本单元,实现空间索引。

无变形的投影

全球网格系统通常至少需要两件事:地图投影和放置在地图上的网格。从地球上的三维位置到地图上的二维点需要地图投影 然后将网格覆盖在地图上,形成一个全球网格系统。在GIS领域中,对空间填充曲线熟悉的同学应该知道,不论是GeoHash, Z2或者Hilbert,虽然看起来都是将空间按照经纬度分割成了一个个大小相等的网格,但实际上这些网格的实际面积并不相等。对于靠近极地的网格,虽然经纬度的间隔没变,但由于地球的曲率,这些网格的实际面积远小于靠近赤道的网格。

这种实际面积不相等的网格索引可能会造成一个问题,那就是由于网格大小不一致导致网格内数据量不一致,造成热网格和冷网格。于是,H3干脆摒弃传统的地图投影,直接在地球上铺满六边形。

想要实现六边形的网格,最大的困难时如何变化。H3 的前身其实是 DDGS(Discrete global grid systems) 中的 ISEA3H,其原理是把无限的不规则但体积相等的六棱柱从二十面体中心延伸,这样任何半径的球体都会穿过棱镜形成相等的面积cell,基于该标准使得每一个地理单元的面积大小就可以保证几乎相同。然而原生的 ISEA3H 方案在任意级别中都存在12个五边形,H3 的主要改进是通过坐标系的调整将其中的五边形都转移到水域上,这样就不影响大多数业务的开展。

H3实现的方法是:将地球当作一个二十面体,这个二十面体的每一个面都是球面三角形,有12个顶点,称为球形二十面体(spherical icosahedron),在这个球形十二面体的每个面上都有相同排列方式的六边形。由于这个球形二十面体的12个顶点每一个都在地球上的水里,可以保证对于每个面做处理时不会遇到边界的edge case,因为Uber还没有轮船快艇业务,只需要保证在陆地上H3好使就行。这个球形二十面体的示意图见图。

在第0层,这个球形二十面体的每一个面长这样:

可以看到,在一个球面三角形的顶角处有个小三角形,这个小三角形就是H3的一种edge case:在二十面体的顶点处,有五个面交于这个顶点,每个面在这个顶点处都有一个小三角形,所以这些小三角形会形成一个五边形。也就是说,H3并不能保证每个空间单元都是六边形,在一些地方还是会存在五边形,但是这样做也不会造成很大影响,因为根据球形二十面体每个顶点都在水里的特性,这种五边形只会出现在水域周围,不会对Uber的打车和外卖业务造成很大的影响。根据这样的索引特性,H3规定在索引的第0层,每个面和上图一样,每个面上有5.5个六边形和3/5个五边形,即第0层一共有110个六边形和12个五边形。H3将这110个六边形称为基网格(base cells)。H3最高可以到15层,也就是说H3有16个层级的空间索引粒度,在粒度最细的第15层中,平均每个网格的大小为0.9平方米,平均边长为0.509713米。在往下一层划分时,每个父网格对应7个子网格,父子网格之间的对应关系是这样的:

 

可以看到,父子网格之间并没有严密的对齐,父网格和其所对应的7个子网格之间会有一点差异,这种差异也导致了H3并不能表现出很好的层级关系。

详细的精度数据如下:

网格编码

对六边形索引,H3使用了一种IJK坐标系来确定六边形的位置。IJK坐标系长这样:

当然,H3不会简单地使用这种坐标系来对所有六边形进行编码,因为在每一个层级六边形的排列方式都不太一样。总的来说,在所有层级中,六边形的排列方式只有两种类型,称为Class II和Class III。在这两种类型的排列方式中,IJK坐标系的三个轴的方向不太一样。对于这种在一个二十面体的面上,根据不同的六边形排列方式使用不同方向坐标轴的坐标系,H3称之为FaceIJK坐标系:

那么对于一个面上所有层级的所有六边形,都使用FaceIJK坐标系编码,再加上面的唯一标识符可不可以?答案是可以,但没必要,因为如果这样做,H3编码会更加表示不了层级之间的关系。H3为了突出层级之间的关联性,使用了一种方法:每个六边形都包含其父六边形的坐标。这样只需要规定好每个网格的子网格坐标的计算方法,对于子网格,只需在父网格的坐标后面追加子网格的坐标即可。这样一来,只需关注一个网格的7个子网格如何计算坐标,所有层级的每个网格的坐标都可以递归得到,这7个网格坐标的计算方法见图,为了表述方便,称之为IJK七网格坐标系。那么还有12个五边形怎么办?那12个五边形随着层级划分仍然是12个,层级越高这些五边形越小,并且都在海里,影响不大,直接不管了。

同时,H3还专门设置了一种叫做网格有向边(directed edges of grid cells)的东西(如图),其目的是为了网格能够快速地找到其某个邻网格。这种网格有向边分为单向网格有向边和双向网格有向边,比如说,对于两个相邻网格A和B,其相交的边是E,如果E是一条由A->B的单向边,那么只能通过A快速找到邻居B,B不能快速找到邻居A;而如果E是一条双向边,那么A可以快速地找到B,B也可以快速地找到A。

H3的索引值最多占63个比特位,可以用一个长整型表示,其结构如下所示:

  • 0-3位:索引模式。0表示无效,1表示普通的网格,2表示单向边,3表示双向边
  • 4-6位:如果索引模式是0或1,则没用;如果是2或3,表示这条边是这个六边形的哪条边,取值范围1-6
  • 7-10位:表示层级,取值范围0-15
  • 11-17位:表示这个网格属于哪个基网格,取值范围0-121
  • 18-21位:表示这个网格的第一代祖宗网格在IJK七网格坐标系下的坐标
  • n-n+2(n<=61)位:同上,表示这个网格的第i代祖宗网格在IJK七网格坐标系下的坐标

可以看到,其实对于层级数比较小的网格,后面很多位是不需要的,如果全置为0则会浪费很多空间,所以H3索引值是可以压缩的。

Uber H3使用示例

Uber平台每天都会发生数以百万计数的事件,包括但不限于:打车、拼车,外卖。 每个事件都发生在特定的位置,如,乘客在自己家中发起乘车订单,而驾驶员则在几英里外的汽车中接受该订单。这些事件使Uber能够更好地了解和优化我们的服务。 例如,可以告诉我们在城市的某个地方需求大于供应,然后就可以做出相应的价格调整,或者通知平台的某个Uber驾驶员附近有两个乘车订单。

要从Uber业务数据中获取信息和知识,需要分析整个城市的数据。 由于城市在地理位置上是多种多样的,因此需要以精细的粒度进行分析。以精细的粒度(事件发生的确切位置)进行分析是非常困难且费力的。对区域(例如城市中的社区)进行分析更加切实可行。

上面的地图描绘了使用H3进行分割的过程。我们使用网格系统将事件存储到六边形区域(即网格)中。我们通过测量所服务的每个城市中六边形的供求来计算高峰价格。 这些六边形构成了我们分析Uber业务的基础。

六边形是一个重要的选择,六边形可以最大程度地减少用户穿越城市时引入的量化误差。还可以轻松地获取近似半径。下图为纽约市(曼哈顿)的按邮政编码划分的地图:

这些区域具有不规则的形状和大小,对分析没有帮助,并且区域划分随时发生改变,而这些改变与我们使用它们的目的完全无关的。 Uber运营团队还可以根据他们对城市的了解来划定区域,但是随着城市的变化,这些区域需要经常更新,并且经常随意定义区域的边缘。

网格系统在Uber运营所在的城市中应该具有可比的形状和大小,并且不受任意更改。 尽管网格系统无法与城市中的街道和社区保持一致,但可以通过对网格单元进行聚合来将其有效地表示社区。 可以使用目标函数完成聚合,从而生成对分析更有用的形状。 确定聚类的成员与设置查找操作也应该是很高效的。随机生成的六边形簇覆盖旧金山市:

Uber H3的安装与使用

H3的安装

H3在使用pip install h3时,容易发生如下错误:

Looking in indexes: http://mirrors.aliyun.com/pypi/simple/
Collecting h3
  Downloading http://mirrors.aliyun.com/pypi/packages/98/2a/704b2db80465ce3a412d6c94f1c1658cffbbb4e76c5c03bb133f588cc131/h3-3.4.3.tar.gz
Building wheels for collected packages: h3
  Building wheel for h3 (setup.py) ... error
  ERROR: Command errored out with exit status 1:
   command: /opt/conda/bin/python -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-h2kzp6ux/h3/setup.py'"'"'; __file__='"'"'/tmp/pip-install-h2kzp6ux/h3/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' bdist_wheel -d /tmp/pip-wheel-kqgcncvz --python-tag cp37
       cwd: /tmp/pip-install-h2kzp6ux/h3/
  Complete output (47 lines):
  running bdist_wheel
  running build
  running build_py
  creating build
  creating build/lib.linux-x86_64-3.7
  creating build/lib.linux-x86_64-3.7/h3
  copying h3/__init__.py -> build/lib.linux-x86_64-3.7/h3
  copying h3/h3.py -> build/lib.linux-x86_64-3.7/h3
  running build_ext
  + VERSION=v3.4.2
  + IS_64BITS=True
  + '[' '' == v3.4.2 ']'
  + command -v cmake
  + echo 'cmake required but not found.'
  cmake required but not found.
  + exit 1
  Traceback (most recent call last):
    File "<string>", line 1, in <module>
    File "/tmp/pip-install-h2kzp6ux/h3/setup.py", line 65, in <module>
      distclass=BinaryDistribution)
    File "/opt/conda/lib/python3.7/site-packages/setuptools/__init__.py", line 145, in setup
      return distutils.core.setup(**attrs)
    File "/opt/conda/lib/python3.7/distutils/core.py", line 148, in setup
      dist.run_commands()
    File "/opt/conda/lib/python3.7/distutils/dist.py", line 966, in run_commands
      self.run_command(cmd)
    File "/opt/conda/lib/python3.7/distutils/dist.py", line 985, in run_command
      cmd_obj.run()
    File "/opt/conda/lib/python3.7/site-packages/wheel/bdist_wheel.py", line 192, in run
      self.run_command('build')
    File "/opt/conda/lib/python3.7/distutils/cmd.py", line 313, in run_command
      self.distribution.run_command(command)
    File "/opt/conda/lib/python3.7/distutils/dist.py", line 985, in run_command
      cmd_obj.run()
    File "/opt/conda/lib/python3.7/distutils/command/build.py", line 135, in run
      self.run_command(cmd_name)
    File "/opt/conda/lib/python3.7/distutils/cmd.py", line 313, in run_command
      self.distribution.run_command(command)
    File "/opt/conda/lib/python3.7/distutils/dist.py", line 985, in run_command
      cmd_obj.run()
    File "/tmp/pip-install-h2kzp6ux/h3/setup.py", line 25, in run
      install_h3(h3_version)
    File "/tmp/pip-install-h2kzp6ux/h3/setup.py", line 18, in install_h3
      subprocess.check_call('bash ./.install.sh {} {}'.format(h3_version, is_64bits), shell=True)
    File "/opt/conda/lib/python3.7/subprocess.py", line 363, in check_call
      raise CalledProcessError(retcode, cmd)
  subprocess.CalledProcessError: Command 'bash ./.install.sh v3.4.2 True' returned non-zero exit status 1.
  ----------------------------------------
  ERROR: Failed building wheel for h3
  Running setup.py clean for h3
Failed to build h3
Installing collected packages: h3
    Running setup.py install for h3 ... error
    ERROR: Command errored out with exit status 1:
     command: /opt/conda/bin/python -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-h2kzp6ux/h3/setup.py'"'"'; __file__='"'"'/tmp/pip-install-h2kzp6ux/h3/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /tmp/pip-record-ccmfladx/install-record.txt --single-version-externally-managed --compile
         cwd: /tmp/pip-install-h2kzp6ux/h3/
    Complete output (49 lines):
    running install
    running build
    running build_py
    creating build
    creating build/lib.linux-x86_64-3.7
    creating build/lib.linux-x86_64-3.7/h3
    copying h3/__init__.py -> build/lib.linux-x86_64-3.7/h3
    copying h3/h3.py -> build/lib.linux-x86_64-3.7/h3
    running build_ext
    + VERSION=v3.4.2
    + IS_64BITS=True
    + '[' '' == v3.4.2 ']'
    + command -v cmake
    + echo 'cmake required but not found.'
    cmake required but not found.
    + exit 1
    Traceback (most recent call last):
      File "<string>", line 1, in <module>
      File "/tmp/pip-install-h2kzp6ux/h3/setup.py", line 65, in <module>
        distclass=BinaryDistribution)
      File "/opt/conda/lib/python3.7/site-packages/setuptools/__init__.py", line 145, in setup
        return distutils.core.setup(**attrs)
      File "/opt/conda/lib/python3.7/distutils/core.py", line 148, in setup
        dist.run_commands()
      File "/opt/conda/lib/python3.7/distutils/dist.py", line 966, in run_commands
        self.run_command(cmd)
      File "/opt/conda/lib/python3.7/distutils/dist.py", line 985, in run_command
        cmd_obj.run()
      File "/opt/conda/lib/python3.7/site-packages/setuptools/command/install.py", line 61, in run
        return orig.install.run(self)
      File "/opt/conda/lib/python3.7/distutils/command/install.py", line 545, in run
        self.run_command('build')
      File "/opt/conda/lib/python3.7/distutils/cmd.py", line 313, in run_command
        self.distribution.run_command(command)
      File "/opt/conda/lib/python3.7/distutils/dist.py", line 985, in run_command
        cmd_obj.run()
      File "/opt/conda/lib/python3.7/distutils/command/build.py", line 135, in run
        self.run_command(cmd_name)
      File "/opt/conda/lib/python3.7/distutils/cmd.py", line 313, in run_command
        self.distribution.run_command(command)
      File "/opt/conda/lib/python3.7/distutils/dist.py", line 985, in run_command
        cmd_obj.run()
      File "/tmp/pip-install-h2kzp6ux/h3/setup.py", line 25, in run
        install_h3(h3_version)
      File "/tmp/pip-install-h2kzp6ux/h3/setup.py", line 18, in install_h3
        subprocess.check_call('bash ./.install.sh {} {}'.format(h3_version, is_64bits), shell=True)
      File "/opt/conda/lib/python3.7/subprocess.py", line 363, in check_call
        raise CalledProcessError(retcode, cmd)
    subprocess.CalledProcessError: Command 'bash ./.install.sh v3.4.2 True' returned non-zero exit status 1.
    ----------------------------------------
ERROR: Command errored out with exit status 1: /opt/conda/bin/python -u -c 'import sys, setuptools, tokenize; sys.argv[0] = '"'"'/tmp/pip-install-h2kzp6ux/h3/setup.py'"'"'; __file__='"'"'/tmp/pip-install-h2kzp6ux/h3/setup.py'"'"';f=getattr(tokenize, '"'"'open'"'"', open)(__file__);code=f.read().replace('"'"'\r\n'"'"', '"'"'\n'"'"');f.close();exec(compile(code, __file__, '"'"'exec'"'"'))' install --record /tmp/pip-record-ccmfladx/install-record.txt --single-version-externally-managed --compile Check the logs for full command output.

通过分析,发现编译环境缺少cmake导致的。

  • 解决方案1:如果安装了anaconda,则可以使用conda install h3-py -c conda-forge进行安装
  • 解决方案2:下载编译好的.whl文件进行安装h3-3.4.3-cp37-cp37m-win_amd64.whl.zip
  • 解决方案3:安装cmake,以Linux系统为例:
sudo apt-get install libssl-dev
wget https://github.com/Kitware/CMake/releases/download/v3.17.0/cmake-3.17.0.tar.gz
tar -zxvf cmake-3.17.0.tar.gz
cd cmake-3.17.0
./bootstrap
make
sudo make install
cmake --version

H3的使用

使用示例:

from h3 import h3

for i in range(0, 16):
    print(h3.geo_to_h3(lat=23.103804, lng=113.370424, res=i))  # 经纬度转化为H3索引
print(h3.h3_to_geo(h3_address="874118b24ffffff"))  # H3索引转化为经纬度(中心点)
print(h3.h3_to_geo_boundary(h3_address="874118b24ffffff"))  # H3索引转化为经纬度(六个顶点)
print(h3.k_ring(h3_address="874118b24ffffff", ring_size=1))  # 相邻网格(k代表拓展的层数)
print(h3.k_ring_distances(h3_address="874118b24ffffff", ring_size=2))  # 讲相邻网格按距离拆分成列表
print(h3.hex_ring(h3_address="874118b24ffffff", ring_size=1))  # 给出相邻的环
print(h3.hex_range_distances(h3_address="874118b24ffffff", ring_size=1))
print(h3.hex_ranges(h3_address="874118b24ffffff", ring_size=1))

可视化工具:

  • Folium
  • desk.gl (Uber’s large-scale WebGL-powered data visualization framework)
  • kepler.gl

Uber H3实战:英国交通事故点聚类

主要流程:

  • 将交通事故的经纬度信息转化为Uber H3
  • 方案一:对交通事故所在的经纬度进行聚类,获取非-1类别的数据,统计在聚集点的每个Uber H3的交通数据数量,并在地图上呈现
  • 方案二:对交通事故所在的H3网格进行汇总统计,针对网格所在区域按照阈值筛选后在地图上进行呈现

具体代码:

import numpy as np
import pandas as pd
import folium
from h3 import h3
from sklearn.cluster import DBSCAN

def create_map(clusters):
    # Create the map object
    map = folium.Map(zoom_start=12)

    # Convert the clusters dictionary items to polygons and add them to the map
    for cluster in clusters.values():
        points = cluster['geom']
        # points = [p[::-1] for p in points]
        tooltip = "{0} accidents".format(cluster['count'])
        polygon = folium.vector_layers.Polygon(locations=points, tooltip=tooltip,
                                               fill=True, 
                                               color='#ff0000', 
                                               fill_color='#ff0000', 
                                               fill_opacity=0.4, weight=3, opacity=0.4)
        polygon.add_to(map)

    # Determine the map bounding box
    max_lat = df.Latitude.max()
    min_lat = df.Latitude.min()
    max_lon = df.Longitude.max()
    min_lon = df.Longitude.min()
    
    # Fit the map to the bounds
    map.fit_bounds([[min_lat, min_lon], [max_lat, max_lon]])
    
    return map

column_types = {'Accident_Index': np.string_, 'LSOA_of_Accident_Location': np.string_}
uk_acc = pd.read_csv("UK_Accidents_2015_2016.csv", dtype=column_types)

#将经纬度转化成H3
h3_level = 11
def lat_lng_to_h3(row):
    return h3.geo_to_h3(row['Latitude'], row['Longitude'], h3_level)
uk_acc['h3'] = uk_acc.apply(lat_lng_to_h3, axis=1)

#使用DBSCAN进行聚类
uk_acc['rad_lng'] = np.radians(uk_acc['Longitude'].to_numpy())
uk_acc['rad_lat'] = np.radians(uk_acc['Latitude'].to_numpy())

eps_in_meters = 50.0
num_samples = 10
EARTH_R = 6370996.8 # 地球半径
uk_acc['cluster'] = DBSCAN(eps=eps_in_meters/EARTH_R, min_samples=num_samples, metric='haversine').fit_predict(uk_acc[['rad_lat', 'rad_lng']])

df = uk_acc[uk_acc.cluster != -1].copy()
clusters = dict()

for index, row in df.iterrows():
    key = row['h3']
    if key in clusters:
        clusters[key]['count'] += 1
    else:
        clusters[key] = {"count": 1,"geom": h3.h3_to_geo_boundary(h3_address=key)}

create_map(clusters)    

# 直接H3代码进行分组
h3_clusters = dict()

for index, row in uk_acc.iterrows():
    key = row['h3']
    if key in h3_clusters:
        h3_clusters[key]["count"] += 1
    else:
        h3_clusters[key] = {"count": 1, "geom": h3.h3_to_geo_boundary(h3_address=key)}

relevant_clusters = { key:value for (key,value) in h3_clusters.items() if value['count'] >= 10}
create_map(relevant_clusters)

参考链接:

发表评论

电子邮件地址不会被公开。 必填项已用*标注