「H3」はUberが開発している「空間インデックス」で、地球上を複数の解像度と六角形（まれに、五角形）で表現します。地理座標（経度・緯度）と解像度を指定すると、その地点を含む六角形セルを一意に表すH3インデックス（セルのID）が取得できます。
H3を利用すると、以下のようなことが簡単に行えます。
- 広域データを「六角形セル単位」で集計する
- 近隣セルの探索（隣接している六角形を簡単に取得できる）
- 解像度の変更（より粗い解像度のセルに変換したり、その逆を行う）
- ポリゴン領域に対して含まれるセルの一覧を取得する
（h3）
4分木のXYZタイルやなどでもよく見られるような正方形を選択しなかったのは以下のような理由があるようです。
- 六角形なら、全ての隣接するセルとの距離が等しい
- 正方形で空間充填するよりも、六角形で埋める方が誤差が少ない
- など
今回は、このH3インデックスを色々触ってみます。
Overture Mapsからデータをダウンロード
Overture Maps Foundationは、Amazon、Meta、Microsoft、TomTomなどの企業が中心となって2022年12月に立ち上げた地図関連のオープンデータプロジェクトです。
建物の形状データや道路のデータなど、多種多様な地理情報が含まれており、自由に利用・解析・可視化などを行うことができます。
この辺りにも色々書いています。
今回の記事では、このOverture Mapsが提供している建物データを使い、特定範囲に含まれる建物を取得してH3インデックスを付与し、可視化と解析を行います。
DuckDBでデータのダウンロード
まず、Overture Mapsが公開している建物データ(Parquet形式) を、DuckDB を使ってダウンロードし、扱いやすいGeoPackage形式に変換してみます。DuckDBは、Azure上のファイルやS3のファイルを直接クエリできるため、膨大な地図データから必要な範囲のみ取り出すことが出来ます。
詳しいことはこちらの記事を参考にしてください。
また、DuckDBを利用したOverture Mapsのダウンロード方法については、公式の記事を参照してください。
https://docs.overturemaps.org/getting-data/duckdb/?query=Places
Overture MapsはWebのビューワーが用意されているので、欲しい地物・属性が存在しているかどうか、などはある程度すぐに確認することができます！
データを取得したい範囲を確認したい場合は、geojson.ioが便利です。
データはAzureかS3のどちらかに配置されているものを選べますが、僕がクエリしたタイミングでは運悪く(なのかは分かりませんが)S3から取得することが出来なかったのでAzureの方を利用しています。
PythonのDuckDBやJupyter Notebookを利用していきますので、適宜環境構築して実行していってみてください。
import duckdb
# メモリ上でDuckDBを起動
con = duckdb.connect(database=':memory:')
# 必要なDuckDB拡張機能をインストール＆ロード
con.execute("INSTALL spatial; LOAD spatial;")
con.execute("INSTALL azure; LOAD azure;")
# Overture Mapsの最新リリースの建物データ
# 最新リリースの情報は定期的に確認する必要があります。
release = "2025-02-19.0"
data_path = f"azure://release/{release}/theme=buildings/type=building/*"
# 解析対象となるBBOX
bbox = {
"type": "Polygon",
"coordinates": [
[
[
139.56721442796606,
35.77552476978808
],
[
139.56721442796606,
35.5781480840387
],
[
139.9653733196339,
35.5781480840387
],
[
139.9653733196339,
35.77552476978808
],
[
139.56721442796606,
35.77552476978808
]
]
]
}
min_lon, max_lon = bbox["coordinates"][0][0][0], bbox["coordinates"][0][2][0]
min_lat, max_lat = bbox["coordinates"][0][1][1], bbox["coordinates"][0][0][1]
ここではmin_lon, max_lon, min_lat, max_latで四隅の座標値を取得しています。次のクエリを利用して、BBOX（長方形範囲）内に完全に収まる建物をフィルタリングしながらダウンロード・保存します。
データはGeoPackage形式で保存します。
query = f"""
COPY (
SELECT id, geometry
FROM read_parquet('{data_path}', filename=true, hive_partitioning=1)
WHERE bbox.xmin > {min_lon} AND bbox.xmax < {max_lon}
AND bbox.ymin > {min_lat} AND bbox.ymax < {max_lat}
) TO './data/overture_data.gpkg'
WITH (FORMAT GDAL, DRIVER 'GPKG');
"""
df = con.execute(query).df()
print(df.head())
今回は1631984件ものデータが得られました。
これで、overture_data.gpkgが生成され、指定した範囲内に含まれる建物データのみがGeoPackageとしてローカルにダウンロードされます。
QGISで建物データを可視化
生成したoverture_data.gpkgは、そのままQGISなどのGISソフトで読み込むことができます。
この時点で、対象エリアにある多数の建物形状を地図上に可視化でき、Overture Mapsの建物データを直接目視で確認することが可能です。
建物の重心点を算出
まずは建物ポリゴンの重心点を計算し、そのポイントを含むGeoDataFrameを作成します。建物の解析をする際、ポリゴンそのものを扱うより、重心点（などの代表点）を使ったほうが扱いやすいケースはよくあります。
import geopandas as gpd
# GeoPackageから建物データを読み込み、CRSをEPSG:4326に変換（必要に応じて）
gdf = gpd.read_file('data/overture_data.gpkg').to_crs("EPSG:4326")
# 重心点を算出し、'centroid'カラムとして追加
gdf['centroid'] = gdf.geometry.centroid
# 建物IDと重心点のみを抽出し、重心をgeometryとして扱うGeoDataFrameを作成
build_points = gpd.GeoDataFrame(
gdf[['id', 'centroid']].rename(columns={'centroid': 'geometry'}),
geometry='geometry',
crs=gdf.crs
)
# ポイント数
print(f"ポイント数: {len(build_points)}")
# ポイントのみをGeoPackageに保存
build_points.to_file('data/overture_buildings_centroid.gpkg', layer='points', driver='GPKG')
build_pointsには、id（オリジナルの建物ID）とgeometry（ポイント）が格納されます。後段でH3の六角形と空間演算を行うために、建物IDやその他属性をポイントと一緒に持っておくとなにかと便利です。
ポイントデータはGeoPackageで保存しておきます。
（黄色いのが建物の重心点です。）
H3の六角形を表示してみる
次に、H3インデックスを生成します。まずは「解像度7で指定したBBOXの範囲をすべてカバーする六角形セル群」を作成し、GeoDataFrame化します。
H3にはgeo_to_h3shapeやh3shape_to_cellsといった関数があり、GeoJSONポリゴン(今回のBBOX)をまとめて六角形セルのリストに変換可能です。以下はPythonのh3-pyを使った例です。
import h3
import pandas as pd
from shapely.geometry import Polygon
# Overture Mapsダウンロード時のものと同一
bbox = {
"type": "Polygon",
"coordinates": [
[
[
139.56721442796606,
35.77552476978808
],
[
139.56721442796606,
35.5781480840387
],
[
139.9653733196339,
35.5781480840387
],
[
139.9653733196339,
35.77552476978808
],
[
139.56721442796606,
35.77552476978808
]
]
]
}
# bboxをH3のポリゴンに変換
poly_shape = h3.geo_to_h3shape(bbox)
# セルの解像度は7を指定
resolution = 7
cells_in_bbox = h3.h3shape_to_cells(poly_shape, resolution)
print(f"セル数: {len(cells_in_bbox)}")
# 六角形セルをShapelyポリゴンに変換
def h3_to_polygon(h3_index):
boundary = h3.cell_to_boundary(h3_index) # [(lat, lng), ...]のリスト
polygon_points = [(lng, lat) for lat, lng in boundary]
return Polygon(polygon_points)
hexagon_polygons = [h3_to_polygon(h) for h in cells_in_bbox]
# gdf化し、列名にH3インデックスを保持する
h3_stats = pd.DataFrame({'h3_7_index': cells_in_bbox})
h3_stats['geometry'] = hexagon_polygons
h3_gdf = gpd.GeoDataFrame(h3_stats, geometry='geometry', crs='EPSG:4326')
# データを保存
h3_gdf.to_file('data/h3_grid_res7.gpkg', layer='grid_res7', driver='GPKG')
ここではカラム名にh3_7_indexを使い、どの解像度で生成したインデックスかを明示的に区別しています。QGIS等で扱う際にも、h3_7_indexが六角形セルを特定する属性としてわかりやすいです。
再度QGISで可視化
出力したはQGISにインポートすれば、解像度7のH3六角形メッシュが、先ほどの建物データと重なり合う形で可視化できます。
インデックス内部に含まれるポイントの数を属性に格納する
GeoPandasのsjoinを利用し、「ポイントがどの多角形（六角形）に属するか」を判定していきます。
# ポイント(建物重心)を六角形へ空間結合
joined_gdf = gpd.sjoin(build_points, h3_gdf, how='left', predicate='within')
# 六角形セルIDごとにポイント数を集計
count_by_hex = joined_gdf.groupby('h3_7_index')['id'].count().reset_index()
count_by_hex.rename(columns={'id': 'points_count'}, inplace=True)
# H3のgdfに集計結果をジョイン
h3_count_gdf = h3_gdf.merge(count_by_hex, on='h3_7_index', how='left')
h3_count_gdf['points_count'] = h3_count_gdf['points_count'].fillna(0).astype(int)
# 結果をGeoPackageに保存
h3_count_gdf.to_file('data/h3_building_count.gpkg', layer='building_count', driver='GPKG')
h3_count_gdf の各レコード（六角形セル）に、内部に含まれる重心ポイント数が points_count として紐づきます。
ポイント数に応じて可視化する
QGISではpoints_countを元に分類表示できるので、各六角形セルごとの建物数を直感的に把握できます。
試しにポイント数が多い六角形ほど色を濃くするとこんな感じになります。
解像度を変更し、集約してみる
H3は六角形で見た目が映え、ビジュアライズに向いていますが、そもそもは空間インデックスなので空間的な集約が得意です。
解像度を柔軟に変更できるのも特徴の一つです。
「解像度7で細かく分割された六角形」を、より粗い「解像度5の六角形」にまとめる、といったことをやってみましょう。cell_to_parentなどの関数を使って簡単に変換できます。
# 解像度7の各セルを解像度5の親セルに変換する
cells_res7 = h3_count_gdf['h3_7_index'].unique().tolist()
cells_res5 = [h3.cell_to_parent(cell, 5) for cell in cells_res7]
# 解像度7のセルから解像度5のセルへの対応関係を辞書として作成
# このマッピングは後で各ポイントの所属する解像度5セルを特定するために使用される
res7_to_5 = dict(zip(cells_res7, cells_res5))
unique_res5 = list(set(cells_res5))
hexagon_polygons_res5 = [h3_to_polygon(h) for h in unique_res5]
# 解像度5セルとそれに対応するジオメトリを含むデータフレームを作成
h3_res5_df = pd.DataFrame({'h3_5_index': unique_res5})
h3_res5_df['geometry'] = hexagon_polygons_res5
h3_res5_gdf = gpd.GeoDataFrame(h3_res5_df, geometry='geometry', crs='EPSG:4326')
# 元の解像度7のデータフレームに、各セルが属する解像度5セルの情報を追加
h3_count_gdf['h3_5_index'] = h3_count_gdf['h3_7_index'].apply(lambda x: res7_to_5[x])
# 解像度5セルごとにポイント数を集計
# groupbyで解像度5セルごとにグループ化し、points_countカラムの合計を計算
sum_by_hex5 = h3_count_gdf.groupby('h3_5_index')['points_count'].sum().reset_index()
sum_by_hex5.rename(columns={'points_count': 'points_count_res5'}, inplace=True)
# 解像度5のジオデータフレームに、集計したポイント数情報をマージ
h3_res5_count_gdf = h3_res5_gdf.merge(sum_by_hex5, on='h3_5_index', how='left')
h3_res5_count_gdf['points_count_res5'] = h3_res5_count_gdf['points_count_res5'].fillna(0).astype(int)
# 結果を保存
h3_res5_count_gdf.to_file('data/h3_building_count_res5.gpkg', layer='building_count_res5', driver='GPKG')
これで、解像度5 の六角形上にポイント数を再集計できます。集計粒度が粗くなるため、1つの六角形あたりのポイント数は解像度7より大きくなります。この h3_res5_count_gdf もGeoPackageに書き出してQGISで可視化すれば、解像度7との比較が可能です。
隣接するセルのポイント数を集計する
個別の地点を指定して「解像度7のH3インデックス」を取得し、そのセルに含まれる建物数を集計するようなこともできます。
tokyo_station = (35.681236, 139.767125)
res7_index_tokyo = h3.latlng_to_cell(tokyo_station[0], tokyo_station[1], 7)
print(f"東京駅付近のH3セルID (解像度7): {res7_index_tokyo}")
# 既に計算済みのh3_count_gdf(解像度7)から該当セルのpoints_countを見る
row = h3_count_gdf[h3_count_gdf['h3_7_index'] == res7_index_tokyo]
if not row.empty:
print("東京駅のある六角形のポイント数:", row.iloc[0]['points_count'])
else:
print("東京駅付近のセルが見つかりませんでした。")
このようにして、任意の地点について「どのH3セルに属するか」「そこに何ポイント含まれるか」を問い合わせられます。
今回は、以下のように表示されました。
東京駅付近のH3セルID (解像度7): 872f5a32dffffff
東京駅のある六角形のポイント数: 6235
最後に、grid_ringやgrid_diskを使って近傍セルを取得し、それらが含むポイント総数を合計してみます。
# 1セル離れた隣接セル
neighbors = h3.grid_ring(res7_index_tokyo, 1)
print("隣接セルID一覧:", neighbors)
# 隣接セルも含めてポイント数を合算
sel_cells = [res7_index_tokyo] + list(neighbors)
sub_gdf = h3_count_gdf[h3_count_gdf['h3_7_index'].isin(sel_cells)]
total_points = sub_gdf['points_count'].sum()
print(f"東京駅セル + 隣接セルのポイント数合計: {total_points}")
出力は以下のようになりました。
隣接セルID一覧: ['872f5aad3ffffff', '872f5a32cffffff', '872f5a328ffffff', '872f5a329ffffff', '872f5aadaffffff', '872f5aadeffffff']
東京駅セル + 隣接セルのポイント数合計: 62545
このように、H3を使えば「あるセルを中心とした一定距離の範囲」に含まれるポイントの数などを容易に集計でき、「近隣検索」や「周辺集計」にとても便利です。
まとめ
今回は、DuckDBを使ってOverture Mapsから建物データをダウンロードしたり、GeoPandasとH3を組み合わせて六角形メッシュ解析を行う一連の流れをやってみました！ポイントは次のとおりです。
H3は全世界規模で一意に定まるセルを利用できるため、可視化や空間分析にとても便利です。
特に近隣関係や解像度変更を簡単に扱えるという大きなメリットがあります。
ぜひ今回の流れを参考に、さまざまな場所やデータセットでH3インデックスを利用した解析を試してみてください！