対象読者
「国土基本図の図郭」をシェープファイルで作成したい人
測量、地図製作、GISなどに携わる方々
pythonユーザ
実行環境
Windows10
Python 3.9.2
pyshp 2.1.3
はじめに
国土地理院の基盤地図情報や統計局のe-Statなどのオープンデータは、２次メッシュや３次メッシュといった地域メッシュというメッシュ体系で整理されていて、図郭もGISデータとして入手しやすくなっています1。
対して、市町村の地形図やLPデータなど公共測量によって作られたデータは、国土基本図の図郭（旧、国土基本図図郭）という図郭で図郭割りされて整理されていて、図郭のデータは一般に流通していないようです。
国土基本図の図郭の図郭割りの方法は公共測量標準図式に明記されているので、今回はこれに従って国土基本図の図郭をpythonで作成します。
出力ファイルの仕様は以下の通りです。
データ形式： シェープファイル（ポリゴン）
図郭の大きさ： 地図情報レベル（以降、レベル）2500
範囲： 40km×30kmの区画
また、シェープファイルを読み書きできるライブラリpyshpを使用します。
国土基本図の図郭について
作業規程の準則の付録7の公共測量標準図式の第84条に定められた図郭割りで切られた図郭のことです。
令和2年3月31日の改定で、なぜか「国土基本図図郭」から「国土基本図 の 図郭」に改名されました。
同第83条には座標軸についての定めがあり、数値地形図データファイルは測量座標系、写真地図データファイルは数学座標系となっています。
今回は横軸をx、縦軸をyとしてGISソフトで扱いたいので、ここでは数学座標系を前提にします。
また、原点は平面直角座標系の原点(0,0)、単位はmとします。
なお、コード番号を付与するために平面直角座標系の系番号を指定しますが、シェープファイルの座標系の定義はここでは行いません。（座標系定義は別途GISソフトで指定することを想定します）
分割法について
公共測量標準図式の第84条の4の一、二、三を引用します。
一
区画名は、各座標系のＹ軸及びＸ軸を基準とし、南北300km、東西160km を含む区域を30km×40kmの長方形に分割して区画を定め、下図によりアルファベット大文字の組合せで表示する。
二
地図情報レベル5000 にあっては座標系内の１区画を100 等分し、下図によりアラビア数字で表示する。
三
地図情報レベル2500 にあっては、地図情報レベル5000 の図郭に相当する区画を各辺で２等分して得られる４個の区画に北西側、北東側、南西側、南東側の順に１～４のアラビア数字で区画番号を定め、地図情報レベル5000 の図郭番号に追加する。
実際は、四にレベル1000、五にレベル500、六に令和2年の改正で追加されたレベル250の分割法についての記載が続きますが、今回はレベル2500の図郭を作るので省略します。
図郭サイズをまとめると下表のとおりです。
なお、図郭名は今回便宜的につけたものです。
| 図郭名 | 東西幅 | 南北幅 |
|---|---|---|
| 区域 | 320km | 600km |
| 区画 | 40km | 30km |
| レベル5000 | 4000m | 3000m |
| レベル2500 | 2000m | 1500m |
pyshpのインストール
本記事投稿時点のバージョンは2.1.3になります。
pipなら
$ pip install pyshp
condaなら
$ conda install pyshp
区画の北西端の座標を取得する関数
図郭番号を振るとき、区画の北西端をスタートとして考えると分かりやすいので、その頂点座標を取得する関数を定義します。
引数に区画のアルファベット２文字を指定します。（実際は１文字目がA~T、２文字目がA~Hに限定されるため、そこをチェックした方が親切です）
戻り値 org_x, org_y が指定の区画の北西端のx座標とy座標です。
def get_org_xy(kukaku: str) -> tuple[int, int]:
str_x = kukaku[1]
str_y = kukaku[0]
shift_x = ord(str_x) - 65
shift_y = ord(str_y) - 65
org_x = (-160 + shift_x*40) * 1000
org_y = (300 - shift_y*30) * 1000
return org_x, org_y
shift_x, shift_y は、区画の上位区域の西端北端それぞれから何区画分シフトさせるかという変数で、ord関数を使って取得しています。
例えば、区画'CD'の北西端座標は、区域の北西端座標(-160km, 300km)から東へ３区画分、南へ２区画分シフトさせた座標になります。'A'のアスキーコードが65です。
str_x = 'D'
str_y = 'C'
shift_x = ord('D') - 65 = 68 - 65 = 3
shift_y = ord('C') - 65 = 67 - 65 = 2
org_x = (-160 + 3*40) * 1000 = -40000
org_y = (300 - 2*30) * 1000 = 24000
1図郭を作成する関数
シェープファイルにポリゴンを作成する関数です。
def make_mesh(w: shapefile.Writer, code: str,
xw: int, xe: int, yn: int, ys: int) -> None:
# ジオメトリの書込
w.poly([
[[xw,yn],[xe,yn],[xe,ys],[xw,ys],[xw,yn]]
])
# レコードの書込
w.record(code)
引数の説明です。
w: shapefile.Writerクラスのインスタンス
code: 図郭番号
xw, xe, yn, ys: 図郭の西端x座標、東端x座標、北端y座標、南端y座標
メインコード
アルファベット２文字の区画番号、平面直角座標系の系番号、ファイルの出力先を指定して、シェープファイルを出力します。
def make_zukaku(kukaku: str, kei: int, fold: str) -> None:
# 地図情報レベル
level = 2500
# 出力ファイル名
filename = f'Level{level:04}{kukaku}.shp'
# 図郭番号用フィールド名
fieldname = f'code{level:04}'
# 区画の北西端座標
org_x, org_y = get_org_xy(kukaku)
# レベル5000の図郭の東西幅と南北幅[m]
dx5000 = int(40/10 * 1000) # 4000
dy5000 = int(30/10 * 1000) # 3000
# レベル2500の図郭の東西幅と南北幅[m]
dx = int(dx5000 / (5000/level)) # 2000
dy = int(dy5000 / (5000/level)) # 1500
# ポリゴンタイプのシェープファイルを新規作成
with shapefile.Writer(os.path.join(fold, filename), shapetype=5) as w:
# 図郭番号格納用のフィールドを追加
w.field(fieldname, 'C', size=9)
# １区画を南北10分割、東西10分割
for y in range(10):
for x in range(10):
# レベル5000の図郭の北西端座標
org_5000_x = org_x + x*dx5000
org_5000_y = org_y - y*dy5000
# レベル5000の図郭を南北2分割、東西2分割
# 北西、北東、南西、南東の順に1,2,3,4で番号付け
for i in range(1,5):
# レベル2500の図郭の西端と東端の座標
if i in (1,3):
xw = org_5000_x
xe = org_5000_x + dx
else:
xw = org_5000_x + dx
xe = org_5000_x + dx*2
# レベル2500の図郭の北端と南端の座標
if i in (1,2):
yn = org_5000_y
ys = org_5000_y - dy
else:
yn = org_5000_y - dy
ys = org_5000_y - dy*2
# 図郭番号
code = f'{kei:02}{kukaku}{y}{x}{i}'
# データ化
make_mesh(w, code, xw, xe, yn, ys)
実行
検証用データとして兵庫県_全域_標高ラスター 2を使いたいので、平面直角座標系第5系で区画'OF'の2500レベルの図郭を作成します。
区画'OF'は、明石海峡と中国自動車道の間あたりの範囲になります。
`インポートと関数のコードはこちら`
import os.path
import shapefile
def get_org_xy(kukaku: str) -> tuple[int, int]:
str_x = kukaku[1]
str_y = kukaku[0]
shift_x = ord(str_x) - 65
shift_y = ord(str_y) - 65
org_x = (-160 + shift_x*40) * 1000
org_y = (300 - shift_y*30) * 1000
return org_x, org_y
def make_mesh(w: shapefile.Writer, code: str,
xw: int, xe: int, yn: int, ys: int) -> None:
w.poly([
[[xw,yn],[xe,yn],[xe,ys],[xw,ys],[xw,yn]]
])
w.record(code)
def make_zukaku(kukaku: str, kei: int, fold: str) -> None:
level = 2500
filename = f'Level{level:04}{kukaku}.shp'
fieldname = f'code{level:04}'
org_x, org_y = get_org_xy(kukaku)
dx5000 = int(40/10 * 1000)
dy5000 = int(30/10 * 1000)
dx = int(dx5000 / (5000/level))
dy = int(dy5000 / (5000/level))
with shapefile.Writer(os.path.join(fold, filename), shapetype=5) as w:
w.field(fieldname, 'C', size=9)
for y in range(10):
for x in range(10):
org_5000_x = org_x + x*dx5000
org_5000_y = org_y - y*dy5000
for i in range(1,5):
if i in (1,3):
xw = org_5000_x
xe = org_5000_x + dx
else:
xw = org_5000_x + dx
xe = org_5000_x + dx*2
if i in (1,2):
yn = org_5000_y
ys = org_5000_y - dy
else:
yn = org_5000_y - dy
ys = org_5000_y - dy*2
code = f'{kei:02}{kukaku}{y}{x}{i}'
make_mesh(w, code, xw, xe, yn, ys)
以下の引数で実行します。
kukaku = 'OF'
kei = 5
fold = r'C:\test'
make_zukaku(kukaku, kei, fold)
出力ファイル名は Level2500OF.shp になります。
結果
標高値をグレースケールで表した主題図（ダウンロードデータ）と、境界と図郭名を赤で示した索引図（作成した図郭）をGISソフトで重ねています。
検証用の標高ラスターは図郭番号05OF891、平面解像度1mのTIFFファイルで、ファイル名は alt_05OF891.tif です。
両ファイルに世界測地系平面直角座標系第5系を定義し、マップ表示も同じ設定しています3。
図の通り、標高ラスターと作成した図郭の境界が一致していることが確認できました。
最後に
40km×30kmの区画からレベル2500の図郭を作成するコードを書きました。
今回は簡略化のためレベル2500のみを対象にしましたが、レベル1000, 500, 250も同様の考えでコード化できます。レベル5000の図郭を基準に分割する点は同じで、図郭番号の決め方がより規則的になります。
また、シェープファイルの作成に使用したpyshpですが、arcpyのデータアクセスモジュールより簡単にシェープファイルを作れるのでおすすめです。
なお、座標値や図郭の幅は単位mで扱う限り整数値しかとらないはずですが、除算(/)するとfloat型になってしまうので明示的にint型としました。除算(//)を使ってint型を維持してもよいです。小数点以下切り捨ての意図はありません。