国土数値情報土砂災害警戒区域などデータを PMTiles + FlatGeobuf + DuckDB で軽く表示してみた

Posted at 2026-09-17

国土数値情報土砂災害警戒区域などデータを PMTiles + FlatGeobuf + DuckDB で軽く表示してみた

前回の記事では、国土数値情報のA33（土砂災害警戒区域）とA46（地すべり防止区域）をDuckDB Spatialへ取り込み、FastAPIでMVTを生成してOpenLayersへ表示しました。

動的MVTは検索条件を変えやすく、更新にも対応しやすい方法です。ただし、最初に開かれたタイルはDuckDBで抽出、クリップ、MVT変換まで行うため、全国規模になると初回表示の負荷が気になります。キャッシュが効いた後は軽くても、利用者が増えたときのサーバー負荷は残ります。

そこで今回は、地図表示をすべて動的MVTへ任せるのをやめました。

低ズーム       PMTiles（20km件数メッシュ）
中ズーム       PMTiles（20m LOD）
高ズーム       FlatGeobuf（元形状）
比較・診断     DuckDB動的MVT（100m / 20m / 元形状）
検索・詳細     FastAPI + DuckDB Spatial

FlatGeobufだけに切り替えるのではなく、縮尺に応じて役割を分けています。低ズームで細かな区域を全部描くことを避けた点が、今回の主な変更です。

今回作ったもの

構成はDocker Composeです。

前処理：Python 3.13、DuckDB Spatial、Tippecanoe
API：FastAPI、DuckDB Spatial
Web：OpenLayers、Vite、Nginx
静的地図：PMTiles、FlatGeobuf
比較用地図：DuckDBから生成する動的MVT

全体は次の流れになります。

通常表示は静的ファイルなので、地図を動かすたびにサーバーで空間SQLを実行しません。住所、告示日、区域番号などの詳細属性はDuckDBに残し、区域をクリックしたときだけFastAPIから取得します。

FlatGeobufだけにしなかった理由

FlatGeobufには空間インデックスを持たせられます。HTTP Range Requestを使い、画面に見えている範囲と交差する地物だけを取得できるため、GeoJSONを丸ごと読む方法よりかなり扱いやすくなります。

一方、FlatGeobufは元形状を保持する形式です。MVTのように、ズームごとに形状を簡略化したタイルがあらかじめ並んでいるわけではありません。

たとえば日本全体を表示しているときにFlatGeobufを使うと、表示範囲は日本全体です。その範囲に入る区域を大量に取得し、ブラウザー側でデコードして描画することになります。通信が減っても、JavaScript処理や描画が重くなる場合があります。

今回の使い分けは次のとおりです。

表示範囲	使用形式	内容
全国・地方	PMTiles	区域重心を20kmセルに集計した件数メッシュ
都道府県・市区町村	PMTiles	20m許容距離で事前簡略化した区域
区域の詳細	FlatGeobuf	元形状を表示範囲だけ取得
性能比較・条件追加の土台	動的MVT	DuckDBでタイル生成しファイルキャッシュ

FlatGeobufは高ズームで使い、PMTilesは広域表示を受け持たせました。

A33とA46の属性で気をつけたところ

A33は国土数値情報の「土砂災害警戒区域データ」です。今回の実装は第2.0版の属性名を基本にしています。

A33_001：現象の種類
A33_002：区域区分
A33_003：都道府県コード
A33_004：区域番号
A33_005：区域名
A33_006：所在地
A33_007：告示日
A33_008：特別警戒区域未指定フラグ

A46は「地すべり防止区域データ」です。A46-a、A46-b、A46-cは災害種別ではなく、所管省庁の違いです。

系列	所管
`A46-a`	国土交通省
`A46-b`	農林水産省農村振興局
`A46-c`	林野庁

入力ファイルによって使われる列系列が異なるため、取り込み時には同じ属性番号をCOALESCEでまとめています。

# A46は所管省庁ごとにa/b/cの列名を持つ。
# 同じ属性番号をまとめ、正規化後は共通列として扱う。
def _a46_candidates(suffix: str) -> tuple[str, ...]:
    """A46のa/b/c列と正規化済み別名を候補として返す。"""

    return (
        f"A46-a_{suffix}",
        f"A46-b_{suffix}",
        f"A46-c_{suffix}",
    )

A33/A46はいずれもJGD2011の緯度経度として読み、DuckDB内では次の2種類を保持します。

geom_4326：FlatGeobuf、GeoJSON、範囲表示用
geom_3857：距離を使うLOD作成、R-tree検索、MVT生成用

COPY ... SRS 'EPSG:4326'は出力ファイルへ座標参照系の情報を付ける指定であり、座標値そのものを変換する処理ではありません。座標変換は、取り込み時のST_Transformで明示しています。

ST_Transform(clean_geom, 'EPSG:6668', 'EPSG:4326', true) AS geom_4326,
ST_Transform(clean_geom, 'EPSG:6668', 'EPSG:3857', true) AS geom_3857

trueは経度・緯度の軸順を意識した指定です。入力データの座標系が異なる場合は、.envのSOURCE_CRSを変更します。

フォルダー構成

参考のため、ZIPファイルを置いておきますので、ダウンロードして活用してください。

sabo_hazard_ol_fastapi_duckdb_hybrid_r002/
├─ api/
│  ├─ app/
│  │  ├─ database.py
│  │  ├─ main.py
│  │  ├─ manifest.py
│  │  ├─ repository.py
│  │  ├─ settings.py
│  │  └─ tiles.py
│  ├─ tests/
│  ├─ Dockerfile
│  └─ requirements.txt
├─ pipeline/
│  ├─ app/
│  │  ├─ archive.py
│  │  ├─ build.py
│  │  ├─ cli.py
│  │  ├─ database.py
│  │  ├─ discovery.py
│  │  ├─ export.py
│  │  ├─ ingest.py
│  │  ├─ optimize.py
│  │  ├─ sample.py
│  │  ├─ settings.py
│  │  ├─ utils.py
│  │  └─ validate.py
│  ├─ Dockerfile
│  └─ requirements.txt
├─ web/
│  ├─ src/
│  │  ├─ layers.js
│  │  ├─ main.js
│  │  ├─ style.css
│  │  └─ styles.js
│  ├─ Dockerfile
│  ├─ index.html
│  ├─ nginx.conf
│  ├─ package.json
│  └─ vite.config.js
├─ scripts/
│  ├─ benchmark_dynamic_mvt.py
│  ├─ benchmark_http.py
│  ├─ build_real.ps1
│  ├─ build_real.sh
│  ├─ build_sample.ps1
│  ├─ build_sample.sh
│  ├─ smoke_test.ps1
│  ├─ smoke_test.sh
│  └─ static_check.py
├─ tests/
├─ data/
│  ├─ raw/A33/
│  ├─ raw/A46/
│  ├─ extracted/
│  ├─ sample/
│  ├─ duckdb/
│  ├─ work/
│  ├─ tile_cache/
│  └─ published/
├─ docs/
├─ .env.example
├─ docker-compose.yml
├─ Makefile
└─ README.md

pipelineは前処理専用です。APIコンテナにTippecanoeやデータ変換処理を入れず、公開中のサービスとビルド処理を分けています。

主なバージョン

依存関係は直接指定するものを固定しています。

Python         3.13
DuckDB         1.5.5
FastAPI        0.141.1
Tippecanoe     2.79.0
OpenLayers     9.2.4
ol-pmtiles     2.0.2
FlatGeobuf JS  4.4.0
Vite           8.2.0

ol-pmtilesとの組み合わせを不用意に変えないよう、OpenLayersは9.2.4へ固定しました。OpenLayersやol-pmtilesを更新する場合は、PMTilesの初回表示、パン、ズーム、Feature選択までまとめて確認します。

Docker Composeの分け方

前処理は常駐させず、toolsプロファイルで必要なときだけ実行します。APIから運用DuckDBを更新できないよう、DBと公開ファイルは読み取り専用でマウントしています。

services:
  pipeline:
    build:
      context: ./pipeline
      args:
        TIPPECANOE_VERSION: "2.79.0"
    env_file:
      - .env
    volumes:
      - ./data:/data
    profiles:
      - tools

  api:
    build:
      context: ./api
    env_file:
      - .env
    ports:
      - "${API_PORT:-8000}:8000"
    volumes:
      # 運用DBと公開マニフェストは前処理だけが更新する。
      - ./data/duckdb:/data/duckdb:ro
      - ./data/published:/data/published:ro

      # 比較用の動的MVTだけはAPIが生成する。
      - ./data/tile_cache:/data/tile_cache

  web:
    build:
      context: ./web
    ports:
      - "${WEB_PORT:-8080}:80"
    volumes:
      - ./data/published:/srv/hazard-data:ro
    depends_on:
      api:
        condition: service_healthy

サンプルデータで起動する

最初は同梱のサンプルを使います。サンプルは性能比較用ではなく、データ変換、公開ファイル、API、画面のつながりを確認するためのものです。

Linux / WSL / Ubuntu

cp .env.example .env
bash scripts/build_sample.sh

Windows PowerShell

Copy-Item .env.example .env
powershell -ExecutionPolicy Bypass -File .\scripts\build_sample.ps1

起動後は次を開きます。

WebMap      http://localhost:8080
Swagger UI http://localhost:8000/docs

処理を分けて実行する場合は次のとおりです。

# DuckDB、FlatGeobuf、PMTiles、manifestを作る。
docker compose --profile tools run --rm pipeline build --sample

# APIとWebを起動する。
docker compose up --build -d api web

国土数値情報のA33/A46を取得する

これまでは、公式ページをブラウザーで開き、ZIPを手作業でdata/rawへ置いていました。1県だけならそれでも困りませんが、複数県を試したり、毎年度更新したりすると、ファイル名と取得元を記録する作業が増えてきます。

そこで、r003ではA33/A46の公式カタログを読み、条件に合うZIPを選ぶスクリプトを追加しました。r004では、これを都道府県コード順の一括処理へ広げています。

scripts/download_ksj_hazard.py

配布URLそのものは年度更新で変わる可能性があります。スクリプトへ最新年度のURLを固定するのではなく、カタログに掲載されたファイル名、年度、地域、形式を読み取ってから取得します。

最初に取得対象だけ確認する

実データを保存する前に、listで選択結果を確認します。次は静岡県のA33/A46を選ぶ例です。

python scripts/download_ksj_hazard.py list \
  --dataset A33 A46 \
  --area-code 22 \
  --year latest \
  --format AUTO

都道府県名でも指定できます。

python scripts/download_ksj_hazard.py list \
  --dataset A33 A46 \
  --prefecture 静岡県

AUTOは既存の前処理に合わせた指定です。

データ	AUTOで優先する形式	理由
A33	Shapefile	属性名と形状を前処理で扱いやすい
A46	GML	公式カタログ上の配布単位に合わせる

形式を固定したい場合は、GML、SHP、GEOJSON、ALLを指定できます。

利用条件を確認してから取得する

A33は都道府県ごとに公開条件が異なる場合があります。A46も、国土数値情報の利用規約だけでなく、都道府県別の利用条件を確認する必要があります。スクリプトは利用可否を自動判定しません。確認後に--accept-termsを付けて実行します。

python scripts/download_ksj_hazard.py download \
  --dataset A33 A46 \
  --area-code 22 \
  --year latest \
  --format AUTO \
  --accept-terms

保存先は、後段の前処理がそのまま読める配置です。

data/raw/
├─ A33/
│  └─ A33-<年度>_<地域>_<形式>.zip
└─ A46/
   └─ A46-<年度>_<地域>_<形式>.zip

--accept-termsは、利用者が条件を確認したことをコマンド上で明示するためのフラグです。規約そのものを置き換えるものではありません。

中断したファイルは`.part`から再開する

大きな全国ZIPでは、途中で接続が切れることもあります。取得中は正式なZIP名にせず、.partを付けて保存します。

A33-25_00_SHP.zip.part

同じコマンドを再実行すると、保存済みサイズを確認してHTTP Rangeによる再開を試みます。サーバーがRangeを受け付けない場合は、先頭から取り直します。

この部分では、途中ファイルと完成ファイルを混同しないようにしています。

def _download_stream(
    url: str,
    *,
    destination_part: Path,
    timeout: float,
    resume: bool,
    referer: str,
) -> None:
    """1つのURLからZIPを``.part``へストリーミング保存する。

    Range要求に対してサーバーが200を返した場合は先頭から保存し直す。
    レスポンス先頭がZIPシグネチャでない場合は、HTMLエラーページなどと
    判断して失敗させる。
    """

ダウンロード後にZIPを検査する

HTTP 200が返っても、その内容が必ずZIPとは限りません。サイト側のエラーページや案内HTMLを保存してしまう場合もあるため、取得後に次を確認します。

ファイル先頭がZIPのシグネチャか
zipfile.is_zipfile()で開けるか
ZIP内のCRC検査を通るか
GISファイルが1件以上入っているか
Shapefileなら.shp、.dbf、.shxが同じ基底名で揃っているか
SHA-256を計算できるか

def validate_zip_file(path: Path, *, expected_format: str) -> ZipValidation:
    """ZIPの構造、CRC、GISファイル構成を検査する。

    Args:
        path: 検査するZIP。
        expected_format: カタログ上の形式。

    Raises:
        DownloaderError: ZIP破損、CRC不一致、GISデータ不在、
            Shapefile構成不足を検出した場合。
    """

検査に失敗したファイルは正式名へ置き換えません。正常なZIPだけをdata/raw/A33またはdata/raw/A46へ残します。

SHA-256と取得元を残す

取得結果はJSONとCSVへ保存します。

data/raw/download_manifests/
├─ ksj_hazard_download_<UTC時刻>.json
├─ ksj_hazard_download_<UTC時刻>.csv
└─ latest.json

マニフェストには、次を記録します。

データセット、地域コード、地域名
データ年度、形式、ファイル名
カタログ上のファイル容量
実際に使用したURL
試したURL候補
ローカルパスと実ファイル容量
SHA-256
ZIP内の項目数
検出したGISファイル
成功、既存利用、失敗などの状態

カタログHTMLも保存します。

data/raw/catalog_cache/
├─ A33/catalog.html
├─ A33/catalog.json
├─ A46/catalog.html
└─ A46/catalog.json

後から「どのカタログを見て、そのZIPを選んだのか」を追いやすくするためです。

OS別のラッパー

Linux、WSL、Ubuntuでは、利用条件を確認したうえで次を実行します。第1引数は地域コードです。

KSJ_ACCEPT_TERMS=YES bash scripts/download_hazard_data.sh 22

Windows PowerShellでは次のとおりです。

powershell -ExecutionPolicy Bypass `
  -File .\scripts\download_hazard_data.ps1 `
  -AreaCode 22 `
  -AcceptTerms

Makefileから実行する場合は、先に一覧を確認します。

make download-list AREA_CODE=22
make download-data AREA_CODE=22 ACCEPT_TERMS=YES

取得後に前処理まで続ける場合は次を使います。

make download-and-build AREA_CODE=22 ACCEPT_TERMS=YES

都道府県コード順に一括処理する

単県指定はそのまま残し、全都道府県とコード範囲の指定を追加しました。国土数値情報の
都道府県コードは、北海道01から沖縄県47までです。

全都道府県の選択結果だけを確認する場合は次を実行します。

python scripts/download_ksj_hazard.py list \
  --dataset A33 A46 \
  --all-prefectures \
  --year latest \
  --format AUTO

実際に取得する場合は、利用条件を確認して--accept-termsを付けます。

python scripts/download_ksj_hazard.py download \
  --dataset A33 A46 \
  --all-prefectures \
  --year latest \
  --format AUTO \
  --accept-terms \
  --continue-on-error

処理順は、データセットをまとめて処理するのではなく、次の順になります。

01 北海道  A33
01 北海道  A46
02 青森県  A33
02 青森県  A46
...
47 沖縄県  A33
47 沖縄県  A46

A46に対象県のデータが掲載されていない場合は、その県のA46だけを警告して飛ばし、
次の都道府県へ進みます。個別ZIPの取得失敗後も処理を続けたい場合は、
--continue-on-errorを付けておきます。

開始・終了コードを指定した部分一括処理もできます。

python scripts/download_ksj_hazard.py download \
  --dataset A33 A46 \
  --prefecture-code-range 02 07 \
  --format AUTO \
  --accept-terms \
  --continue-on-error

Linux、WSL、Ubuntu用ラッパーでは、allまたは01-07のような範囲を第1引数へ渡します。

KSJ_ACCEPT_TERMS=YES \
  bash scripts/download_hazard_data.sh all --continue-on-error

KSJ_ACCEPT_TERMS=YES \
  bash scripts/download_hazard_data.sh 02-07 --continue-on-error

PowerShellではスイッチまたは開始・終了コードを使います。

powershell -ExecutionPolicy Bypass `
  -File .\scripts\download_hazard_data.ps1 `
  -AllPrefectures `
  -AcceptTerms `
  -ContinueOnError

powershell -ExecutionPolicy Bypass `
  -File .\scripts\download_hazard_data.ps1 `
  -StartAreaCode 02 `
  -EndAreaCode 07 `
  -AcceptTerms `
  -ContinueOnError

01 49のように公式コード外まで指定しても、48・49は警告して除外します。
沖縄県の公式コードは47なので、49へ置き換えてダウンロードURLを組み立てることはしません。

全国版と全都道府県

A33全国版は--nationwideで選べます。

python scripts/download_ksj_hazard.py download \
  --dataset A33 \
  --nationwide \
  --format SHP \
  --accept-terms

全国版は大きいため、最初に1県分でDuckDB取り込み、PMTiles生成、画面表示まで確認してから切り替えます。

47都道府県を個別に取得する場合は次のとおりです。

python scripts/download_ksj_hazard.py download \
  --dataset A33 A46 \
  --all-prefectures \
  --format AUTO \
  --accept-terms \
  --continue-on-error

--continue-on-errorを付けると、一部地域で取得に失敗しても残りを続けます。失敗内容はマニフェストへ残ります。

A46の過年度をどう扱うか

A46の公式ページには複数年度が掲載されています。公式の注意事項では、過年度分も合わせて使い、重複地点は最新年度を使うよう案内されています。

ただし、過年度ZIPをそのまま全部取り込むと、同じ区域を二重登録する可能性があります。今回のダウンローダーは、安全側として既定値をlatestにしました。

python scripts/download_ksj_hazard.py download \
  --dataset A46 \
  --area-code 22 \
  --a46-history latest \
  --accept-terms

履歴を確認する目的で全年度を取得する場合だけ、allへ変更します。

python scripts/download_ksj_hazard.py download \
  --dataset A46 \
  --area-code 22 \
  --a46-history all \
  --accept-terms

この場合、前処理側で最新版優先の重複整理を追加してから運用データへ使います。

再取得、ドライラン、オフライン確認

既存の正常ZIPは再利用します。取り直す場合は--overwriteを付けます。

python scripts/download_ksj_hazard.py download \
  --dataset A33 A46 \
  --area-code 22 \
  --accept-terms \
  --overwrite

選択内容と候補URLだけを確認する場合は--dry-runを使います。ZIPは取得しません。

python scripts/download_ksj_hazard.py download \
  --dataset A33 A46 \
  --area-code 22 \
  --dry-run

一度保存したカタログだけで選択結果を再現する場合は--offlineを指定します。

python scripts/download_ksj_hazard.py list \
  --dataset A33 A46 \
  --area-code 22 \
  --offline

公式ページの構成が変わったときは、--catalog-url A33=...と--url-templateで一時的にURL規則を上書きできます。通常は、HTML内で見つかった公式URLを優先します。

ダウンロード後に前処理を実行する

ZIPは展開せず、次へ置かれています。

data/raw/A33/*.zip
data/raw/A46/*.zip

そのまま既存のビルドを実行します。

bash scripts/build_real.sh

PowerShellでは次のとおりです。

powershell -ExecutionPolicy Bypass -File .\scripts\build_real.ps1

ZIPはファイルごとに別フォルダーへ展開し、完了時に.completeを置きます。再実行時は展開済みデータを再利用します。

ZIP内に../などが含まれていても展開先の外へ書き出さないよう、展開前にパスを確認しています。

def _safe_extract(zip_path: Path, destination: Path) -> None:
    """ZIP Slipを防ぎながらアーカイブを展開する。

    Args:
        zip_path: 展開するZIP。
        destination: 展開先フォルダー。

    Raises:
        ValueError: ZIP内のパスが展開先の外側を指す場合。
    """

    destination.mkdir(parents=True, exist_ok=True)
    destination_root = destination.resolve()

    with zipfile.ZipFile(zip_path) as archive:
        for member in archive.infolist():
            output_path = (destination / member.filename).resolve()
            if not output_path.is_relative_to(destination_root):
                raise ValueError(
                    f"展開先の外側を指すZIP項目を検出しました: {member.filename}"
                )
        archive.extractall(destination)

入力形式はShapefile、FlatGeobuf、GeoJSON、GMLに対応しています。同じフォルダーに同名の複数形式がある場合は、二重登録を避けるため、次の優先順で1つだけ読みます。

Shapefile → FlatGeobuf → GeoJSON → GML

国土数値情報に付属する説明用JSONを誤ってGISデータとして拾わないよう、.jsonは対象にせず、.geojsonだけを読んでいます。

高速化1：元形状と表示用LODを分ける

DuckDBには元形状を残したまま、100mと20mのLODテーブルを作ります。

a33_hazard       元形状
├─ a33_lod100    低めのズームで使う比較用LOD
└─ a33_lod20     PMTilesと中ズーム動的MVT用

同じ構成をA46にも作成

LOD作成ではST_SimplifyPreserveTopologyを使います。ただし、小さなポリゴンや短い線は許容距離によって空になることがあります。防災区域を速度のために黙って消すのは避けたいので、空になった場合は元形状へ戻します。

def _create_lod_table(
    connection,
    *,
    source_table: str,
    output_table: str,
    columns: tuple[str, ...],
    tolerance_m: float,
) -> None:
    """元テーブルから表示専用の簡略形状テーブルを作成する。

    ``ST_SimplifyPreserveTopology`` で小さなポリゴンや短い線が空になった場合、
    そのFeatureを消さずに元形状へ戻す。防災区域では表示速度のために区域を
    黙って間引かないことを優先する。
    """

    selected_columns = _column_list(columns)
    source_identifier = quote_identifier(source_table)
    output_identifier = quote_identifier(output_table)
    tolerance_literal = format(float(tolerance_m), ".15g")

    connection.execute(
        f"""
        CREATE OR REPLACE TABLE {output_identifier} AS
        WITH candidates AS (
            SELECT
                {selected_columns},
                geom_3857 AS original_geom,
                CASE
                    WHEN geometry_kind = 'point' OR {tolerance_literal} <= 0
                        THEN geom_3857
                    ELSE ST_SimplifyPreserveTopology(
                        geom_3857,
                        {tolerance_literal}
                    )
                END AS simplified_geom
            FROM {source_identifier}
        )
        SELECT
            {selected_columns},
            CASE
                WHEN simplified_geom IS NULL OR ST_IsEmpty(simplified_geom)
                    THEN original_geom
                ELSE simplified_geom
            END AS geom_3857,
            (
                simplified_geom IS NULL OR ST_IsEmpty(simplified_geom)
            )::BOOLEAN AS simplification_fallback
        FROM candidates
        """
    )

LODテーブルには、地図表示とクリック判定に必要な属性だけを残します。

id
データセット名
区域区分・現象種別・所管
区域名
geometry_kind

所在地や告示番号などは元テーブルへ残します。全属性をMVTへ詰め込むと、タイルサイズが大きくなるためです。

高速化2：低ズームは件数メッシュへ切り替える

全国表示で個々の区域をすべて描いても、画面上では重なって見えます。そこで、低ズームは20kmセルの件数メッシュにしました。

各区域の重心をセルへ割り当て、次を集計します。

A33の件数
A33警戒区域の件数
A33特別警戒区域の件数
A46の件数
合計件数

これは区域形状ではなく、分布を見るための概観レイヤーです。区域の有無や境界を判断する用途には使いません。

低ズームの負荷を区域の間引きで逃がさず、別の表現へ切り替えています。

高速化3：元形状とLODへR-treeを作る

元形状、100m LOD、20m LOD、件数メッシュへR-treeを作ります。

CREATE INDEX IF NOT EXISTS idx_a33_hazard_geom
ON a33_hazard
USING RTREE (geom_3857);

CREATE INDEX IF NOT EXISTS idx_a33_lod100_geom
ON a33_lod100
USING RTREE (geom_3857);

CREATE INDEX IF NOT EXISTS idx_a33_lod20_geom
ON a33_lod20
USING RTREE (geom_3857);

A46も同様です。

R-treeを作っただけで、すべての空間SQLが必ずインデックス検索になるわけではありません。今回の動的MVTでは、タイル範囲をSQLの定数式として組み立てています。

def build_explain_sql(
    *,
    dataset: DynamicDataset,
    zoom: int,
    tile_x: int,
    tile_y: int,
    policy: DynamicMvtPolicy,
) -> tuple[str, str]:
    """R-tree利用状況を確認する ``EXPLAIN`` SQLを返す。"""

    validate_tile_coordinate(zoom, tile_x, tile_y, policy)
    spec = DYNAMIC_DATASET_SPECS[dataset]
    table_name = select_source_table(spec, zoom, policy)

    # タイル範囲をプレースホルダーへ渡さず、計画時に値が分かる式にする。
    envelope = f"ST_TileEnvelope({zoom}, {tile_x}, {tile_y})"
    sql = f"""
        EXPLAIN
        SELECT id
        FROM {_quote_identifier(table_name)}
        WHERE ST_Intersects(geom_3857, {envelope})
    """
    return sql, table_name

APIには実行計画の確認用URLも用意しました。

GET /api/debug/tiles/a33/12/3622/1612/explain?build_id=<build_id>

レスポンスのuses_rtree_index_scanがtrueか確認できます。タイル番号はデータ範囲に合わせて変更します。

FlatGeobufは元形状で出力する

FlatGeobufには、簡略化していないgeom_4326を出力します。

def _copy_to_flatgeobuf(connection, select_sql: str, output_path: Path) -> None:
    """SELECT結果を空間インデックス付きFlatGeobufとして保存する。"""

    output_path.parent.mkdir(parents=True, exist_ok=True)
    if output_path.exists():
        output_path.unlink()

    connection.execute(
        f"""
        COPY (
            {select_sql}
        )
        TO {sql_literal(output_path)}
        WITH (
            FORMAT GDAL,
            DRIVER 'FlatGeobuf',
            SRS 'EPSG:4326',
            LAYER_CREATION_OPTIONS 'SPATIAL_INDEX=YES'
        )
        """
    )

A33/A46は形状種別ごとに分けています。

fgb/
├─ a33_polygon.fgb
├─ a33_line.fgb
├─ a33_point.fgb
├─ a46_polygon.fgb
├─ a46_line.fgb
└─ a46_point.fgb

0件の形状種別は出力しません。FlatGeobufは1ファイル1レイヤーとして扱う方が分かりやすく、画面側でも必要なレイヤーだけ作れます。

PMTilesは20m LODから作る

PMTilesの入力は20m LODです。DuckDBから一度FlatGeobufへ出力し、Tippecanoeへ渡します。

DuckDB a33_lod20
    ↓
作業用FlatGeobuf
    ↓
Tippecanoe
    ↓
a33_polygon.pmtiles

Tippecanoeのコマンドは次のように組み立てています。

def _run_tippecanoe(
    settings: Settings,
    *,
    input_path: Path,
    output_path: Path,
    layer_name: str,
    title: str,
    min_zoom: int,
    max_zoom: int,
) -> None:
    """FlatGeobufから1レイヤーのPMTilesを作成する。

    防災区域を件数制限で間引かないよう、タイル内のFeature数とファイルサイズの
    上限を解除する。全国版ではタイルが大きくなる可能性があるため、開始ズームと
    LOD許容距離は計測結果を見ながら調整する。
    """

    command = [
        settings.tippecanoe_binary,
        "--force",
        f"--output={output_path}",
        f"--layer={layer_name}",
        f"--name={title}",
        "--projection=EPSG:4326",
        f"--minimum-zoom={min_zoom}",
        f"--maximum-zoom={max_zoom}",
        "--read-parallel",
        "--no-feature-limit",
        "--no-tile-size-limit",
        # 小さな区域を見た目だけの都合で縮小しない。
        "--no-tiny-polygon-reduction",
        # 隣接ポリゴンの共有境界を別々に単純化しない。
        "--no-simplification-of-shared-nodes",
        str(input_path),
    ]
    subprocess.run(command, check=True)

--no-feature-limitと--no-tile-size-limitを付けているため、低いズームから個別区域を入れるとタイルが大きくなります。既定値では個別区域PMTilesをズーム9から始め、さらに低いズームは件数メッシュへ任せています。

この設定は「必ず速い値」ではありません。対象地域、区域数、形状の細かさによって変わるので、実データで調整します。

公開ファイルをビルドIDで揃える

PMTilesやFlatGeobufへ長期キャッシュを付ける場合、同じURLの内容を上書きすると、旧データと新データが混ざることがあります。

そこで、公開先をビルドIDごとに分けました。

data/published/
├─ manifest.json
├─ builds/
│  └─ 20260807T001234Z/
│     ├─ manifest.json
│     ├─ file_catalog.json
│     ├─ fgb/
│     └─ pmtiles/
└─ reports/

画面が最初に読むのはdata/published/manifest.jsonです。そこから、ビルドID付きのURLを参照します。

{
  "build_id": "20260807T001234Z",
  "layers": [
    {
      "id": "a33_polygon",
      "flatgeobuf": {
        "url": "/data/builds/20260807T001234Z/fgb/a33_polygon.fgb"
      },
      "pmtiles": {
        "url": "/data/builds/20260807T001234Z/pmtiles/a33_polygon.pmtiles"
      }
    }
  ]
}

公開順序も決めています。

作業用DuckDBを作る
ビルドID別の静的ファイルを作る
件数、形状、インデックス、ハッシュを検査する
運用DuckDBを差し替える
最後にルートmanifest.jsonを差し替える

生成途中のファイルを画面から参照しないようにしています。厳密な分散トランザクションではありませんが、単一ホストの静的公開として、更新途中の混在をかなり減らせます。

動的MVTは比較用として残す

前の記事の動的MVTは削除せず、性能比較と将来の条件検索に使えるよう残しました。

ズームごとの参照テーブルは次のとおりです。

ズーム	参照テーブル
9～10	100m LOD
11～13	20m LOD
14～16	元形状

MVT生成SQLは次の形です。

def build_mvt_sql(
    *,
    dataset: DynamicDataset,
    zoom: int,
    tile_x: int,
    tile_y: int,
    policy: DynamicMvtPolicy,
) -> tuple[str, str]:
    """動的MVT生成SQLと参照するLODテーブル名を返す。

    タイル番号は整数として検査した後に定数式へ埋め込む。テーブル名、列名、
    MVTレイヤー名は固定定義から選び、URL文字列をSQL識別子には使わない。
    """

    validate_tile_coordinate(zoom, tile_x, tile_y, policy)
    spec = DYNAMIC_DATASET_SPECS[dataset]
    table_name = select_source_table(spec, zoom, policy)
    table_identifier = _quote_identifier(table_name)
    extent = int(policy.extent)
    buffer = int(policy.buffer)
    envelope = f"ST_TileEnvelope({zoom}, {tile_x}, {tile_y})"
    properties = ",\n                ".join(
        _quote_identifier(column) for column in spec.property_columns
    )
    struct_expression = _mvt_struct_expression(spec)

    sql = f"""
        WITH tile_rows AS (
            SELECT
                {properties},
                ST_AsMVTGeom(
                    geom_3857,
                    ST_Extent({envelope}),
                    {extent},
                    {buffer},
                    true
                ) AS geom
            FROM {table_identifier}
            WHERE ST_Intersects(geom_3857, {envelope})
        ),
        clipped_rows AS (
            SELECT *
            FROM tile_rows
            WHERE geom IS NOT NULL
              AND NOT ST_IsEmpty(geom)
        )
        SELECT ST_AsMVT(
            {struct_expression},
            '{spec.layer_name}',
            {extent},
            'geom',
            'id'
        ) AS tile
        FROM clipped_rows
    """
    return sql, table_name

同じビルドID、データセット、XYZタイルはファイルへ保存し、2回目以降はDuckDBで再生成しません。

data/tile_cache/dynamic/
└─ <build_id>/
   └─ a33/
      └─ 12/3622/1612.pbf

空タイルも0バイトファイルとして保存します。該当地物がない場所を何度開いても、同じ空間検索を繰り返さないためです。

OpenLayersでPMTilesとFlatGeobufを切り替える

FlatGeobufソース

flatgeobufパッケージのOpenLayers向けローダーを使い、VectorTileSourceとして扱います。

/**
 * 1つのFlatGeobufをOpenLayersのVectorTileSourceとして扱う。
 *
 * FlatGeobuf側のPacked Hilbert R-treeを使い、表示タイルの範囲に交差する
 * FeatureだけをRange Requestで取得する。URLはビルドID付きなので、公開後は
 * 内容が変わらない静的ファイルとして長期キャッシュできる。
 */
function createFlatGeobufSource(url) {
  const source = new VectorTileSource({
    projection: "EPSG:3857",
    tileUrlFunction,
    wrapX: false,
  });

  source.setTileLoadFunction(
    createTileLoadFunction(source, url, "EPSG:4326"),
  );
  return source;
}

PMTilesソース

const pmtilesLayer = new VectorTileLayer({
  ...commonLayerOptions(layerConfig, "PMTiles"),
  source: new PMTilesVectorSource({
    url: layerConfig.pmtiles.url,
  }),
});

自動切替

既定では次のように切り替えます。

if (zoom < Number(policy.pmtiles_min_zoom)) {
  // 低ズームは件数メッシュ。
  overviewActive = true;
} else if (zoom >= Number(policy.flatgeobuf_min_zoom)) {
  // 高ズームは元形状を範囲取得。
  useFlatGeobuf = true;
} else {
  // 中間は事前生成した20m LODのPMTiles。
  usePmtiles = true;
}

画面では、自動切替のほか、PMTiles固定、FlatGeobuf固定、動的MVTを選べます。同じ場所、同じズームで方式を切り替えられるため、ブラウザーのNetworkとPerformanceで比較しやすくしています。

StyleはFeatureごとに作らない

大量のFeatureを表示するときは、スタイル関数内で毎回new Style()を呼ぶのも負荷になります。

今回は区域区分ごとのStyleを先に作り、全Featureで共有しています。

/** A33の区域区分ごとに再利用するスタイル。 */
const A33_STYLES = {
  warning: new Style({
    fill: new Fill({ color: "rgba(255, 213, 0, 0.36)" }),
    stroke: new Stroke({ color: "rgba(166, 118, 0, 0.95)", width: 1.1 }),
  }),
  special: new Style({
    fill: new Fill({ color: "rgba(230, 44, 55, 0.38)" }),
    stroke: new Stroke({ color: "rgba(155, 18, 27, 0.98)", width: 1.2 }),
  }),
};

/**
 * A33/A46の表示用スタイルを返す。
 *
 * StyleオブジェクトをFeatureごとに作ると、全国データではGC負荷が目立つ。
 * ここでは少数のStyleを先に作り、全Featureで共有する。
 */
export function hazardStyle(feature) {
  const dataset = feature.get("dataset");
  const geometryKind = feature.get("geometry_kind") ?? "polygon";

  if (dataset === "A33") {
    const zoneClass = feature.get("zone_class") ?? "unknown";
    return A33_STYLES[zoneClass] ?? A33_STYLES.unknown;
  }

  return A46_POLYGON_STYLE;
}

実際のソースでは面、線、点を分けています。

詳細属性はクリック時にAPIから読む

MVTとFlatGeobufには、表示と選択に必要な属性だけを入れています。区域をクリックすると、datasetとidを使って元テーブルから1件取得します。

GET /api/features/A33/12345?build_id=<build_id>

返す内容は次のようなものです。

{
  "dataset": "A33",
  "id": 12345,
  "phenomenon_name": "土石流",
  "zone_name": "土砂災害警戒区域（指定済）",
  "area_no": "...",
  "area_name": "...",
  "address": "...",
  "notice_date": "...",
  "geometry_kind": "polygon",
  "geometry": {
    "type": "Polygon",
    "coordinates": []
  }
}

APIはリクエストのbuild_idとDuckDB内のbuild_idを確認します。画面が古いmanifestを保持している間にDBだけ新しくなった場合は、別ビルドの属性を返さずHTTP 409にします。

NginxでRange Requestとキャッシュを設定する

PMTilesとFlatGeobufはRange Requestを使うため、Nginxでは静的ファイルをそのまま配信します。

location /data/builds/ {
    alias /srv/hazard-data/builds/;

    sendfile on;
    gzip off;

    add_header Accept-Ranges bytes always;
    add_header Cache-Control "public, max-age=31536000, immutable" always;
}

ルートのmanifest.jsonは更新されるため、長期キャッシュしません。

location = /data/manifest.json {
    alias /srv/hazard-data/manifest.json;
    add_header Cache-Control "no-store" always;
}

APIは同一オリジンの/api/へプロキシします。FlatGeobufの範囲取得で余分なCORS設定を増やさないためです。

公開前に検査する

ビルド処理では、静的ファイルを公開する前に次を確認します。

A33/A46の元件数
LODの件数が元件数と一致するか
IDの重複
NULL、空ジオメトリ
未対応のGeometryCollectionなどが残っていないか
必要なR-treeが存在するか
manifestとDuckDBのビルドID
FlatGeobufとPMTilesのファイルヘッダー
公開ファイルのSHA-256
manifestに動的MVT設定が揃っているか

検査結果はJSONで残します。

data/published/reports/validation_<build_id>.json
data/published/reports/build_<build_id>.json

現在公開中のデータを再検査する場合は次を実行します。

docker compose --profile tools run --rm pipeline validate --runtime

起動後の疎通確認

Linux / WSLでは次を実行します。

bash scripts/smoke_test.sh

PowerShellでは次のとおりです。

powershell -ExecutionPolicy Bypass -File .\scripts\smoke_test.ps1

このスクリプトは次を確認します。

Web側のmanifestを取得できる
APIとmanifestのビルドIDが一致する
PMTilesまたはFlatGeobufへRange Requestを送り、HTTP 206が返る
指定した1024バイトだけ取得できる

手作業で確認する場合は次のようにします。

curl -i \
  -H "Range: bytes=0-1023" \
  http://localhost:8080/data/builds/<build_id>/pmtiles/a33_polygon.pmtiles

期待する応答は次です。

HTTP/1.1 206 Partial Content
Content-Range: bytes 0-1023/...
Accept-Ranges: bytes

静的チェックと単体テスト

Pythonの構文、docstring、JavaScript構文、前処理/API/Web間の主要な契約を確認するスクリプトを入れています。

make check

個別に実行する場合は次のとおりです。

python scripts/static_check.py
PYTHONPATH=pipeline python -m pytest -q tests
PYTHONPATH=api python -m pytest -q api/tests

static_check.pyでは、Pythonモジュールと関数にdocstringがあるかも確認します。実装を後で修正したときに、コメントのない補助関数が増えないようにしました。

今回の配布物では、静的チェック、前処理の入力探索テスト2件、API周辺の単体テスト12件まで確認しています。作成環境にはDockerとTippecanoeがなく、外部npm/PyPIレジストリも名前解決できなかったため、Docker Composeによる一括ビルド、DuckDB Spatialの実データ変換、Tippecanoe実行、Vite本番ビルドは未実行です。実際に使う環境では、最初にサンプルビルドとスモークテストを実行してください。

速度を比較する方法

この記事では、実測していない数値は書きません。サンプルデータは小さいため、全国データの性能を判断する材料にもなりません。

実データでは、場所とズームを固定し、同じ条件で比較します。

ケース	目安の表示範囲	確認する表示
z6	全国	件数メッシュ
z9	都道府県	PMTiles
z12	市区町村	PMTilesと動的MVT
z15	区域詳細	FlatGeobuf、PMTiles、動的MVT

ブラウザーでは次を記録します。

HTTPリクエスト数
転送バイト数
初回表示までの時間
2回目の表示時間
Main threadの処理時間
メモリー使用量
パン操作後の再描画時間
表示Feature数

結果表は次のようにしておくと比較しやすくなります。

方式	zoom	備考
PMTiles	12
FlatGeobuf	15
動的MVT	12	cache miss
動的MVT	12	cache hit

HTTP経路だけを繰り返し測るスクリプトも入れています。

python scripts/benchmark_http.py \
  --base-url http://localhost:8080

動的MVTは次で確認します。

python scripts/benchmark_dynamic_mvt.py \
  --base-url http://localhost:8080 \
  --dataset A33 \
  --zoom 12 \
  --repeat 5

動的MVTのレスポンスには、キャッシュ状態と参照テーブルをヘッダーへ付けています。

X-Tile-Cache: HIT / MISS
X-Tile-Source-Table: a33_lod100 / a33_lod20 / a33_hazard

主な設定値

.envで調整します。

設定	既定値	内容
`OVERVIEW_CELL_SIZE_M`	20000	件数メッシュの一辺
`LOD_COARSE_TOLERANCE_M`	100	動的MVT低ズーム用LOD
`LOD_DETAIL_TOLERANCE_M`	20	PMTiles・動的MVT中ズーム用LOD
`OVERVIEW_MAX_ZOOM`	9	件数メッシュの作成上限
`PMTILES_MIN_ZOOM`	9	個別区域PMTilesの開始ズーム
`PMTILES_MAX_ZOOM`	14	PMTilesの最大ズーム
`FLATGEOBUF_MIN_ZOOM`	14	FlatGeobufへ切り替えるズーム
`DYNAMIC_MVT_MAX_ZOOM`	16	動的MVTの最大ズーム
`PUBLISH_RETAIN_BUILDS`	3	静的ファイルとキャッシュの保持世代

全国データで最初に調整するのは、PMTILES_MIN_ZOOMとLOD_DETAIL_TOLERANCE_Mです。

タイルが大きすぎる場合：PMTilesの開始ズームを上げる
形状が細かすぎる場合：20m LODの値を少し大きくする
小さな区域が元形状へ戻りすぎる場合：simplification_fallback件数を確認する
低ズームで分布が粗すぎる場合：件数メッシュのセル寸法を小さくする

一度に値を大きく変えず、ファイルサイズと画面表示を見ながら調整します。

運用上の注意

PMTilesは更新時に作り直す

PMTilesは読み取り向けの単一アーカイブです。区域が更新されたときは、該当タイルだけを直接書き換えるのではなく、ビルドを実行して新しいPMTilesを作ります。

ビルドID付きの別URLで公開するため、古いファイルをすぐ消す必要はありません。既定では3世代を残します。

FlatGeobufを低ズームから表示しない

FlatGeobufは範囲検索できますが、全国範囲は全国データを意味します。高ズーム用という役割を崩さない方が安定します。

動的MVTは条件検索が必要な場面へ使う

たとえば次のような条件は、静的PMTilesだけで増やし続けるより動的MVTが向いています。

特定の告示日以降
市区町村指定
現象種別指定
利用者ごとの権限
更新直後の暫定表示

通常表示を静的ファイルへ任せ、条件が必要な場面だけ動的処理を使う方が、構成を整理しやすくなります。

公示図書の代わりにはしない

国土数値情報のA33/A46は、Web表示や分析には便利ですが、区域の厳密な位置を確定する資料ではありません。重要事項説明、法的判断、境界付近の確認では、各都道府県が公開する告示図書や原典資料を確認します。

今回の低ズーム件数メッシュは、さらに概観用の表示です。セルが塗られていても、セル全体が区域という意味ではありません。

まとめ

FlatGeobufを使うこと自体は難しくありません。ただし、全国規模の区域をすべてFlatGeobufへ置き換えるだけでは、低ズームの描画負荷が残ります。

今回は次の形にしました。

広域の概観        件数メッシュPMTiles
通常の区域表示    20m LOD PMTiles
詳細な元形状      FlatGeobuf
検索・属性取得    FastAPI + DuckDB Spatial
比較・条件表示    DuckDB動的MVT

効果が大きいのは、ファイル形式の変更だけではなく、縮尺ごとに表示内容を変えたことです。低ズームでは件数分布を見せ、高ズームになってから元形状を読むようにしました。

この構成なら、通常の地図操作は静的配信で軽くしつつ、DuckDBを検索、詳細属性、動的抽出へ残せます。前の記事の実装を捨てずに、負荷のかかる部分だけを静的ファイルへ移す形です。

参考資料

国土数値情報 A33 土砂災害警戒区域データ第2.0版
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A33-v2_0.html
国土数値情報 A46 地すべり防止区域データ第1.1版
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A46-v1_1.html
DuckDB Spatial: R-tree Indexes
https://duckdb.org/docs/stable/core_extensions/spatial/r-tree_indexes.html
DuckDB Spatial Functions
https://duckdb.org/docs/stable/core_extensions/spatial/functions.html
DuckDB Spatial: GDAL Integration
https://duckdb.org/docs/stable/core_extensions/spatial/gdal.html
FlatGeobuf
https://flatgeobuf.org/
GDAL FlatGeobuf driver
https://gdal.org/en/stable/drivers/vector/flatgeobuf.html
Tippecanoe
https://github.com/felt/tippecanoe
PMTiles
https://docs.protomaps.com/pmtiles/
PMTiles + OpenLayers
https://docs.protomaps.com/pmtiles/openlayers

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

国土数値情報 土砂災害警戒区域などデータを PMTiles + FlatGeobuf + DuckDB で軽く表示してみた

国土数値情報 土砂災害警戒区域などデータを PMTiles + FlatGeobuf + DuckDB で軽く表示してみた

今回作ったもの

FlatGeobufだけにしなかった理由

A33とA46の属性で気をつけたところ

フォルダー構成

主なバージョン

Docker Composeの分け方

サンプルデータで起動する

Linux / WSL / Ubuntu

Windows PowerShell

国土数値情報のA33/A46を取得する

最初に取得対象だけ確認する

利用条件を確認してから取得する

中断したファイルは.partから再開する

ダウンロード後にZIPを検査する

SHA-256と取得元を残す

OS別のラッパー

都道府県コード順に一括処理する

全国版と全都道府県

A46の過年度をどう扱うか

再取得、ドライラン、オフライン確認

ダウンロード後に前処理を実行する

高速化1：元形状と表示用LODを分ける

高速化2：低ズームは件数メッシュへ切り替える

高速化3：元形状とLODへR-treeを作る

FlatGeobufは元形状で出力する

PMTilesは20m LODから作る

公開ファイルをビルドIDで揃える

動的MVTは比較用として残す

OpenLayersでPMTilesとFlatGeobufを切り替える

FlatGeobufソース

PMTilesソース

自動切替

StyleはFeatureごとに作らない

詳細属性はクリック時にAPIから読む

NginxでRange Requestとキャッシュを設定する

公開前に検査する

起動後の疎通確認

静的チェックと単体テスト

速度を比較する方法

主な設定値

運用上の注意

PMTilesは更新時に作り直す

FlatGeobufを低ズームから表示しない

動的MVTは条件検索が必要な場面へ使う

公示図書の代わりにはしない

まとめ

参考資料

国土数値情報土砂災害警戒区域などデータを PMTiles + FlatGeobuf + DuckDB で軽く表示してみた

国土数値情報土砂災害警戒区域などデータを PMTiles + FlatGeobuf + DuckDB で軽く表示してみた

中断したファイルは`.part`から再開する