全国版：国土基本図1/2500図郭をPythonで作成する

JGD2011対応：国土基本図1/2500図郭をPythonで作成する

平面直角座標系ごとに正式な図郭コード（HD***形式）付きの 1/2500 図郭データを自動生成する手順メモ。

1. 参考にした資料

本記事の図郭コード計算ロジックは、公共測量標準図式および下記記事の内容を参考にしています。

@shiba54 (Y S) さん： 「国土基本図の図郭」を作成する

実装は上記アイデアをもとに、”全国対応版”JGD2011の各平面直角座標系に対応するよう再構成したものです。

2. 前提知識のざっくり整理

2.1 平面直角座標系（JGD2011）

日本の平面直角座標系（Gauss–Krüger 投影）は I〜XIX の 19 系。 各系には中央経線と EPSG コードが割り当てられています。

• VII系：EPSG 6675（石川・富山・岐阜・愛知 ほか）
• X系：EPSG 6678（青森・秋田・岩手・山形・宮城 ほか）
• …といった形で全国をカバー

2.2 国土基本図の図郭体系（レベル2500）

公共測量標準図式に基づき、図郭は階層的に構成されています。

• 1/50,000：40km × 30km（2文字コード：AA〜TT）
• 1/5,000：4km × 3km（上記 50,000 を 10×10 分割 → 2桁追加）
• 1/2,500：2km × 1.5km（5,000 を 2×2 分割 → 1桁（1〜4）追加）

例：HD433 のようなコードが 1/2500 の図郭記号となります。

3. 必要な環境

Python 3.x と、次のライブラリを使用します。

pip install geopandas shapely pyproj

QGIS を使っている環境であれば、ほぼそのまま導入できるはずです。

4. 実装のポイント

4.1 系一覧テーブル

スクリプト内で、JGD2011 の各系ごとに EPSG と中央経線、そして「おおよその都道府県」を定義しておきます。 これをそのまま「プルダウン風」の選択肢として利用します。

4.2 図郭コードの算出

図郭コードは、@shiba54 さんの記事と同様に、以下のステップで計算します。

1. 1/50,000 図郭：40km × 30km 単位で行・列を求め、A〜T の 2文字に変換
2. 1/5,000 図郭：4km × 3km のブロック番号（0〜9, 0〜9）を追加
3. 1/2,500 図郭：2km × 1.5km で 2×2 分割し、1〜4 を付与

これを関数 kokudo2500_code(x_left, y_top) として実装し、 各セルの西端X座標と北端Y座標から HD*** 形式のコードを生成します。

5. コード全文

以下のスクリプトを 1 ファイルとして保存して実行します。

import math
import geopandas as gpd
from shapely.geometry import box
from pyproj import Transformer

# --- JGD2011 各系の定義（抜粋・備考付き） ---
SYSTEMS = [
    {"roman": "I",  "num": 1,  "epsg": 6669, "lon0": 129.5,      "pref": "長崎, 鹿児島西部"},
    {"roman": "II", "num": 2,  "epsg": 6670, "lon0": 131.0,      "pref": "福岡, 佐賀, 熊本, 大分, 宮崎, 鹿児島東部"},
    {"roman": "III","num": 3,  "epsg": 6671, "lon0": 132.1666667,"pref": "山口, 島根, 広島"},
    {"roman": "IV", "num": 4,  "epsg": 6672, "lon0": 133.5,      "pref": "香川, 愛媛, 徳島, 高知"},
    {"roman": "V",  "num": 5,  "epsg": 6673, "lon0": 134.3333333,"pref": "兵庫, 鳥取, 岡山"},
    {"roman": "VI", "num": 6,  "epsg": 6674, "lon0": 136.0,      "pref": "京都, 大阪, 奈良, 和歌山, 三重西部"},
    {"roman": "VII","num": 7,  "epsg": 6675, "lon0": 137.1666667,"pref": "石川, 富山, 岐阜, 愛知"},
    {"roman": "VIII","num":8,  "epsg": 6676, "lon0": 138.5,      "pref": "新潟, 長野, 山梨, 静岡"},
    {"roman": "IX", "num": 9,  "epsg": 6677, "lon0": 139.8333333,"pref": "東京, 神奈川, 埼玉, 茨城西部"},
    {"roman": "X",  "num":10,  "epsg": 6678, "lon0": 140.8333333,"pref": "青森, 秋田, 岩手, 山形, 宮城"},
    {"roman": "XI", "num":11,  "epsg": 6679, "lon0": 140.25,     "pref": "北海道西部"},
    {"roman": "XII","num":12,  "epsg": 6680, "lon0": 142.25,     "pref": "北海道中部〜東部"},
    {"roman": "XIII","num":13, "epsg": 6681, "lon0": 144.25,     "pref": "北海道東端"},
    {"roman": "XIV","num":14,  "epsg": 6682, "lon0": 142.0,      "pref": "伊豆諸島・小笠原付近"},
    {"roman": "XV", "num":15,  "epsg": 6683, "lon0": 127.5,      "pref": "沖縄本島周辺"},
    {"roman": "XVI","num":16, "epsg": 6684, "lon0": 124.0,      "pref": "宮古・八重山"},
    {"roman": "XVII","num":17, "epsg": 6685, "lon0": 131.0,      "pref": "大東諸島"},
    {"roman": "XVIII","num":18,"epsg": 6686, "lon0": 136.0,      "pref": "小笠原諸島"},
    {"roman": "XIX","num":19,  "epsg": 6687, "lon0": 154.0,      "pref": "南鳥島"},
]

CRS_GEOG = "EPSG:6668"  # JGD2011 緯度経度
CELL_W = 2000           # 2500図郭 東西 2km
CELL_H = 1500           # 2500図郭 南北 1.5km

def kokudo2500_code(x_left: int, y_top: int) -> str:
    # 1/50,000
    col50000 = int((x_left + 160000) // 40000)
    row50000 = int((300000 - y_top) // 30000)
    c1 = chr(ord("A") + row50000)
    c2 = chr(ord("A") + col50000)
    # 1/5,000
    n5000 = int(((300000 - y_top) // 3000) % 10)
    e5000 = int(((x_left + 160000) // 4000) % 10)
    # 1/2,500
    rn = int(((300000 - y_top) // 1500) % 2)
    cn = int(((x_left + 160000) // 2000) % 2)
    quad = rn * 2 + cn + 1
    return f"{c1}{c2}{n5000}{e5000}{quad}"

# --- 系の選択（疑似プルダウン） ---
print("==== JGD2011 平面直角座標系を選択 ====")
for s in SYSTEMS:
    print(f"{s['num']:>2}: 系{s['roman']:<3} EPSG:{s['epsg']}  備考:{s['pref']}")

choice = input("\\n系の番号 or EPSGコード を入力（例: 7 or 6675）: ").strip()

selected = None
for s in SYSTEMS:
    if choice == str(s["num"]) or choice == str(s["epsg"]):
        selected = s
        break

if not selected:
    raise SystemExit("不正な入力です。1〜19の番号か EPSGコードを指定してください。")

epsg = selected["epsg"]
lon0 = selected["lon0"]
print(f"--> {selected['roman']}系 (EPSG:{epsg}) / 想定エリア: {selected['pref']}")

# --- 対象経緯度範囲（中央経線 ±4度, 緯度 24〜46.5度 をデフォルト） ---
span = 4.0
min_lon = lon0 - span
max_lon = lon0 + span
min_lat = 24.0
max_lat = 46.5

CRS_PLANE = f"EPSG:{epsg}"
transformer = Transformer.from_crs(CRS_GEOG, CRS_PLANE, always_xy=True)

x1, y1 = transformer.transform(min_lon, min_lat)
x2, y2 = transformer.transform(max_lon, max_lat)

min_x = math.floor(min(x1, x2) / CELL_W) * CELL_W
max_x = math.ceil (max(x1, x2) / CELL_W) * CELL_W
min_y = math.floor(min(y1, y2) / CELL_H) * CELL_H
max_y = math.ceil (max(y1, y2) / CELL_H) * CELL_H

geoms = []
codes = []

y = min_y
while y < max_y:
    x = min_x
    while x < max_x:
        left   = int(x)
        right  = int(x + CELL_W)
        bottom = int(y)
        top    = int(y + CELL_H)
        codes.append(kokudo2500_code(left, top))
        geoms.append(box(left, bottom, right, top))
        x += CELL_W
    y += CELL_H

out_name = f"kokudo2500_epsg{epsg}.shp"
gdf = gpd.GeoDataFrame({"code2500": codes}, geometry=geoms, crs=CRS_PLANE)
gdf.to_file(out_name, driver="ESRI Shapefile", encoding="utf-8")

print(f"作成図郭数: {len(codes)}")
print(f"出力完了: {out_name}")

※範囲を系全域より狭くしたい場合は、span や min_lat / max_lat を適宜調整してください。

6. QGISでの確認

1. QGIS でプロジェクトCRSを、選択した系の EPSG に設定
2. 出力した kokudo2500_epsgXXXX.shp を読み込み
3. レイヤのラベルに code2500 を指定
4. 既存の 1/2500 図郭や 1/1250 図郭と重ねてコード一致を確認

同じセル位置でコードが一致していれば、図郭生成は成功です。

※本記事内の実装は執筆時点の仕様に基づくものであり、正式な定義については公共測量標準図式および国土地理院の最新資料をご確認ください。