ふと思ったけど、日本全国250ｍメッシュだとGISフォーマットごとの容量としてはどのくらいなるの？WEB表示できる？
普段GISデータを触っていると、「この範囲を全部メッシュ化したら、ファイルはどのくらい大きくなるのだろう」と気になることがあります。GeoJsonだと無理っぽいし、shapefileは使いたくないし。悩む部分です。
今回は、日本全国の陸域を250mメッシュのポリゴンで持たせた場合を考えてみました。属性は10項目で、すべて文字列型、1項目あたり最大50文字とします。
単に容量を比べるだけでは実際の使い方が見えにくいので、後半ではOpenLayers、PMTiles、FastAPI、DuckDB Spatialを使ってWeb地図へ表示する構成も整理します。
今回の条件
最初に条件をそろえておきます。
| 項目 | 条件 |
|---|---|
| 対象範囲 | 日本全国の陸域 |
| メッシュ数 | 約600万件 |
| 形状 | 250mメッシュ相当の長方形ポリゴン |
| 座標参照系 | EPSG:4326 |
| 属性数 | 文字列10項目 |
| 文字数 | 1項目あたり最大50文字 |
| ベンチマーク時の文字数 | 上限側を見るため、全件・全項目を50文字で作成 |
| 空間インデックス | GeoPackage、FlatGeobuf、Shapefileで作成 |
| 容量表記 | 1GB＝1,000,000,000バイト |
ここでの「最大50文字」は、平均50文字という意味ではありません。
実際の値は0～50文字の範囲に収まる想定です。ただし、容量の上限側をつかむため、簡易ベンチマークでは10項目すべてを50文字まで埋めています。そのため、実データの平均文字数が短ければ、可変長形式の容量は表より小さくなります。
250mメッシュは何件くらいになるのか
国土地理院が公表している2026年4月1日時点の国土面積は、377,974.87km²です。
250m四方を単純な正方形として計算すると、1セルの面積は0.0625km²になります。
0.25km × 0.25km = 0.0625km²
国土面積を割ると、約604.8万セルです。
377,974.87km² ÷ 0.0625km²
≒ 6,047,598セル
そこで、この記事では分かりやすく600万件を基準にします。
ただし、標準地域メッシュの4分の1地域メッシュは、厳密な250m正方形ではありません。総務省統計局の定義では、緯度方向7.5秒、経度方向11.25秒の区画です。緯度によって東西方向の実距離が変わります。
陸域の抽出方法でも件数は変わります。
- メッシュ中心が陸地にあるセルだけを残す
- 陸域ポリゴンと少しでも交差するセルを残す
- 湖沼や河川を陸域側へ含める
- 統計値が入っているセルだけを残す
したがって、600万件は容量比較のための概算です。本番データを作るときは、使用する陸域マスクで最終件数を数えます。
文字列だけで何GBになるのか
全600万件について、10項目すべてが最大の50文字まで入った場合、文字数の合計は次のとおりです。
6,000,000件 × 10項目 × 50文字
= 3,000,000,000文字
30億文字です。
半角英数字が中心で、UTF-8上で1文字1バイトと考えられる場合、文字列本体だけで約3GBになります。
3,000,000,000文字 × 1バイト
= 3.0GB
日本語が中心で、1文字をおおむね3バイトとして考えると、文字列本体だけで約9GBです。
3,000,000,000文字 × 3バイト
= 9.0GB
同じ50文字でも、半角英数字と日本語ではバイト数が大きく変わります。
実際には、メッシュコードや年月のような半角文字、自治体名や説明文のような日本語が混在するはずです。その場合は、3GBと9GBの間に収まると考えるのが自然です。
ポリゴン形状はどのくらいか
単純な長方形ポリゴンは、始点へ戻る座標を含めて5点を持ちます。
一般的な2次元PolygonをWKBで保存すると、1件あたり約93バイトです。
93バイト × 6,000,000件
= 558,000,000バイト
≒ 0.56GB
実際のファイルには、地物ごとの管理情報、属性のオフセット、ファイルヘッダー、空間インデックスなどが加わります。
それでも今回の条件では、ポリゴン形状より10項目の文字列属性の方が容量へ大きく効きます。
簡易ベンチマークの結果
全国600万件を直接作ると時間とディスクをかなり使うため、1万件の規則格子を各形式へ保存し、ファイルサイズを600倍しました。
6,000,000件 ÷ 10,000件 = 600倍
文字列は、行ごと、列ごとに内容が変わるようにしています。自治体名や区分名のような繰り返しの多い実データより、圧縮が効きにくい条件です。
半角英数字50文字を10項目へ入れた場合
| フォーマット | 全国600万件への換算 | 備考 |
|---|---|---|
| Shapefile一式 | 約4.1GB | DBFを最適幅へ縮小し、QIX空間インデックスを含む |
| FlatGeobuf | 約4.4GB | 空間インデックスを含む |
| GeoPackage | 約4.6GB | R-tree空間インデックスを含む |
| GeoJSON | 約5.4GB | 非圧縮、空間インデックスなし |
| GeoParquet＋ZSTD | 約2～4GBを想定 | 文字列の重複率と並び順で変動 |
日本語50文字を10項目へ入れた場合
| フォーマット | 全国600万件への換算 | 備考 |
|---|---|---|
| Shapefile一式 | 約10.1GB | UTF-8では50文字がおおむね150バイトになる |
| FlatGeobuf | 約10.4GB | 空間インデックスを含む |
| GeoJSON | 約11.4GB | 非圧縮、空間インデックスなし |
| GeoPackage | 約12.8GB | R-tree空間インデックスを含む |
| GeoParquet＋ZSTD | 約5～9GBを想定 | 同じ値が多ければさらに小さくなる |
ここで意外だったのは、半角英数字のShapefileが比較的小さいことです。ただし、これは文字列フィールドを実データに合わせて縮小した結果です。全国600万件を1ファイルへ入れる運用を勧めているわけではありません。
この数値は上限側の試算
今回のベンチマークでは、全600万件について、10項目すべてを50文字で埋めた状態を想定しています。
実データでは、例えば次のような列が混ざります。
mesh_code : 10文字
prefecture_name : 数文字
city_name : 数文字から十数文字
category : 数文字
source_name : 数文字から数十文字
description : 0～50文字
このように、多くの列が50文字未満であれば、GeoPackage、FlatGeobuf、GeoJSON、GeoParquetは今回の表より小さくなります。
一方、ShapefileのDBFは固定幅です。値が短くても、設定したフィールド幅に基づいてレコード領域を確保します。Shapefileだけは「実際に入っている文字数」と「DBFへ設定した幅」を分けて考える必要があります。
ShapefileはDBFの幅に注意する
Shapefileの属性はDBFへ保存されます。
GDALのShapefileドライバーでは、幅を指定していない文字列フィールドは既定で80文字として作られます。半角50文字の列をそのまま10列出力すると、必要以上の余白が残ります。
今回の測定では、次のレイヤー作成オプションを指定しました。
RESIZE=YES
ENCODING=UTF-8
SPATIAL_INDEX=YES
RESIZE=YESを付けると、DBFの文字列幅を実データに合わせて縮められます。
半角50文字の場合、RESIZE=YESを指定した全国換算は約4.1GBでした。指定しない場合は約5.7GBまで増えました。
日本語50文字はUTF-8でおおむね150バイトになるため、DBFだけで約9GBになります。
さらに、GDALのドキュメントでは、他ソフトとの互換性を考え、.shpと.dbfをそれぞれ2GB以下に抑えることが推奨されています。今回の規模では全国1ファイルを避け、都道府県別や1次メッシュ別に分割する必要があります。
容量だけを見るとShapefileも健闘しますが、全国規模の保管形式としては選びにくい、というのが実際のところです。
GeoParquetは文字列の中身でかなり変わる
GeoParquetは、分析用の原本として使いやすい形式です。
Parquetは列指向なので、列ごとに圧縮されます。同じ値が何度も現れる文字列列では、辞書エンコーディングや圧縮が効きます。
圧縮しやすい列の例です。
都道府県名
市区町村名
土地利用区分
データ種別
出典名
更新年月
反対に、次のような列は圧縮しにくくなります。
メッシュごとに異なる説明文
長いURL
ランダムな識別文字列
暗号化済み文字列
今回の最大50文字条件では、GeoParquet＋ZSTDの容量は次の範囲を見ておくとよさそうです。
| 属性の内容 | 全国600万件の目安 |
|---|---|
| 区分名や自治体名など、繰り返しが多い | 約1～3GB |
| 半角中心で、各行の内容がかなり異なる | 約2～4GB |
| 日本語中心で、各行の内容がかなり異なる | 約5～9GB |
GeoParquetは文字列の中身で差が出やすいため、全国版を作る前に、1県分または1次メッシュ数個分で実測するのが確実です。
フォーマットごとの使い分け
GeoParquet
分析用の原本として使います。
- 10項目すべてを保持する
- ZSTDで圧縮する
- 1次メッシュや地域ブロックで分割する
- DuckDBから直接読む
600万件を1ファイルにまとめることもできますが、更新、差し替え、並列処理を考えると、分割データセットの方が扱いやすくなります。
data/master/mesh_250m/
├── first_mesh=3036/
├── first_mesh=3037/
├── first_mesh=3622/
├── first_mesh=3623/
└── ...
DuckDB Spatial
今回の構成では、DuckDB Spatialを正式にデータ基盤へ組み込みます。
DuckDBはGISファイル形式そのものではありませんが、GeoParquetなどの原本を読み込み、検索、集計、空間抽出、Web API向けのデータ作成を行う処理エンジンとして使いやすいです。Spatial拡張を読み込むと、GEOMETRY型、空間関数、R-tree空間インデックス、GDAL経由のGISファイル入出力を利用できます。
この記事では、役割を次のように分けます。
GeoParquet : 保管・交換用
DuckDB : 検索・集計・空間解析・Web API用
PMTiles : Web地図の描画用
全国250mメッシュをDuckDBへ取り込む場合は、例えば次のようなテーブルを用意します。
SQL全文を表示する
-- Spatial拡張を一度インストールする。
-- Dockerイメージのビルド時や初期構築時に実行しておくと、
-- 実運用時に外部ネットワークへ接続せずにLOADできる。
INSTALL spatial;
LOAD spatial;
-- 全国250mメッシュ本体。
-- geometryと10項目の文字列属性を同じテーブルへ保持する例。
CREATE TABLE mesh_250m (
mesh_code VARCHAR PRIMARY KEY,
attr01 VARCHAR,
attr02 VARCHAR,
attr03 VARCHAR,
attr04 VARCHAR,
attr05 VARCHAR,
attr06 VARCHAR,
attr07 VARCHAR,
attr08 VARCHAR,
attr09 VARCHAR,
attr10 VARCHAR,
geometry GEOMETRY
);
-- BBOX検索や点選択など、geometryとの交差判定を多く使う場合は、
-- geometry列へR-tree空間インデックスを作成する。
CREATE INDEX mesh_250m_rtree
ON mesh_250m USING RTREE (geometry);
DuckDBは単一の.duckdbファイルとして保持できるため、Docker環境でも扱いやすいです。ただし、データ更新のたびに複数ユーザーが同時編集するような業務システムを想定するのであれば、後述するPostGISを検討します。
GeoPackage
QGISで確認したり、一部を手作業で修正したりする用途に向いています。
全国1ファイルではなく、都道府県別、地方別、1次メッシュ別へ分割した方が安全です。
FlatGeobuf
単一ファイルで空間インデックスを持てるため、BBOXを指定した部分取得に向いています。
ただし、全国表示をFlatGeobufだけで行うより、Web表示用にはMVTへ変換した方が軽くなります。
Shapefile
受け渡し先から指定された場合に使う形式だと思います。
- DBFの固定幅で容量が増えやすい
- 2GBを超えるファイルは互換性に不安がある
- フィールド名に制約がある
- 複数ファイルを一式で管理する必要がある
- prjファイルがない場合が多い
全国分を出力する場合は、属性を減らし、地域ごとに分割します。
GeoJSON
小範囲のAPI応答やデバッグには便利です。
一方、全国600万件の一括配信には向きません。ファイルが大きいだけでなく、ブラウザ側で大量のFeatureをJavaScriptオブジェクトへ展開する負荷もかかります。
測定用コード
測定環境は次のとおりです。
Python : 3.13
GeoPandas : 1.1.2
Pyogrio : 0.12.1
Shapely : 2.1.2
GDAL : 3.11.4
CRS : EPSG:4326
必要なライブラリを入れます。
# GeoPandasは表形式データとgeometryをまとめて扱うために使用する。
# PyogrioはGDAL経由で各GIS形式へ書き出すために使用する。
# Shapelyは測定用の長方形ポリゴンを作るために使用する。
python -m pip install geopandas pyogrio shapely
以下が測定用スクリプトです。コードだけを見ても条件が分かるように、コメントを多めに入れています。
長いソースは記事を読み進めやすくするため折りたたんでいます。コード全文を表示するをクリックすると展開できます。
Pythonコード全文を表示する
"""全国陸域の250mメッシュを想定し、GIS形式ごとの容量を比較する。
全国約600万件を直接生成すると、測定だけで大きなディスク容量と時間を使う。
そこで、このスクリプトでは1万件の規則格子を作成して各形式へ保存し、
得られたファイルサイズを全国600万件へ線形換算する。
属性条件
----------
- 文字列属性は10項目
- 1項目あたりの最大文字数は50文字
- 上限側の容量を見るため、測定データは全項目をちょうど50文字で作る
- 半角英数字と日本語の2種類を切り替えて測定できる
注意
----
線形換算は概算である。全国規模になると、空間インデックスの階層、
SQLiteページ、ファイルシステム、地物の並び順などの影響が加わる。
最終判断では、実データの一部を使った追加測定を行うこと。
"""
from __future__ import annotations
import argparse
import base64
import hashlib
import shutil
from pathlib import Path
from typing import Literal
import geopandas as gpd
from shapely.geometry import box
# 測定結果の出力先。
# 実行のたびに中身を削除して作り直すため、必要なファイルは置かないこと。
OUTPUT_DIR = Path("data/benchmark_250m_max50")
# 全国600万件を直接作らず、まず1万件だけ保存してサイズを測る。
SAMPLE_FEATURE_COUNT = 10_000
# 全国の陸域250mメッシュ数を約600万件として換算する。
TARGET_FEATURE_COUNT = 6_000_000
# 文字列属性は10項目とする。
ATTRIBUTE_COUNT = 10
# 今回の条件は「平均50文字」ではなく「最大50文字」。
# ベンチマークでは上限側を見るため、すべての値を50文字で生成する。
MAX_TEXT_LENGTH = 50
# 4分の1地域メッシュに相当する緯度・経度方向の幅。
# 緯度方向は7.5秒、経度方向は11.25秒である。
MESH_WIDTH_DEGREES = 11.25 / 3600.0
MESH_HEIGHT_DEGREES = 7.5 / 3600.0
# Shapefileは複数ファイルで1レイヤを構成する。
# 容量比較では、形状、属性、座標系、文字コード、空間インデックスを合算する。
SHAPEFILE_COMPONENTS = {
".shp", # ポリゴン形状
".shx", # 形状インデックス
".dbf", # 属性テーブル
".prj", # 座標参照系
".cpg", # 文字コード
".qix", # GDALが作成する空間インデックス
}
# 日本語モードで使用する文字集合。
# 同じ文章を全行へ入れると圧縮が効きすぎるため、ハッシュ値を使って
# 行ごと・列ごとに異なる文字列を組み立てる。
JAPANESE_CHARACTERS = (
"あいうえおかきくけこさしすせそたちつてとなにぬねの"
"はひふへほまみむめもやゆよらりるれろわをん"
"日本全国陸域格子属性情報防災人口地形雨量道路河川"
"土地利用更新年月出典名称区分説明"
)
# コマンドライン引数で選べる文字列モードを型として定義する。
TextMode = Literal["ascii", "japanese"]
def parse_arguments() -> argparse.Namespace:
"""コマンドライン引数を読み込む。
Returns
-------
argparse.Namespace
text_mode属性に"ascii"または"japanese"を持つ名前空間。
"""
parser = argparse.ArgumentParser(
description="全国250mメッシュ相当のGISファイル容量を概算します。"
)
parser.add_argument(
"--text-mode",
choices=("ascii", "japanese"),
default="ascii",
help="属性文字列の種類。既定値はasciiです。",
)
return parser.parse_args()
def make_ascii_text(row_index: int, column_index: int) -> str:
"""行と列ごとに異なる50文字の半角文字列を作る。
全レコードへ同じ文字列を入れると、圧縮可能な形式で結果が
不自然に小さくなる。そこで、行番号と列番号からSHAKE-256の
ダイジェストを作り、URL-safe Base64へ変換して使用する。
Parameters
----------
row_index : int
地物の0始まり行番号。
column_index : int
属性列の1始まり列番号。
Returns
-------
str
ちょうど50文字の半角文字列。
"""
seed = f"{column_index}:{row_index}".encode("utf-8")
# Base64変換後に50文字以上を確保できるよう、48バイト生成する。
digest = hashlib.shake_256(seed).digest(48)
encoded = base64.urlsafe_b64encode(digest).decode("ascii")
# 今回は最大容量側を見るため、常に上限の50文字を返す。
return encoded[:MAX_TEXT_LENGTH]
def make_japanese_text(row_index: int, column_index: int) -> str:
"""行と列ごとに異なる50文字の日本語文字列を作る。
SHA-256の各バイトを日本語文字集合の位置へ割り当てる。
必要な50文字へ達するまでハッシュ生成を繰り返す。
Parameters
----------
row_index : int
地物の0始まり行番号。
column_index : int
属性列の1始まり列番号。
Returns
-------
str
ちょうど50文字の日本語文字列。
"""
characters: list[str] = []
counter = 0
while len(characters) < MAX_TEXT_LENGTH:
seed = f"{column_index}:{row_index}:{counter}".encode("utf-8")
digest = hashlib.sha256(seed).digest()
# 0～255の各バイト値を、文字集合の範囲へ剰余で割り当てる。
characters.extend(
JAPANESE_CHARACTERS[value % len(JAPANESE_CHARACTERS)]
for value in digest
)
counter += 1
return "".join(characters[:MAX_TEXT_LENGTH])
def make_text(
row_index: int,
column_index: int,
text_mode: TextMode,
) -> str:
"""指定された文字種で測定用文字列を作る。"""
if text_mode == "ascii":
return make_ascii_text(row_index, column_index)
return make_japanese_text(row_index, column_index)
def build_grid(text_mode: TextMode) -> gpd.GeoDataFrame:
"""規則格子ポリゴンと10個の文字列属性を作成する。
Parameters
----------
text_mode : TextMode
"ascii"または"japanese"。
Returns
-------
geopandas.GeoDataFrame
EPSG:4326の長方形ポリゴンと文字列属性10項目を持つデータ。
"""
# 1行200セルとし、1万件なら50行の規則格子になる。
grid_column_count = 200
geometries = []
for feature_index in range(SAMPLE_FEATURE_COUNT):
row_index, column_index = divmod(feature_index, grid_column_count)
# 測定用なので、東京付近を起点とする単純な規則格子を作る。
min_x = 139.0 + column_index * MESH_WIDTH_DEGREES
min_y = 35.0 + row_index * MESH_HEIGHT_DEGREES
max_x = min_x + MESH_WIDTH_DEGREES
max_y = min_y + MESH_HEIGHT_DEGREES
geometries.append(box(min_x, min_y, max_x, max_y))
# attr01～attr10を作り、各値を上限の50文字で埋める。
attributes = {
f"attr{attribute_index:02d}": [
make_text(row_index, attribute_index, text_mode)
for row_index in range(SAMPLE_FEATURE_COUNT)
]
for attribute_index in range(1, ATTRIBUTE_COUNT + 1)
}
grid = gpd.GeoDataFrame(
attributes,
geometry=geometries,
crs="EPSG:4326",
)
# 条件違いのデータを誤って測らないよう、文字数を明示的に検査する。
for column_name in attributes:
maximum_length = grid[column_name].str.len().max()
if maximum_length != MAX_TEXT_LENGTH:
raise ValueError(
f"{column_name}の最大文字数が条件と一致しません: "
f"{maximum_length}"
)
return grid
def shapefile_size(shapefile_path: Path) -> int:
"""Shapefileを構成する関連ファイルの合計バイト数を返す。"""
return sum(
candidate.stat().st_size
for candidate in shapefile_path.parent.iterdir()
if candidate.stem == shapefile_path.stem
and candidate.suffix.lower() in SHAPEFILE_COMPONENTS
)
def output_size(output_path: Path) -> int:
"""単一ファイル形式またはShapefile一式の容量を返す。"""
if output_path.suffix.lower() == ".shp":
return shapefile_size(output_path)
return output_path.stat().st_size
def write_outputs(grid: gpd.GeoDataFrame) -> dict[str, Path]:
"""同じGeoDataFrameを4種類のGIS形式へ保存する。
Parameters
----------
grid : geopandas.GeoDataFrame
測定対象のポリゴンデータ。
Returns
-------
dict[str, pathlib.Path]
フォーマット名と出力先パスの対応。
"""
outputs = {
"FlatGeobuf": OUTPUT_DIR / "mesh_flatgeobuf.fgb",
"GeoPackage": OUTPUT_DIR / "mesh_geopackage.gpkg",
"Shapefile": OUTPUT_DIR / "mesh_shapefile.shp",
"GeoJSON": OUTPUT_DIR / "mesh_geojson.geojson",
}
# FlatGeobufは単一ファイル内に空間インデックスを持つ。
grid.to_file(
outputs["FlatGeobuf"],
driver="FlatGeobuf",
engine="pyogrio",
index=False,
)
# GeoPackageはR-tree空間インデックスを明示的に作成する。
grid.to_file(
outputs["GeoPackage"],
layer="mesh",
driver="GPKG",
engine="pyogrio",
index=False,
layer_options={"SPATIAL_INDEX": "YES"},
)
# Shapefileは次の3点を明示する。
# 1. UTF-8で属性を書き出す。
# 2. RESIZE=YESでDBFの既定幅80文字を実データ幅へ縮める。
# 3. QIX空間インデックスを作り、他形式と条件をそろえる。
grid.to_file(
outputs["Shapefile"],
driver="ESRI Shapefile",
engine="pyogrio",
index=False,
layer_options={
"ENCODING": "UTF-8",
"RESIZE": "YES",
"SPATIAL_INDEX": "YES",
},
)
# GeoJSONは非圧縮で保存する。
# 全国一括配信用ではなく、容量比較の基準として出力する。
grid.to_file(
outputs["GeoJSON"],
driver="GeoJSON",
engine="pyogrio",
index=False,
)
return outputs
def print_estimates(outputs: dict[str, Path]) -> None:
"""実測値と全国600万件への線形換算値を表示する。"""
scale = TARGET_FEATURE_COUNT / SAMPLE_FEATURE_COUNT
print("format sample_MB national_GB")
print("--------------- ---------- -----------")
for format_name, output_path in outputs.items():
measured_bytes = output_size(output_path)
estimated_bytes = measured_bytes * scale
print(
f"{format_name:15s} "
f"{measured_bytes / 1_000_000:10.2f} "
f"{estimated_bytes / 1_000_000_000:11.2f}"
)
def main() -> None:
"""測定データを作成し、各形式の容量を全国規模へ換算する。"""
arguments = parse_arguments()
text_mode: TextMode = arguments.text_mode
# 前回の測定結果が混ざらないよう、出力先を作り直す。
if OUTPUT_DIR.exists():
shutil.rmtree(OUTPUT_DIR)
OUTPUT_DIR.mkdir(parents=True)
grid = build_grid(text_mode)
outputs = write_outputs(grid)
print_estimates(outputs)
if __name__ == "__main__":
main()
半角英数字を測る場合です。
# 10項目すべてへ半角英数字50文字を入れて測定する。
python benchmark_formats.py --text-mode ascii
結果は次のようになりました。
format sample_MB national_GB
--------------- ---------- -----------
FlatGeobuf 7.39 4.43
GeoPackage 7.73 4.64
Shapefile 6.76 4.06
GeoJSON 8.95 5.37
日本語を測る場合です。
# 10項目すべてへ日本語50文字を入れて測定する。
python benchmark_formats.py --text-mode japanese
format sample_MB national_GB
--------------- ---------- -----------
FlatGeobuf 17.39 10.43
GeoPackage 21.41 12.84
Shapefile 16.76 10.06
GeoJSON 18.95 11.37
この方法は単純な線形換算です。全国分を実際に作ると、空間インデックスの階層やデータの並び順で差が出ます。容量計画では、表の値へ20～30％程度の余裕を持たせておく方が安全です。
Web表示では「原本・処理・描画」を分ける
全国250mメッシュをWeb地図へ表示する場合、600万ポリゴンと最大50文字の文字列10項目を、そのままブラウザへ送る構成にはしません。そんなことをしたら、固まってしまいます。
今回の構成では、役割を3段階に分けます。
もう少し詳しく書くと
ここでは、PostGISは考えません。
ここで大事なのは、PMTilesへ10項目すべてを入れないことです。
地図描画に必要なのは、例えば次の程度です。
feature_id : MVT内の数値ID
mesh_code : 250mメッシュコード
class : 色分け用の区分値
value : 凡例や簡易ポップアップ用の値
attr01～attr10の長い文字列は、利用者がメッシュをクリックしたときにFastAPIから1件だけ取得します。
DuckDBへGeoParquetを取り込む
DuckDB SpatialはGeoParquetをデータ基盤へ取り込む処理にも使えます。
次は、GeoParquetを読み込み、DuckDBのテーブルへ保存する例です。
SQL全文を表示する
-- Spatial拡張を読み込む。
-- INSTALLは初回だけでよく、運用コンテナでは事前にインストール済みにしておく。
INSTALL spatial;
LOAD spatial;
-- GeoParquetから全国250mメッシュを読み込む。
-- Spatial拡張が有効な環境では、GeoParquetのgeometry列を
-- GEOMETRY型として扱える構成にできる。
CREATE OR REPLACE TABLE mesh_250m AS
SELECT
mesh_code,
attr01,
attr02,
attr03,
attr04,
attr05,
attr06,
attr07,
attr08,
attr09,
attr10,
geometry
FROM read_parquet('data/master/mesh_250m/**/*.parquet');
-- 全国一括検索を毎回フルスキャンしないよう、geometry列へR-treeを作る。
-- R-treeはST_Intersectsなどの空間述語を使った絞り込みで利用される。
CREATE INDEX IF NOT EXISTS mesh_250m_rtree
ON mesh_250m USING RTREE (geometry);
-- mesh_codeによる詳細属性検索も頻繁に行うため、
-- 一意性を保証できる場合は通常のインデックスや主キーも検討する。
CREATE UNIQUE INDEX IF NOT EXISTS mesh_250m_code_idx
ON mesh_250m (mesh_code);
データを毎回DuckDBへコピーしたくない場合は、GeoParquetをDuckDBから直接問い合わせる方法もあります。
-- 原本をParquetのまま保持し、必要な列だけを直接読む例。
-- Parquetは列指向なので、SELECTした列だけを読む運用と相性がよい。
SELECT
mesh_code,
attr01,
attr02
FROM read_parquet('data/master/mesh_250m/**/*.parquet')
WHERE attr01 = '対象区分';
ただし、Web APIで繰り返しBBOX検索するのであれば、GEOMETRY列を持つDuckDBテーブルへ取り込み、R-treeを作った方が扱いやすい場面があります。原本はGeoParquet、実行用はDuckDB、と分けておくと構成が分かりやすくなります。
DuckDB SpatialでBBOX検索する
Web地図では、現在表示している範囲だけを検索する処理がよくあります。
SQL全文を表示する
-- :min_lon、:min_lat、:max_lon、:max_latには、
-- OpenLayersから送られてきた表示範囲を渡す想定。
-- EPSG:4326でgeometryを保存している場合の例。
SELECT
mesh_code,
attr01,
attr02,
geometry
FROM mesh_250m
WHERE ST_Intersects(
geometry,
ST_MakeEnvelope(
:min_lon,
:min_lat,
:max_lon,
:max_lat
)
);
R-treeを利用するときは条件があります。DuckDBの公式ドキュメントでは、GEOMETRY型に対するR-treeと、ST_Intersects、ST_Within、ST_Containsなどの空間述語によるインデックススキャンが説明されています。
全国600万件を対象にする場合でも、毎回全件をブラウザへ返すのではなく、地図表示範囲、ズーム、属性条件でできるだけ先に絞ることが重要です。
全国表示で250mポリゴンをそのまま描かない
全国が見える縮尺で600万個の250mポリゴンを描いても、1セルはほとんど見えません。
ズームに応じて表示単位を変えます。
| ズーム | 表示するデータの例 |
|---|---|
| z0～6 | 都道府県、市区町村、20km集計 |
| z7～9 | 5km集計 |
| z10～11 | 1km集計 |
| z12以上 | 250mメッシュ |
低ズームでは、平均値、最大値、件数、代表区分などへ集約します。
例えばDuckDB側で1km集計テーブルを事前に作っておけば、Web表示時の処理を軽くできます。
-- ここではmesh_1km_codeがあらかじめ付与されている想定。
-- 実運用では、250mメッシュコードから1kmメッシュコードを生成してもよい。
CREATE OR REPLACE TABLE mesh_1km_summary AS
SELECT
mesh_1km_code,
-- 250mセル数を保持しておくと、欠損セルの確認にも使える。
COUNT(*) AS mesh_count,
-- 数値列がある場合はAVGやMAXなどを計算する。
-- 文字列列は用途に応じて代表値や件数へ置き換える。
MAX(class) AS representative_class
FROM mesh_250m
GROUP BY mesh_1km_code;
250mメッシュは十分に拡大したときだけ表示し、現在の画面に必要なタイルだけを取得します。
Web描画はPMTilesを基本にする
全国250mメッシュの基本表示は、DuckDBから毎回ポリゴンを返すより、あらかじめベクタータイルを作ってPMTilesとして配信する方が扱いやすいです。
PMTilesには、地図を描くための短い属性だけを持たせます。
mesh_code
class
value
最大50文字の10項目をすべて各タイルへ入れると、タイル容量と通信量が増えます。詳細属性はDuckDBへ残し、クリック時にAPIから取得します。
DuckDBだけでMVTを完結させない
DuckDB SpatialにはST_AsMVTGeomがあり、geometryをMVT向けのタイル座標へ変換・クリップする処理に利用できます。
一方、この記事ではDuckDBだけでMVTファイル生成まで完結することを前提にしません。
静的配信ならTippecanoeなどのタイル生成ツールへ渡してPMTilesを作る方が構成を分離しやすく、運用も分かりやすくなります。
動的MVTが必要な場合は、次のように役割を分けます。
OpenLayers
↓ Z/X/Y
FastAPI
↓ タイルのBBOXを計算
DuckDB Spatial
↓ BBOX内の地物だけ取得
FastAPI側のMVTエンコーダ
↓ PBF
OpenLayers
更新頻度の低い全国250mメッシュでは、まずPMTilesを基本にし、動的MVTは本当に必要になった段階で追加する方が無理がありません。
FastAPI＋DuckDB Spatialで詳細属性を返す
メッシュをクリックしたときは、MVTではなくDuckDBから詳細属性を1件取得します。
次の例では、mesh_codeを受け取り、最大50文字の10属性をJSONで返します。
Pythonコード全文を表示する
"""全国250mメッシュの詳細属性を返すFastAPIアプリ。
地図描画そのものはPMTilesへ任せ、FastAPIはクリックされたメッシュの
詳細属性だけをDuckDBから取得する。
この分離により、地図移動のたびに10個の文字列属性を大量送信せずに済む。
"""
from __future__ import annotations
from contextlib import closing
from pathlib import Path
from typing import Annotated
import duckdb
from fastapi import FastAPI, HTTPException, Path as ApiPath
# 実行用DuckDBファイル。
# 元のGeoParquetとは分けておくと、再構築や差し替えがしやすい。
DATABASE_PATH = Path("data/runtime/japan_250m_mesh.duckdb")
app = FastAPI(
title="Japan 250m Mesh API",
version="1.0.0",
)
DETAIL_SQL = """
SELECT
mesh_code,
attr01,
attr02,
attr03,
attr04,
attr05,
attr06,
attr07,
attr08,
attr09,
attr10
FROM mesh_250m
WHERE mesh_code = ?
LIMIT 1
"""
def open_read_only_connection() -> duckdb.DuckDBPyConnection:
"""読み取り専用DuckDB接続を作成する。
Returns
-------
duckdb.DuckDBPyConnection
Web APIから参照するための読み取り専用接続。
Notes
-----
ここでは分かりやすさを優先し、リクエストごとに接続を作る。
高負荷環境ではFastAPIのワーカー数、接続の持ち方、キャッシュを
実測しながら調整する。
"""
if not DATABASE_PATH.exists():
raise RuntimeError(
f"DuckDBファイルが見つかりません: {DATABASE_PATH}"
)
connection = duckdb.connect(
str(DATABASE_PATH),
read_only=True,
)
# 空間関数をこのAPIで使う場合に備えてSpatial拡張を読み込む。
# INSTALLはDockerイメージ作成時などに済ませておく。
connection.execute("LOAD spatial")
return connection
@app.get("/api/v1/meshes/{mesh_code}")
def get_mesh_detail(
mesh_code: Annotated[
str,
ApiPath(
min_length=10,
max_length=10,
pattern=r"^[0-9]{10}$",
description="250mメッシュコード",
),
],
) -> dict[str, str | None]:
"""指定した250mメッシュの10属性を返す。
Parameters
----------
mesh_code : str
10桁の250mメッシュコード。
Returns
-------
dict[str, str | None]
メッシュコードとattr01～attr10。
Raises
------
HTTPException
対象メッシュが存在しない場合は404を返す。
"""
# SQLへ文字列を直接連結しない。
# ?プレースホルダーへ値を渡すことで、検索条件を安全に扱う。
with closing(open_read_only_connection()) as connection:
row = connection.execute(
DETAIL_SQL,
[mesh_code],
).fetchone()
if row is None:
raise HTTPException(
status_code=404,
detail="mesh not found",
)
# SELECTの列順と合わせてJSONを組み立てる。
# 将来属性が増える場合はPydanticモデルへ切り出してもよい。
column_names = [
"mesh_code",
"attr01",
"attr02",
"attr03",
"attr04",
"attr05",
"attr06",
"attr07",
"attr08",
"attr09",
"attr10",
]
return dict(zip(column_names, row, strict=True))
このAPIは、地図表示用の大量ポリゴンを返すためではなく、クリックされた1メッシュの詳細情報を返すために使います。
OpenLayersから詳細属性を取得する
PMTiles上の250mメッシュをクリックしたら、タイルに入っているmesh_codeを使ってFastAPIへ問い合わせます。
JavaScriptコード全文を表示する
/**
* HTMLへ文字列を表示する前に、特殊文字をエスケープする。
*
* APIから取得した文字列をそのままinnerHTMLへ渡すと、
* データ中にHTMLが含まれていた場合に意図しない表示になる。
* そのため、詳細パネルへ入れる前に最低限のエスケープを行う。
*
* @param {unknown} value APIから取得した属性値
* @returns {string} HTML表示用にエスケープした文字列
*/
function escapeHtml(value) {
const text = String(value ?? "");
return text
.replaceAll("&", "&")
.replaceAll("<", "<")
.replaceAll(">", ">")
.replaceAll('"', """)
.replaceAll("'", "'");
}
/**
* mesh_codeを使って詳細属性APIを呼び出す。
*
* @param {string} meshCode 250mメッシュコード
* @returns {Promise<Object>} attr01～attr10を含むJSON
*/
async function fetchMeshDetail(meshCode) {
// パスへ埋め込む値なのでencodeURIComponentを通す。
const encodedMeshCode = encodeURIComponent(meshCode);
const response = await fetch(
`/api/v1/meshes/${encodedMeshCode}`,
{
headers: {
Accept: "application/json",
},
},
);
if (!response.ok) {
throw new Error(
`メッシュ詳細の取得に失敗しました: HTTP ${response.status}`,
);
}
return await response.json();
}
/**
* OpenLayersのクリックイベントから250mメッシュを取得し、
* 詳細属性を右側パネルへ表示する。
*/
map.on("singleclick", async (event) => {
// クリック位置にある最初のFeatureだけを使う。
// レイヤーが複数ある場合はlayerFilterを追加して対象を絞る。
const feature = map.forEachFeatureAtPixel(
event.pixel,
(candidate) => candidate,
);
if (!feature) {
return;
}
const meshCode = feature.get("mesh_code");
// 背景地図など別Featureをクリックした場合を考慮する。
if (!meshCode) {
return;
}
const panel = document.getElementById("mesh-detail");
panel.textContent = "読み込み中...";
try {
const detail = await fetchMeshDetail(meshCode);
// attr01～attr10は最大50文字だが、すべての属性を
// ベクタータイルへ入れず、このタイミングで初めて取得する。
panel.innerHTML = `
<h3>${escapeHtml(detail.mesh_code)}</h3>
<dl>
<dt>属性01</dt><dd>${escapeHtml(detail.attr01)}</dd>
<dt>属性02</dt><dd>${escapeHtml(detail.attr02)}</dd>
<dt>属性03</dt><dd>${escapeHtml(detail.attr03)}</dd>
<dt>属性04</dt><dd>${escapeHtml(detail.attr04)}</dd>
<dt>属性05</dt><dd>${escapeHtml(detail.attr05)}</dd>
<dt>属性06</dt><dd>${escapeHtml(detail.attr06)}</dd>
<dt>属性07</dt><dd>${escapeHtml(detail.attr07)}</dd>
<dt>属性08</dt><dd>${escapeHtml(detail.attr08)}</dd>
<dt>属性09</dt><dd>${escapeHtml(detail.attr09)}</dd>
<dt>属性10</dt><dd>${escapeHtml(detail.attr10)}</dd>
</dl>
`;
} catch (error) {
console.error(error);
panel.textContent = "詳細情報を取得できませんでした。";
}
});
この構成なら、全国を表示している間に長い属性を何度も通信せずに済みます。
DuckDBを使う場合の容量はどう考えるか
DuckDBも単一ファイルなので、GeoPackageと同じように「何GBになるか」が気になります。
ただし、DuckDBはGIS交換フォーマットではなくデータベースです。内部の列圧縮、データ型、文字列の重複、インデックス、テーブル構成によってファイルサイズが変わるため、ShapefileやGeoJSONと同じ条件で単純比較するのは少し違います。
今回のようなデータなら、容量計画では次の順で考えます。
- GeoParquetの実測値を原本容量の基準にする
- 同じデータをDuckDBへロードして
.duckdbファイルを実測する - R-treeや通常インデックスを作成した後でもう一度測る
- 更新用の一時領域とバックアップ分を加える
つまり、DuckDBについては「全国で○GB」と固定値を書くより、実際の属性分布を持つ1県分または数十万件で試し、その倍率を見る方が確実です。
今回の記事のベンチマークスクリプトへDuckDB測定を追加するなら、次のようにします。
Pythonコード全文を表示する
"""GeoParquetをDuckDBへ取り込み、DBファイル容量を確認する例。"""
from pathlib import Path
import duckdb
DATABASE_PATH = Path("data/benchmark/mesh_250m.duckdb")
GEOPARQUET_PATH = Path("data/benchmark/mesh_250m.parquet")
def build_duckdb() -> int:
"""GeoParquetをDuckDBへ取り込み、ファイルサイズを返す。
Returns
-------
int
DuckDBファイルのバイト数。
"""
# 古い測定ファイルが残っていると正しい比較にならないため削除する。
if DATABASE_PATH.exists():
DATABASE_PATH.unlink()
connection = duckdb.connect(str(DATABASE_PATH))
try:
# Spatial拡張を読み込む。
# INSTALLは初回だけ必要なので、測定環境では事前実行でもよい。
connection.execute("INSTALL spatial")
connection.execute("LOAD spatial")
# GeoParquetをDuckDBの実テーブルへ取り込む。
# read_parquetを使うことで、Parquetの列をそのまま読み込める。
connection.execute(
"""
CREATE TABLE mesh_250m AS
SELECT *
FROM read_parquet(?)
""",
[str(GEOPARQUET_PATH)],
)
# WebのBBOX検索を想定し、geometry列へR-treeを作成する。
# インデックス作成前後の容量を別々に測ると、
# R-treeがどの程度増分になるか確認できる。
connection.execute(
"""
CREATE INDEX mesh_250m_rtree
ON mesh_250m USING RTREE (geometry)
"""
)
# CHECKPOINTで変更内容をデータベースファイルへ反映してから測る。
connection.execute("CHECKPOINT")
finally:
connection.close()
return DATABASE_PATH.stat().st_size
全国規模へ進む前に、この実測を入れておくとディスク計画がかなり現実的になります。
PostGISは最初から入れない
PostGISを使わない、という意味ではありません。
今回の全国250mメッシュは、基本的に次の使い方を想定しています。
原本を作る
↓
定期的に更新する
↓
DuckDBで集計・検索する
↓
PMTilesを作る
↓
Webで参照する
この段階では、PostgreSQL＋PostGISまで入れると、構成要素と運用項目が増えます。
そのため、初期構成はDuckDB Spatialまでにします。
一方、次のような要件が出てきたらPostGISを検討します。
- 複数利用者が同じGISデータを編集する
- Web画面から地物や属性を頻繁に更新する
- 複数のAPIやバッチ処理が同時に書き込む
- トランザクションを使った更新管理が必要になる
- ユーザー権限や業務データとGISデータを一体で管理する
その段階では、例えば次の構成へ広げます。
PostGISを追加してもDuckDBを捨てる必要はありません。
DuckDB Spatial
└─ 大量データ処理・ETL・分析・タイル生成前処理
PostgreSQL + PostGIS
└─ 多人数編集・同時更新・業務運用
という分担にできます。
フォルダー構成
今回の記事の構成をDocker Composeへ載せるなら、例えば次のように分けます。
japan_250m_mesh_web/
├── compose.yaml
├── .env
├── README.md
│
├── data/
│ ├── raw/
│ │ └── README.md
│ │
│ ├── master/
│ │ └── mesh_250m/
│ │ └── *.parquet
│ │
│ ├── runtime/
│ │ └── japan_250m_mesh.duckdb
│ │
│ └── published/
│ └── japan_250m.pmtiles
│
├── scripts/
│ ├── build_geoparquet.py
│ ├── build_duckdb.py
│ ├── build_lod.py
│ ├── build_pmtiles.sh
│ └── benchmark_formats.py
│
├── api/
│ ├── Dockerfile
│ ├── requirements.txt
│ └── app/
│ ├── main.py
│ ├── database.py
│ ├── schemas.py
│ └── routers/
│ └── meshes.py
│
└── web/
├── Dockerfile
├── nginx.conf
└── src/
├── main.js
└── style.css
役割が分かれているので、後からPostGISを追加する場合も、db/サービスと同期処理を追加すれば済みます。
同じメッシュ形状を何度も保存しない
全国250mメッシュを今後いろいろな用途へ使うなら、データセットごとに同じ600万ポリゴンを複製しない方がよいです。
例えば人口、防災、降雨、地形、土地利用を追加するときは、mesh_codeを共通キーにします。
mesh_250m
├─ mesh_code
└─ geometry
population
├─ mesh_code
└─ population
hazard
├─ mesh_code
└─ hazard_class
rainfall
├─ mesh_code
└─ rain_mm
DuckDBなら、このような表をSQLで結合できます。
-- geometryを持つ共通メッシュと人口属性をmesh_codeで結合する。
-- 各業務テーブルへ同じポリゴンを複製しないため、
-- ストレージと更新処理を減らせる。
SELECT
m.mesh_code,
m.geometry,
p.population
FROM mesh_250m AS m
LEFT JOIN population AS p
USING (mesh_code);
メッシュコードから矩形geometryを再構築できる設計なら、さらに形状の重複を減らす方法もあります。
今回の結論
全国陸域の250mメッシュを約600万件、ポリゴン、文字列属性10項目、各項目最大50文字という条件で考えると、非圧縮系のGIS形式では数GBから10GBを超える規模になります。
今回の試算では、10項目をすべて50文字まで埋めた上限側の条件で、次の程度になりました。
| フォーマット | 半角英数字中心 | 日本語中心 |
|---|---|---|
| Shapefile一式 | 約4.1GB | 約10.1GB |
| FlatGeobuf | 約4.4GB | 約10.4GB |
| GeoPackage | 約4.6GB | 約12.8GB |
| GeoJSON | 約5.4GB | 約11.4GB |
| GeoParquet＋ZSTD | 約2～4GB | 約5～9GB |
そして、実際にWebで使う構成は次の形にします。
PostGISは最初から入れません。
まずはDuckDB Spatialまでで運用し、将来、多人数による編集・同時更新・トランザクション管理が必要になった段階でPostgreSQL＋PostGISを追加します。
今回のような全国規模の固定メッシュでは、すべてを1つの仕組みへ押し込むより、
- GeoParquet：保管
- DuckDB Spatial：処理
- PMTiles：配信
- FastAPI：詳細検索
- OpenLayers：表示
- PostGIS：必要になったら運用編集DBとして追加
という役割分担にした方が、後から構成を変更しやすそうです。
次は、この構成をDocker Compose上に実装して、実際の250mメッシュをGeoParquetへ変換し、DuckDBへ取り込み、OpenLayersで表示するところまで試してみたいと思います。
参考資料
- 総務省統計局「地域メッシュ統計」
https://www.stat.go.jp/data/mesh/ - 国土地理院「全国都道府県市区町村別面積調」
https://www.gsi.go.jp/KOKUJYOHO/MENCHO-title.htm - GDAL GeoPackage driver
https://gdal.org/en/stable/drivers/vector/gpkg.html - GDAL ESRI Shapefile driver
https://gdal.org/en/stable/drivers/vector/shapefile.html - GeoParquet specification
https://geoparquet.org/ - DuckDB Spatial Extension
https://duckdb.org/docs/stable/core_extensions/spatial/overview.html - DuckDB Spatial Functions
https://duckdb.org/docs/stable/core_extensions/spatial/functions.html - DuckDB R-Tree Indexes
https://duckdb.org/docs/stable/core_extensions/spatial/r-tree_indexes.html - DuckDB Parquet
https://duckdb.org/docs/stable/data/parquet/overview.html - PMTiles
https://docs.protomaps.com/pmtiles/ - Tippecanoe
https://github.com/felt/tippecanoe - OpenLayers
https://openlayers.org/ - FastAPI
https://fastapi.tiangolo.com/ - PostGIS
https://postgis.net/docs/