Google Earth Engine（GEE）——使用 GeoPandas 和 Uber 的 H3 空间索引进行快速多边形点分析

空间索引方法有助于加速空间查询。大多数 GIS 软件和数据库都提供了一种机制来计算和使用数据图层的空间索引。QGIS 和 PostGIS 使用基于 R-Tree 数据结构的空间索引方案 - 它使用几何边界框创建分层树。这是非常有效的，并在某些类型的空间查询中产生了很大的加速。查看我的高级 QGIS 课程的空间索引部分，我将展示如何在 QGIS 中使用基于 R 树的空间索引。

如果您使用 Python 进行地理处理，GeoPandas 库还提供了使用 .sidex 属性的基于 R-Tree 的空间索引的易于使用的实现。赫尔辛基大学的 AutoGIS 课程有一个很好的例子，将空间索引与 geopandas 一起使用。

在这篇文章中，我想谈谈另一个名为H3 的空间索引系统。

这个开源索引系统由 Uber 创建，使用六边形网格单元。该系统类似于另一个名为S2 的基于单元格的索引系统——它是在谷歌开发的。这两个系统都提供了一种将地球上的坐标转换cell id为以特定分辨率映射到六边形或矩形网格单元的方法。这些单元格 id 具有独特的属性，例如附近的单元格具有相似的 id，您可以通过截断它们的长度来找到父单元格。这些属性使得诸如聚合数据、查找附近对象、测量距离之类的操作非常快速。

在这篇文章中，我将向你展示如何创建使用点密度图geopandas和h3-py库在Python。

国家地理空间情报局的海事安全信息门户以反航运活动消息的形式提供所有海盗事件的形状文件。

数据以 zip 文件形式提供ASAM_shp.zip。实际的数据层是一个ASAM_events.shp位于文件夹内的shapefile ASAM_data_download。该数据集包含全球 8000 多个已记录盗版事件的点位置。这是原始点图层在 QGIS 中的可视化效果。

我们将通过在 H3 提供的六边形网格上聚合事件点来创建密度图。我们从导入库开始。

import geopandas as gpd
from h3 import h3

GeoPandas 允许直接从 zip 文件中读取数据层。

zipfile = 'zip://ASAM_shp.zip!ASAM_data_download/ASAM_events.shp'
incidents = gpd.read_file(zipfile)

H3 库有 15 级分辨率。 此表 显示了每个级别的详细信息。我们选择级别 3，这导致网格大小约为 100 公里。该函数lat_lng_to_h3将位置的坐标转换为所选级别的 H3 id。我们h3为级别 3 的点添加一个名为H3 网格 ID的列。

h3_level = 3
 
def lat_lng_to_h3(row):
    return h3.geo_to_h3(
      row.geometry.y, row.geometry.x, h3_level)
 
incidents['h3'] = incidents.apply(lat_lng_to_h3, axis=1)

现在是有趣的部分。由于落在网格单元中的所有点都具有相同的 id，我们可以简单地聚合具有相同网格 id 的所有行，以找到落在网格多边形中的所有点。因此，通过使用基于网格的索引系统 - 复杂的空间“多边形点”操作变成了对表的简单聚合。我们groupby在h3列上使用 Panda 的函数，并count在输出中添加一个新列，其中包含每个 H3 id 的行数。

counts = incidents.groupby(['h3']).h3.agg('count').to_frame('count').reset_index()

我们现在知道每个网格单元的盗版事件数量。要将结果可视化或将其导出到 GIS，我们需要将 H3 单元 ID 转换为几何图形。该 h3_to_geo_boundary 函数采用 H3 键并返回形成六边形单元格的坐标列表。由于 GeoPandas 使用 shapely 库来构建几何，我们将坐标列表转换为一个匀称的 Polygon 对象。请注意h3_to_geo_boundary 我们设置的函数 的可选第二个参数， 与默认(lat,lon)相比，True 它返回(x,y)顺序中 的坐标

from shapely.geometry import Polygon
 
def add_geometry(row):
    points = h3.h3_to_geo_boundary(
      row['h3'], True)
    return Polygon(points)
 
counts['geometry'] = counts.apply(add_geometry, axis=1)

我们将数据框转换为带有 CRS EPSG:4326（WGS84 纬度/经度）的 GeoDataframe 并将其写入地理包。

gdf = gpd.GeoDataFrame(counts, crs='EPSG:4326')
output_filename = 'gridcounts.gpkg'
gdf.to_file(driver='GPKG', filename=output_filename)

gridcounts.gpkg 可以在 QGIS 中打开和可视化。这是显示生成的 hexbin 地图的图层，其中显示了世界各地的盗版热点。

从读取输入到创建聚合网格层的整个过程只需 2 秒多一点。将其与使用空间索引的 QGIS 模型进行比较，该模型至少需要 5 倍。H3 特别适合这种空间聚合并且速度非常快。

这篇文章中使用的代码和数据集可以在我的Github 存储库中找到。您还可以在 Binder 中实时运行 Jupyter Notebook 。