はじめに

PMTilesやFlatGeobufなどの現代的な地理空間データフォーマットでは、ヒルベルト曲線(Hilbert curve) が空間インデックスとして利用されています。

他方、PMTilesやFlatGeobufを使われたことがあるという方でも、インデックスの仕組みまで理解しているという方は多くはないのではないでしょうか（そんなことないかもですが。。）。

かくいう私自身、インデックスの詳細を理解せずにこれらのフォーマットを使っていたのですが、仕事上の関係で調べる機会があり、実際に調べてみると結構興味深い内容だったので、今更かもしれませんが、記事としてまとめてみることにしました。

この記事を通して、私と同じようによく知らなかった方が、少しでも面白いと感じていただけたら幸いです。なお、専門家というわけではないので、記述に数学的な厳密性に欠けている部分もあるかと思いますが、ご容赦ください（あくまでイメージだけ掴んでいただければと思います）。

ヒルベルト曲線とは

定義

ヒルベルト曲線は、ドイツの数学者ダフィット・ヒルベルトによって1891年に考案された空間充填曲線(Space-filling curve) の一種です。1次元の線が2次元空間を隙間なく埋め尽くすという、一見すると直感に反する性質を持っています。

空間充填曲線の特徴

空間充填曲線には以下の性質があります。

連続性: 曲線は途切れることなく連続している
充填性: 反復を繰り返すことで、2次元空間のすべての点を通過する
自己相似性: 分割すれば同じ形が現れる（フラクタル特性）
局所性: 2次元空間で近い点は、1次元の曲線上でも比較的近くに配置される

自己相似性により実装がアルゴリズム的に扱いやすくなるため、情報技術分野で有用となります。また、局所性があることから、地理空間データのインデックスとして有効といえます。

ヒルベルト曲線の構造

基本パターン

ヒルベルト曲線は、再帰的な構造を持っています（自己相似性）。最も単純な1次のヒルベルト曲線から始めて、次数を増やしていくことでより複雑な曲線を生成します。

1次ヒルベルト曲線

  0   3
  |   |
  1---2

4つの正方形領域を訪れる、コの字型の単純なパターンです（向きに決まりはありませんが、ここでは上向きのコの字で説明します）。

2次ヒルベルト曲線

  0---1  14---15
      |   |
  3---2  13---12
  |           |
  4   7---8   11
  |   |   |   |
  5---6   9---10

1次の曲線を4つ組み合わせ、それぞれを回転・反転させて接続します。

n次ヒルベルト曲線の生成規則

n次の曲線は、(n-1)次の曲線を4つ使って構成されます。

左上: (n-1)次の曲線を90度左回転
左下: (n-1)次の曲線をそのまま配置
右下: (n-1)次の曲線をそのまま配置
右上: (n-1)次の曲線を90度右回転

これらを適切に接続することで、より高次の曲線が生成されます。

座標とインデックスの変換

ヒルベルト曲線の最も重要な応用は、2次元座標(x, y)を1次元のインデックスに変換する(写像する) ことです。

変換の仕組み

n次のヒルベルト曲線では、2^n × 2^n のグリッドを 0 から 2^(2n) - 1 までの整数にマッピングします。

例: 2次の曲線(4×4グリッド)

(0, 0) → 0
(1, 0) → 1
(1, 1) → 2
(0, 1) → 3
...
(3, 0) → 15

PMTilesにおけるヒルベルト曲線

PMTilesとは

PMTilesは、地図タイルを単一ファイルにまとめて効率的に配信するためのアーカイブ形式です。
HTTPのRange Requestを活用し、必要なタイルだけを取得できます。

PMTilesの仕様の詳細はPMTiles Specificationを参照ください。

PMTilesのタイルIDにヒルベルト曲線を利用

PMTilesのディレクトリインデックス（タイルID）でヒルベルト曲線が利用されています。

ディレクトリとは、タイルデータへのポインタです。PMTilesからタイルデータを取り出す際には、まず、ディレクトリを取得して、どこに目的のタイルデータがあるのかを調べる必要があります。

おおまかに以下のような手順でタイルデータを取得することになります（リーフディレクトリがある場合はもうひと手間必要ですがここでは割愛します）。

ディレクトリの取得: ディレクトリをHTTP Range Requestで取得
タイルIDを算出: 取得したいタイルのタイル座標（z/x/y）からタイルIDを計算
タイルエントリの位置取得: ディレクトリから該当のタイルIDのエントリを探索して、タイルデータのOffsetとLengthを取得
タイルエントリの取得: 目的のタイルデータをHTTP Range Requestで取得

ズームレベルzのタイルIDは、z次のヒルベルト曲線の順序でタイルIDとして付番され、z0からの通し番号になっています。

例として、z2までのXYZタイル座標とタイルIDの関係は以下のとおりです。z1のタイル座標は1次のヒルベルト曲線順、z2のタイル座標は2次のヒルベルト曲線順に付番されています。

| z | X | Y | TileID |
|---|---|---|--------|
| 0 | 0 | 0 | 0 |
| 1 | 0 | 0 | 1 |
| 1 | 0 | 1 | 2 |
| 1 | 1 | 1 | 3 |
| 1 | 1 | 0 | 4 |
| 2 | 0 | 0 | 5 |
| 2 | 1 | 0 | 6 |
| 2 | 1 | 1 | 7 |
| 2 | 0 | 1 | 8 |
| 2 | 0 | 2 | 9 |
| 2 | 0 | 3 | 10 |
| 2 | 1 | 3 | 11 |
| 2 | 1 | 2 | 12 |
| 2 | 2 | 2 | 13 |
| 2 | 2 | 3 | 14 |
| 2 | 3 | 3 | 15 |
| 2 | 3 | 2 | 16 |
| 2 | 3 | 1 | 17 |
| 2 | 2 | 1 | 18 |
| 2 | 2 | 0 | 19 |
| 2 | 3 | 0 | 20 |

なお、TileID >> (z,x,y)の変換はpythonで以下のとおり簡単に確認することができます。

pip install pmtiles
python
>> from pmtiles.tile import tileid_to_zxy
>> tileid_to_zxy(20)
(2, 3, 0)

なぜヒルベルト曲線なのか

PMTilesでは、タイルをタイルID順にソートしたときに、連続したタイルが同一タイルデータになる場合、それらを一つのエントリでまとめて管理する仕組みが採用されています（Run-Length Encoding）。

タイルIDにヒルベルト曲線を用いることで、空間的に近いタイルが連続することになり、同じタイルデータが連続している領域では、一組のディレクトリインデックスとタイルデータのみで管理することができます。

そのため、海ポリゴンなど同じタイルデータが連続する領域が多く含まれるタイルセットでは、データ容量及びデータアクセス効率が大幅に向上します。

例えば、筆者が、Natural Earthの1:10m Oceanについて、z0-8のPMTileを生成したところ、以下のディレクトリエントリが含まれていました。

TileID=40192,Offset=193094,Length=108,RunLength=1782

このように複数のタイルIDで同じタイルデータを参照する場合、RunLengthは2以上の値になっています（つまりRunLength分だけ同じタイルデータをみてねの意味です）。具体的には以下の意味になります。

タイルIDが40192-41973までは同じタイルデータを参照
- 40192=TileID
- 41973=TileID+RunLength-1
タイルデータは、タイルデータセクションの193094-193201バイトに格納
- 193094=Offset
- 193201=Offset+Length-1

これによりTileIDが40193-41973のエントリはディレクトリ及びタイルデータともに省略されます。
より詳細についてはprotomapsの以下のブログ記事を参照ください。

PMTiles version 3: Hilbert Tile IDs and Run-Length Encoding

FlatGeobufにおけるヒルベルト曲線

ヒルベルト曲線が利用されているもう一つの例としてFlatGeobufでの活用方法を紹介したいと思います。

FlatGeobufとは

FlatGeobufは、GeoJSONやShapefileに代わる、Cloud-Optimizedなベクトル地理空間データフォーマットです。HTTPのRange Requestを活用し、大規模なデータセットから必要な部分だけを効率的に取得できます。

FlatGeobufの仕様の詳細はFlatGeobuf#Specificationを参照ください。

Packed Hilbert R-Treeの構造

FlatGeobufでは、Packed Hilbert R-Treeという空間インデックスを採用しています。これは、ヒルベルト曲線とR木の利点を組み合わせた手法です。

R木の説明は割愛しますが、ご存知ない方は、PostGISの空間インデックス解説ページを参照ください。

基本的な仕組み

地物の並び替え: すべての地物の境界ボックス（bbox）の中心点を計算し、その中心点をヒルベルト曲線に沿ってソートします。
リーフノード構築: ソートされた地物のbboxを元にR木のリーフノードを作成します。
内部ノード～ルートノード構築: リーフノードを作成順に固定数（インデックスノードサイズ）ずつグループ化して、R木の内部ノードを作成します。同様に下位ノードから上位ノードを作成していき、ノード数がインデックスノードサイズ以下になったところで、その上位ノードがルートノードになります。
静的インデックス: 一度構築されたツリーは変更しません（Packed）。

地物のヒルベルト順序:
Feature 0 (hilbert_value: 5)
Feature 1 (hilbert_value: 12)
Feature 2 (hilbert_value: 18)
Feature 3 (hilbert_value: 22)
...

階層構造:
         [Root Node]
        /    |    \
    [Node0][Node1][Node2] ...
     /  \    /  \   /  \
   [F0][F1][F2][F3][F4][F5] ...

なぜPacked Hilbert R-Treeなのか

地物をヒルベルト曲線順にソートした上でR-Treeインデックスを構築することで、各リーフノードが包含するフィーチャ同士が近接するため、ノード間のbboxの重複を小さく抑えられます。

これにより、効率的に対象のフィーチャを探索することが可能になります。

実際に、インデックス付きのFlatGeobufに対して、特定のbbox内の地物を検索する場合、以下のような流れになります。

インデックスの取得: ファイルのインデックスセクションをHTTP Range Requestで取得
ツリーの探索: ルートノードから開始し、クエリボックスと交差するノードのみを辿る
地物データの取得: リーフノードに到達したら、該当する地物のオフセットを取得
地物の取得: 必要な地物データのみをHTTP Range Requestで取得

なお、FlatGeobufの公式の仕様では、インデックスはオプション扱いとなっているため、実際にインデックスが作成されるかは、使用するツールごとに確認が必要です。例えば、GDALのogr2ogrで変換した場合はデフォルトでインデックスが作成されます。

日本の行政区域で構築するとどうなるか

実際に日本の行政区域データ（国土数値情報）を使ってFlatGeobufを作成するとどうなるか試してみました。

わかりやすくするため、都道府県のデータを、小さな離島などを削除して簡略化したデータに加工し、ogr2ogrを使ってFlatGeobufに変換します。

ogr2ogr -f FlatGeobuf prefecture.fgb prefecture.geojson

pythonを使ってインデックスを解析すると以下のようになりました。

緑で網掛けになっているところが各インデックスノードのbboxです。

白で縁取りした番号がインデックスセクション内での格納順で、ヒルベルト曲線の逆順になっています。なお、インデックスノードサイズは16です。

ルートノード：1個

内部ノード：3個

リーフノード（＝フィーチャのbbox）：47個

おおよそ北東から南西にソートされていることがわかります。ヒルベルト曲線順でいうと南西から北東です。三重県が例外的に最後（格納順でいうと最初）になっているなど少し不可解な順番になっているところもありますが、全体的に近い番号同士で近接はしているのかなと思います。

なお、DuckDBでもヒルベルトインデックスを計算できます。結果は以下のとおりですが、計算に使われる正規化範囲やビット深度（○次のヒルベルト曲線を使うか）の違いによるものか上記と順番に差異があります（こちらの方がきれいにソートされているように見えます）。

INSTALL spatial;
LOAD spatial;

CREATE TABLE prefecture AS SELECT * FROM ST_Read('prefecture.geojson');

WITH bounds AS (
  SELECT ST_Extent(ST_MakeEnvelope(122, 20, 154, 46)) AS b
)
SELECT
  name,
  ST_Hilbert(geom, (SELECT b FROM bounds)) AS hilbert_index
FROM prefecture
ORDER BY hilbert_index DESC;

┌──────────┬───────────────┐
│   name   │ hilbert_index │
│ varchar  │    uint32     │
├──────────┼───────────────┤
│ 北海道   │    2610209600 │
│ 青森県   │    2446149408 │
│ 岩手県   │    2437809110 │
│ 秋田県   │    2433158149 │
│ 新潟県   │    2397555018 │
│ 群馬県   │    2389472842 │
│ 栃木県   │    2387357905 │
│ 福島県   │    2378561404 │
│ 茨城県   │    2376528597 │
│ 山形県   │    2358893222 │
│ 宮城県   │    2357525002 │
│ 千葉県   │    2185530984 │
│ 神奈川県 │    2178129160 │
│ 埼玉県   │    2174581310 │
│ 東京都   │    2173293723 │
│ 山梨県   │    2172160600 │
│ 長野県   │    2169822273 │
│ 静岡県   │    2167859539 │
│ 三重県   │    2136066287 │
│ 愛知県   │    2128867184 │
│   ·      │         ·     │
│   ·      │         ·     │
│   ·      │         ·     │
│ 奈良県   │    2091252816 │
│ 和歌山県 │    2089862587 │
│ 高知県   │    2078347388 │
│ 香川県   │    2075111694 │
│ 広島県   │    2070326293 │
│ 愛媛県   │    2066155682 │
│ 大分県   │    2060012896 │
│ 佐賀県   │    2057782560 │
│ 福岡県   │    2055865975 │
│ 山口県   │    2047532265 │
│ 島根県   │    2022068621 │
│ 岡山県   │    2018649279 │
│ 鳥取県   │    2017596435 │
│ 石川県   │    1905076365 │
│ 富山県   │    1902207763 │
│ 長崎県   │    1165215354 │
│ 熊本県   │     801850623 │
│ 宮崎県   │     795981768 │
│ 鹿児島県 │     783015464 │
│ 沖縄県   │     113487712 │
├──────────┴───────────────┤
│    47 rows (40 shown)    │
└──────────────────────────┘

他の空間充填曲線との比較

最後に他の空間充填曲線にも触れておきたいと思います。ここでは、しばしばヒルベルト曲線と比較されるZ階級曲線（モートン順序ともいう）を取り上げます。

Z階級曲線(Z-order curve)

Z階級曲線も空間充填曲線の一種で、平面を4分割してZになる順序で付番します。

Z階級:                 ヒルベルト:

1次
  0---1                 0   3
  |   |                 |   |
  2---3                 1---2

2次
  0---1   4---5         0---1  14---15
  |   |   |   |             |   |   
  2---3   6---7         3---2  13---12
                        |           |
  8---9  12---13        4   7---8   11
  |   |   |   |         |   |   |   |
  10--11 14---15        5---6   9---10

ビットに規則性があり実装が比較的簡単にできます。ビット演算方法の詳細はZ曲線での空間充填とモートンオーダーについてを参照ください。

Z階級曲線の特徴

利点: ビット演算で簡単に計算できる
欠点: ヒルベルト曲線に比べて局所性が劣る

空間局所性の差が、大規模な地図データでは顕著なデータ効率性の差を生むため、地理空間データフォーマットではヒルベルト曲線の方が採用されているものと考えられます。

都道府県をZ階級でソートするとどうなるか

さきほど使った国土数値情報の行政区域データ由来の都道府県ポリゴンを、Z階級でもソートしてみようと思います。DuckDBのコミュニティ拡張のlindelを使うと簡単に計算できます。

INSTALL lindel FROM community;
LOAD lindel;

WITH bounds AS (
    SELECT 122.0 AS minx, 20.0 AS miny,
           154.0 AS maxx, 46.0 AS maxy
)
SELECT
  name,
  morton_encode([
    ST_X(ST_Centroid(geom)) - bounds.minx
         / (bounds.maxx - bounds.minx) * 65535,
    ST_Y(ST_Centroid(geom)) - bounds.miny
         / (bounds.maxy - bounds.miny) * 65535
  ]::integer[2]) AS zcode
FROM prefecture, bounds
ORDER BY zcode DESC;

┌──────────┬──────────────────────┐
│   name   │        zcode         │
│ varchar  │        uint64        │
├──────────┼──────────────────────┤
│ 北海道   │ 18446744026721128970 │
│ 青森県   │ 18446744026721126227 │
│ 岩手県   │ 18446744026721126226 │
│ 秋田県   │ 18446744026721126224 │
│ 宮城県   │ 18446744026721126215 │
│ 山形県   │ 18446744026721126213 │
│ 福島県   │ 18446744026721126212 │
│ 茨城県   │ 18446744026721126209 │
│ 栃木県   │ 18446744026721126209 │
│ 千葉県   │ 18446744026721126208 │
│ 新潟県   │ 18446744026721125870 │
│ 群馬県   │ 18446744026721125867 │
│ 埼玉県   │ 18446744026721125866 │
│ 東京都   │ 18446744026721125866 │
│ 神奈川県 │ 18446744026721125866 │
│ 長野県   │ 18446744026721125865 │
│ 山梨県   │ 18446744026721125864 │
│ 富山県   │ 18446744026721125859 │
│ 石川県   │ 18446744026721125859 │
│ 岐阜県   │ 18446744026721125858 │
│   ·      │           ·          │
│   ·      │           ·          │
│   ·      │           ·          │
│ 愛知県   │ 18446744026721125815 │
│ 三重県   │ 18446744026721125813 │
│ 奈良県   │ 18446744026721125813 │
│ 大阪府   │ 18446744026721125791 │
│ 和歌山県 │ 18446744026721125790 │
│ 岡山県   │ 18446744026721125789 │
│ 香川県   │ 18446744026721125789 │
│ 徳島県   │ 18446744026721125788 │
│ 広島県   │ 18446744026721125783 │
│ 愛媛県   │ 18446744026721125782 │
│ 高知県   │ 18446744026721125782 │
│ 山口県   │ 18446744026721125695 │
│ 大分県   │ 18446744026721125694 │
│ 福岡県   │ 18446744026721125692 │
│ 佐賀県   │ 18446744026721125692 │
│ 熊本県   │ 18446744026721125691 │
│ 宮崎県   │ 18446744026721125691 │
│ 鹿児島県 │ 18446744026721125688 │
│ 長崎県   │ 18446744026721125686 │
│ 沖縄県   │ 18446744026721125470 │
├──────────┴──────────────────────┤
│ 47 rows (40 shown)    2 columns │
└─────────────────────────────────┘

結構いい感じになりましたw。このくらいの粒度のデータだとそこまでの差異はなさそうです。

まとめ

ヒルベルト曲線は、以下の理由でPMTilesやFlatGeobufなどの現代的（Cloud Optimized）な地理空間データフォーマットに利用されています。

空間局所性の保持: 2次元空間で近い点を1次元の曲線上でも近くに配置
効率的なデータ管理: 地理的に近いデータのエントリが連続することで効率的にデータを管理
効率的なデータアクセス: 空間的に近接する地物をグルーピングすることで効率的なデータ探索が可能

いかがだったでしょうか。空間的な2次元配置を1次元にソートできるというところで、様々な応用の可能性を感じられたのではないでしょうか。

参考資料

Hilbert curve - Hilbert curveの概要（Wikipedia）
PMTiles Specification - PMTiles v3のフォーマット仕様
PMTiles version 3: Hilbert Tile IDs and Run-Length Encoding - PMTilesにおけるヒルベルトタイルID及びRunLengthの詳細
FlatGeobuf Specification - FlatGeobufのフォーマット仕様
Packed Hilbert R-Trees - Packed Hilbert R-Treeの概要（Wikipedia）
Space-filling curve - 空間充填曲線の概要（Wikipedia）