e-Statの全国1kmメッシュをGeoParquetにまとめる ― 防災Web地図 Step 1
前回の実装フローに従って、今回はStep 1を実装していきます。
まずは、全国の防災データを1kmメッシュで整理するため、最初に共通のメッシュマスターを作ります。
メッシュを作るだけなら、緯度経度の範囲から格子を発生させる方法もあります。ただ、後で国勢調査や地域メッシュ統計と結び付けるなら、独自に作った1,000m角の格子ではなく、e-Statが公開している3次メッシュ境界を使う方が扱いやすいです。
今回のStep 1では、e-Statから取得した第1次メッシュ単位のShapefileを全国分まとめ、検査済みのGeoParquetへ変換します。
Step 1のフロー
e-Stat 3次メッシュ境界 Shapefile
↓
入力ファイル、属性、CRS、geometryを検査
↓
メッシュコードを8桁の文字列へ統一
↓
全国分を結合、重複検査、並べ替え
↓
GeoParquet 1.1へ保存
↓
読み戻し検査
↓
manifest.jsonとvalidation_report.jsonを作成
この記事は、次の5段階で防災Web地図を作るうちのStep 1です。
防災Web地図を作る
Step 1 全国1kmメッシュをGeoParquetへ変換
Step 2 洪水・土砂・津波などを正規化し、1kmメッシュへ集計
Step 3 TippecanoeでPMTilesを作成
Step 4 NginxとOpenLayersで静的表示
Step 5 DuckDBとFastAPIで詳細・統計APIを追加
今回作る成果物
Step 1の出力は3ファイルです。
workspace/artifacts/step1/
├─ mesh_master_1km.parquet
├─ manifest.json
└─ validation_report.json
mesh_master_1km.parquetが、Step 2以降で使う全国共通のメッシュマスターです。
manifest.jsonには、入力したShapefile、SHA-256、実行環境、成果物版を記録します。validation_report.jsonには、GeoParquetを書き出した後に読み戻して確認した結果を残します。
全国分を一度作った後は、GeoParquetだけ見ても、どのZIPを使ったのか分かりません。元データを更新したときに差を追えるよう、来歴と検査結果も成果物として扱います。
元データ
使用するのは、e-Statの統計地理情報システムで公開されている次のデータです。
境界データ
3次メッシュ（1kmメッシュ）
世界測地系緯度経度
Shapefile
ダウンロード画面は次です。
e-Statの説明では、3次メッシュ境界は第1次地域区画単位で提供されています。作成範囲は、国土にかかる第2次地域区画で、第2次地域区画の全域が水面の場合は除かれます。
したがって、ここでいう「全国1kmメッシュ」は、日本周辺の海域を無制限に埋めた格子ではありません。e-Statが提供する国土側のメッシュ境界を全国分まとめたものです。
1kmメッシュの定義書では、Shapefile名と属性が次のように示されています。
| 列名 | 内容 |
|---|---|
KEY_CODE |
8桁のメッシュコード |
MESH1_ID |
4桁の1次メッシュコード |
MESH2_ID |
2桁の2次メッシュ部分コード |
MESH3_ID |
2桁の3次メッシュ部分コード |
OBJ_ID |
9桁の通し番号 |
ファイル名は MESH0xxxx.shp です。xxxxには1次メッシュコードが入ります。
1,000m角の格子を自作しない理由
3次メッシュは、地域メッシュコードで定められた緯度経度方向の区画です。Web地図で使うEPSG:3857上に1,000m間隔の正方形を並べても、e-Statの3次メッシュとは一致しません。
独自格子でも表示はできますが、後から人口統計や事業所統計を結合するときに、公式のメッシュコードをそのまま使えなくなります。今回は防災データだけでなく、人口、世帯、建物なども重ねる予定なので、最初から公式境界を基準にします。
Step 1の処理範囲
ここでは防災区域との空間結合は行いません。洪水や土砂災害のポリゴンを重ねる処理はStep 2へ分けます。
Step 1の役割は、今後の処理で何度も使うメッシュ境界を一度きちんと整えることです。
実行環境
Docker Composeで実行します。
参考のためにZIPファイルを置いときますので、ダウンロードして使ってください。
Python 3.13
GeoPandas 1.1.4
Pyogrio 0.13.0
PyArrow 25.0.0
Shapely 2.1系
PyProj 3.7系
主要なI/Oライブラリは版を固定しています。GeoParquetはファイル形式の再現性が重要なので、後でライブラリが更新されても、同じStep 1イメージを使えば同じ条件で作り直せるようにしています。
GeoPandasのto_parquet()では、WKB、ZSTD、GeoParquet 1.1、bbox coveringを指定します。
frame.to_parquet(
temporary_path,
index=False,
compression="zstd",
compression_level=6,
geometry_encoding="WKB",
write_covering_bbox=True,
schema_version="1.1.0",
row_group_size=row_group_size,
)
WKBを選んだのは、GeoArrowよりも読み込み側の対応範囲を優先したためです。bbox coveringは、地物ごとのxmin、ymin、xmax、ymaxをParquet列として持たせます。Step 2で範囲を絞って読む場合や、GeoPandasのBBOX読み込みを試す場合に使えます。
フォルダー構成
japan_disaster_pmtiles_step1_r001/
├─ README.md
├─ VALIDATION.md
├─ docker-compose.yml
├─ Makefile
├─ .env.example
├─ pyproject.toml
│
├─ builder/
│ ├─ Dockerfile
│ └─ requirements.txt
│
├─ src/mesh_step1/
│ ├─ archives.py
│ ├─ cli.py
│ ├─ config.py
│ ├─ demo.py
│ ├─ errors.py
│ ├─ hashing.py
│ ├─ inspect_output.py
│ ├─ json_io.py
│ ├─ logging_utils.py
│ ├─ mesh_code.py
│ ├─ metadata.py
│ ├─ pipeline.py
│ ├─ sources.py
│ ├─ transform.py
│ ├─ validation.py
│ └─ writer.py
│
├─ tests/
├─ contracts/
├─ docs/
├─ scripts/
└─ workspace/
├─ downloads/step1/
├─ raw/mesh/
├─ raw/mesh_demo/
├─ artifacts/step1/
├─ artifacts/step1_demo/
└─ logs/
原本、展開後データ、成果物を分けています。
workspace/downloads/step1
e-Statから取得したZIP
workspace/raw/mesh
ZIPを展開したShapefile一式
workspace/artifacts/step1
Step 2へ渡す検査済み成果物
原本ZIPを加工済みデータと混ぜないため、処理をやり直すときにも判断しやすくなります。
まずデモデータで確認する
全国分のZIPをそろえる前に、小さなデモで処理を通します。
cp .env.example .env
make demo
make demoでは、次を順番に実行します。
Dockerイメージ作成
↓
架空の1kmメッシュShapefileを生成
↓
2つの1次メッシュ単位ファイルを結合
↓
GeoParquetを作成
↓
読み戻し検査
↓
manifestの要点を表示
デモは32セルです。経度139度付近の1次メッシュ境界をまたぐように配置し、複数ファイルの結合を確認できるようにしています。
workspace/raw/mesh_demo/
├─ MESH05238/
│ ├─ MESH05238.shp
│ ├─ MESH05238.shx
│ ├─ MESH05238.dbf
│ └─ MESH05238.prj
└─ MESH05239/
├─ MESH05239.shp
├─ MESH05239.shx
├─ MESH05239.dbf
└─ MESH05239.prj
成果物は実データと分けて保存します。
workspace/artifacts/step1_demo/
├─ mesh_master_1km.parquet
├─ manifest.json
└─ validation_report.json
デモデータは処理確認用の架空データです。地域判定や防災判断には使いません。
e-StatのZIPを配置する
ダウンロードしたZIPは、展開せずに次へ置きます。
workspace/downloads/step1/
├─ MESH0xxxx.zip
├─ MESH0yyyy.zip
└─ ...
実際のZIP名は配布時の名称を保って構いません。展開後にMESH0xxxx.shpを再帰検索します。
全国分はファイル数が多いため、一度に手作業で展開すると、同じZIPを二重に展開したり、一部だけ別フォルダーへ置いたりしやすくなります。そこで、ZIPの展開もStep 1のコマンドに含めました。
make build
make prepare
既存の展開先を削除して作り直す場合は、次を実行します。
./scripts/prepare_archives.sh --overwrite
ZIPをそのまま展開しない
archives.pyでは、ZIP内のパスを確認してから展開します。
def _safe_member_path(destination: Path, member_name: str) -> Path:
"""ZIP内パスを検査し、展開先の絶対パスを返す。"""
member = PurePosixPath(member_name)
if member.is_absolute() or ".." in member.parts:
raise SourceDataError(
f"ZIP内に安全でないパスがあります: {member_name}"
)
target = (destination / Path(*member.parts)).resolve()
target.relative_to(destination.resolve())
return target
公的機関から取得したZIPだけを使う予定でも、展開処理は今後ほかのデータにも流用します。../を含むパスやシンボリックリンクを許可しない形にしてあります。
Shapefileの構成を先に検査する
Shapefileは.shpだけでは完結しません。
Step 1では、通常次の4ファイルを必須としています。これは気をつけてくださいね。特にprjファイルがない場合がありますので。必ず4ファイルが必須です。
.shp geometry
.shx geometryの索引
.dbf 属性
.prj 座標参照系
.cpgがある場合は、ハッシュ対象へ含めます。.cpgは文字コード指定ファイルです。
REQUIRED_SIDECAR_EXTENSIONS = (
".shp",
".shx",
".dbf",
".prj",
)
OPTIONAL_SIDECAR_EXTENSIONS = (".cpg",)
Windowsで展開すると、拡張子の大文字・小文字が混在することがあります。Path.with_suffix()で決め打ちせず、同じstemのファイルを走査して拡張子を小文字へそろえて判定します。
読み込み前にはpyogrio.read_info()で次を確認します。
CRS
geometry型
属性列名
地物数
info = pyogrio.read_info(shapefile)
crs = info.get("crs")
geometry_type = str(info.get("geometry_type") or "")
feature_count = int(info.get("features", -1))
fields = [str(value) for value in info.get("fields", [])]
この段階で.prjがない、KEY_CODEがない、Polygonではない、といった入力は止めます。
必要な属性だけを読む
全国分のメッシュ境界は、ファイルごとに順番に読みます。
frame = pyogrio.read_dataframe(
source.shapefile,
columns=actual_columns,
use_arrow=True,
force_2d=True,
)
必要なのは、公式定義書の5列とgeometryです。使わない列を読み込まないようにしています。
DockerイメージにはPyArrowを入れているため、通常はArrow経由で読みます。一方、メッシュコードの単体テストだけをローカルで行う場合にPyArrowがないこともあるため、その場合は通常のGDAL読み込みへ切り替えます。
use_arrow = importlib.util.find_spec("pyarrow") is not None
本番処理の仕様を緩めるためではなく、GeoParquetを書かないテストまでPyArrow必須にしないための切り分けです。
メッシュコードは文字列で扱う
KEY_CODEは8桁の数字ですが、計算値ではなく識別コードです。
DBFの読み込み方によって、文字列、整数、浮動小数点として見える可能性があるため、Step 1では固定桁の文字列へ統一します。
def normalize_fixed_width_identifier(
values: pd.Series,
*,
width: int,
field_name: str,
) -> pd.Series:
"""識別コードを数字だけの固定桁文字列へ正規化する。"""
normalized = (
values.astype("string")
.str.strip()
.str.replace(r"\.0$", "", regex=True)
.str.zfill(width)
)
invalid_mask = ~normalized.str.fullmatch(
rf"\d{{{width}}}",
na=False,
)
if invalid_mask.any():
raise SourceDataError(
f"{field_name}を{width}桁へ正規化できません。"
)
return normalized
8桁へそろえるだけではなく、部分コードとの対応も確認します。
MESH1_ID == KEY_CODEの先頭4桁
MESH2_ID == KEY_CODEの5～6桁
MESH3_ID == KEY_CODEの7～8桁
ファイル名のxxxx == MESH1_ID
ZIPの取り違えや属性の不整合を、全国結合後まで持ち越さないためです。
CRSを決め打ちしない
入力ShapefileのCRSは、.prjをGDALで読み取ります。
if frame.crs is None:
raise SourceDataError(
f"読み込み後のCRSが未設定です: {source.shapefile}"
)
出力はEPSG:4326へ統一します。
if not frame.crs.equals(target_crs):
frame = frame.to_crs(target_crs)
ここでset_crs()は使いません。set_crs()は座標値を変えずにCRS情報だけを付ける処理です。入力と出力のCRSが違う場合は、座標値を変換するto_crs()を使います。
CRSの文字列は、同じ座標参照系でもWKTとEPSG表記で異なることがあります。文字列を直接比較せず、PyProjが解釈したCRS同士で比較しています。
geometryは自動修復しない
メッシュ境界は、Step 2の空間集計で基準になります。
不正geometryを見つけたときにmake_valid()で自動修復すると、公式境界をどのように変えたのか分かりにくくなります。Step 1では次を確認し、異常があれば入力ファイルを特定して停止します。
geometry欠損
空geometry
不正geometry
Polygon / MultiPolygon以外
invalid_count = int((~frame.geometry.is_valid).sum())
if invalid_count:
raise SourceDataError(
f"不正geometryがあります: {source_name} / {invalid_count:,}件"
)
防災区域側では、作成機関や年度によってgeometryの状態が異なるため、Step 2で修復方針を別に決めます。基盤メッシュと主題データで同じ扱いにしないようにしています。
全国分を結合する
検査済みのShapefileをGeoDataFrameへ読み、最後に結合します。
merged = gpd.GeoDataFrame(
pd.concat(frames, ignore_index=True, copy=False),
geometry="geometry",
crs=target_crs,
)
結合後にmesh_codeの重複を確認します。
duplicate_mask = merged["mesh_code"].duplicated(keep=False)
if duplicate_mask.any():
raise SourceDataError(
"全国結合後に重複するメッシュコードがあります。"
"入力ZIPの二重展開も確認してください。"
)
同じZIPを別フォルダーへ二重展開していると、ファイル単位の検査は通っても、全国結合後に同じKEY_CODEが現れます。このため、重複検査は結合後にも行います。
メッシュコード順に並べる
結合後はmesh_code順に並べます。
merged = (
merged.sort_values("mesh_code", kind="stable")
.reset_index(drop=True)
)
メッシュコードは、1次、2次、3次メッシュの階層を含みます。同じ地域のセルが近くに並ぶため、完全にランダムな順序より、ParquetのRow Group単位で範囲を絞るときに扱いやすくなります。
GeoHashやHilbert曲線による並べ替えも考えられますが、Step 1では公式コードを見れば順序が分かる単純な構成にしています。PMTiles用の表示順序や低ズーム集計は、後のStepで別に調整します。
出力列
GeoParquetには次の列を保存します。
| 列名 | 内容 |
|---|---|
feature_id |
MVTやAPIでも使える安定した数値ID |
mesh_code |
8桁の3次メッシュコード |
mesh1_code |
4桁の1次メッシュコード |
mesh2_code |
6桁の2次メッシュコード |
mesh3_code |
8桁の3次メッシュコード |
center_lon |
メッシュ境界の東西中点 |
center_lat |
メッシュ境界の南北中点 |
source_obj_id |
e-Stat原本のOBJ_ID
|
source_mesh_file |
原本Shapefileの相対パス |
geometry |
EPSG:4326のメッシュポリゴン |
bbox |
GeoParquet 1.1のcovering bbox |
feature_idには8桁メッシュコードの数値を使います。
merged.insert(
0,
"feature_id",
merged["mesh_code"].astype("int64"),
)
処理のたびに1から連番を振ると、入力範囲を追加しただけで既存地物のIDが変わることがあります。メッシュコードから作れば、同じセルは同じIDになります。
mesh2_codeは、原本の2桁部分コードではなく、1次メッシュを含む6桁の完全コードとして保存します。Step 2以降でprefix集計をしやすくするためです。
一時ファイルを経由して保存する
GeoParquetの保存途中で処理が止まった場合、出力先に不完全なファイルを残したくありません。
そこで、一時ファイルへ書き出した後、完了したものだけをos.replace()で正式名へ置き換えます。
temporary_path = output_path.with_name(
f".{output_path.stem}.{uuid.uuid4().hex}.tmp.parquet"
)
try:
frame.to_parquet(
temporary_path,
index=False,
compression="zstd",
compression_level=6,
geometry_encoding="WKB",
write_covering_bbox=True,
schema_version="1.1.0",
row_group_size=row_group_size,
)
os.replace(temporary_path, output_path)
finally:
temporary_path.unlink(missing_ok=True)
既存成果物は、初期状態では上書きしません。意図して作り直すときだけ--overwriteを付けます。
GeoParquetを書いた後に読み戻す
ファイルが存在し、容量が0より大きいだけでは、Step 2へ渡しません。
PyArrowでParquetとGeoParquetメタデータを確認し、GeoPandasでgeometryを読み戻します。
主な検査項目は次です。
ファイルが存在する
行数が1件以上ある
書き込み前後の行数が一致する
必須列がある
GeoParquet versionが1.1.0
primary geometryがgeometry列
geometry encodingがWKB
bbox coveringメタデータがある
mesh_codeが8桁
mesh_codeとfeature_idが一意
階層コードが一致する
CRSがEPSG:4326
geometryに欠損、空、不正形状がない
検査結果はJSONへ保存します。
{
"name": "mesh_code_unique",
"status": "passed",
"message": "mesh_codeに重複はありません。",
"value": 0
}
1項目でもfailedになった場合、処理は異常終了します。失敗した検査報告は残るため、原因を確認できます。
manifestに記録する内容
manifest.jsonには次を残します。
成果物版
開始・終了時刻
処理時間
出力CRS
Row Groupサイズ
Python・GDAL・パッケージ版
入力Shapefileごとの相対パス
入力地物数
Shapefile一式の容量とSHA-256
出力GeoParquetの容量とSHA-256
検査結果の要約
Shapefileは複数ファイルで1データセットなので、.shpだけではなく、.shx、.dbf、.prj、存在する場合は.cpgもまとめてハッシュします。
実データで実行する
.env.exampleをコピーします。
cp .env.example .env
成果物版を設定します。
STEP1_IMAGE_TAG=r001
STEP1_ARTIFACT_VERSION=20260804_r001
STEP1_ROW_GROUP_SIZE=100000
Dockerイメージを作ります。
make build
ZIPを展開します。
make prepare
GeoParquetを作ります。
make step1
作成済み成果物を別コマンドで再検査します。
make validate
要約を表示します。
make inspect
テストを実行します。
make test
再実行するとき
既存成果物がある場合、通常のmake step1は停止します。
内容を確認して置き換える場合だけ、次を使います。
./scripts/run_step1.sh --overwrite
ZIPの展開先も作り直す場合は、次です。
./scripts/prepare_archives.sh --overwrite
実データとデモデータは別フォルダーなので、デモだけ削除できます。
make clean-demo
よくあるエラー
MESH0*.shp が見つかりません
ZIPの配置場所と展開処理を確認します。
find workspace/raw/mesh -name 'MESH0*.shp' | head
見つからなければ、次を再実行します。
make prepare
.prjがありません
.shpだけを別フォルダーへ移していないか確認します。Shapefile一式を同じフォルダーへ置きます。
MESH0xxxx.shp
MESH0xxxx.shx
MESH0xxxx.dbf
MESH0xxxx.prj
全国結合後に重複するメッシュコードがあります
同じZIPを複数の場所へ展開している可能性があります。
find workspace/raw/mesh -name 'MESH0xxxx.shp'
入力を整理した後、展開先を作り直します。
./scripts/prepare_archives.sh --overwrite
メモリ不足になる
全国分は一度GeoDataFrameへ結合するため、Dockerへ割り当てたメモリが少ないと停止することがあります。
まずDocker DesktopやWSL2のメモリ上限を確認します。それでも厳しい場合は、Step 1の次版で1次メッシュ単位のGeoParquet Datasetへ分割し、最後にDuckDBで参照する構成へ変更できます。
今回はStep間の受け渡しを単純にするため、1つのGeoParquetへまとめています。
Step 2へ渡すもの
Step 2が読むのは次です。
workspace/artifacts/step1/mesh_master_1km.parquet
workspace/artifacts/step1/manifest.json
Step 1のPythonソースや展開途中のShapefileを、Step 2から直接参照しません。
Step 1の入力・処理方法が変わる
↓
同じデータ契約でGeoParquetを出す
↓
Step 2は変更しない
この分け方にしておくと、全国メッシュの作り直しと、防災区域の空間集計を別々に検証できます。
まとめ
Step 1では、e-Statの3次メッシュ境界を全国共通のGeoParquetへまとめました。
ポイントは、Shapefileを読んで保存する部分より、前後の検査です。
公式属性の桁数と対応を確認する
CRSを決め打ちしない
不正geometryを黙って修復しない
全国結合後の重複を確認する
保存途中のファイルを正式成果物にしない
書き出したGeoParquetを読み戻す
入力と出力のハッシュを残す
次のStep 2では、このmesh_master_1km.parquetへ洪水、土砂災害、津波などの公式GISデータを重ねます。元ポリゴンを残しながら、1kmメッシュごとの最大クラス、影響面積、面積率を作る予定です。
参考リンク
- e-Stat 統計地理情報システム データダウンロード
- e-Stat ダウンロードデータについて
- e-Stat メッシュ境界（1km）定義書
- GeoPandas
GeoDataFrame.to_parquet - Pyogrio Introduction
- GeoParquet 1.1.0 specification