GISPython開発
国土数値情報を使ってCADデータを作成する【Python】【DXF】【Shapefile】
2024.09.03
国土数値情報を使ってCADデータを作成する【Python】【DXF】【Shapefile】
2024.09.03
こんにちは。
IU BIM STUDIOの原田です。
今回は国土交通省が提供する日本の地形等に関するデータ、国土数値情報を使ってみたいと思います。
国土数値情報とは
国土数値情報とは国土交通省が行政区域、河川などの情報をGISデータで整備しているものです。
無料でダウンロードして使うことができます。
ダウンロードは下記のリンクからすることができます。
データのダウンロード
上記のサイトにアクセスしたら、画像の赤枠の部分、国土数値情報をクリックします。
データのダウンロードページに遷移しますので、今回は行政区域を選択します。
ダウンロードしたい都道府県を選びます。
今回は大阪府を選びます。
一覧に飛びますのでダウンロードボタンを押します。
ダウンロードできたら、zipファイルを展開してください。
中身のデータはGML、GeoJSON、Shapefileだそうですが、今回はshapefileを試してみます。
ライブラリのインストール
Shapefileを読み込むためのpythonライブラリはpyshpというのがあるそうです。また、緯度経度と直角平面座標系を変換するのにpyprojというライブラリを使います。
これらをインストールしましょう。
pip install pyshp
pip install pyprojまた、今回は行政区域のデータをDXFファイルに書き出すのですが、DXFをpythonで扱うためのライブラリ、ezdxfもインストールします。
pip install ezdxf準備は以上で完了です。
Shapefileの読み込み
Shapefileは下記のような拡張子のファイルを使うそうです。
N03-20240101_27.dbf
N03-20240101_27.prj
N03-20240101_27.shp
N03-20240101_27.shx
とはいえpyshpを使えば読み込み自体は簡単にできるようなので、あまり深く考えずにやっていきます。
Import時はpyshpではなく、shapefileのようなので気を付けてください。
ファイルを読み込みます。
import shapefile
shp_path = "N03-20240101_27.shp"
sf = shapefile.Reader(shp_path, encoding="utf-8")次のようにして情報を取り出すことができます。
shapes = sf.shapes() # 位置・形状情報
records = sf.records() # 属性情報形状の頂点座標は下記のように取り出します。
for shape in sf.shapes():
for point in shape.points:
print(point)
> (135.514807224, 34.720568892)
> (135.514894449, 34.720739441)
> (135.51531, 34.721553054)
> (135.515490558, 34.721906667)
> ...座標情報は経度、緯度のタプルとして入っていますのでこれを加工していきます。
座標の変換
座標は緯度経度で入っていましたので、そのままCADで使おうとすると正しい形状、大きさになりません。そのため形状を平面直角座標に変換していきます。変換には先ほどインストールしたpyprojを使います。
pyprojではトランスフォーマーというのを作って座標を変換するそうです。そしてトランスフォーマーを作るのにEPSGコードというのを入力します。
JGD2011という2011年に測量されたデータの場合、EPSGコードは緯度経度がEPSG6668、大阪の平面直角座標がEPSG6674ということなのでこれを使います。あまり詳しくないので違っていたらすみません。
import pyproj
src_proj = "EPSG:6668"
dest_proj = "EPSG:6674"
transformer = pyproj.Transformer.from_crs(src_proj, dest_proj, always_xy=True)
トランスフォーマーを使って座標変換します。
vertices_list = []
for shape in sf.shapes():
points_transfer = list(
map(lambda p: transformer.transform(p[0], p[1]), shape.points)
)
vertices_list.append(points_transfer)これで直角平面座標の座標が取得できました。
DXFファイルの作成
次にezdxfを使い座標からポリラインを作成していきます。
import ezdxf
doc = ezdxf.new()
msp = doc.modelspace()
for vertices in vertices_list:
msp.add_lwpolyline(vertices)
doc.saveas("polylines.dxf")これで境界線がポリラインとして保存されました。
では、保存されたdxfを確認します。
何か変な線がありますね…
何なんだこれは…
Googleマップで調べてみます。
どうやらこれは飛び地のようでした…
こういうのがあるから処理が大変ですよね…
データで扱いやすいようにまとめてくれればいいのに…
とは言え、そうもいかないのでコードを修正します。
shapefileのデータを見てみると、各shapeの中のpartsというリストにポリラインの区切りが入っているようでした。
なので、partsの長さが2以上の時はポリラインを分割します。
vertices_list = []
for shape in sf.shapes():
part_count = len(shape.parts)
if part_count >= 2:
for i in range(part_count - 1):
start = shape.parts[i]
end = shape.parts[i + 1]
part = shape.points[start:end]
points_transfer = list(
map(lambda p: transformer.transform(p[0], p[1]), part)
)
vertices_list.append(points_transfer)
else:
points_transfer = list(
map(lambda p: transformer.transform(p[0], p[1]), shape.points)
)
vertices_list.append(points_transfer)再度DXFファイルを確認してみます。
これで大阪府の形状をdxfにすることができました。
今回のまとめ
これで国土数値情報のデータからCADで使えるデータが作成できました。
他にもいろいろなデータがあるので、何が使えそうか調べてみたいと思います。
今回は以上です。
最後までお読みいただきありがとうございました。
コードの全文を以下に載せておきます。
import ezdxf
import pyproj
import shapefile
def create_transformer() -> pyproj.Transformer:
src_proj = "EPSG:6668"
dest_proj = "EPSG:6674"
transformer = pyproj.Transformer.from_crs(src_proj, dest_proj, always_xy=True)
return transformer
def convert_points(
sf: shapefile.Reader, transformer: pyproj.Transformer
) -> list[tuple[float, float]]:
vertices_list = []
for shape in sf.shapes():
part_count = len(shape.parts)
if part_count >= 2:
for i in range(part_count - 1):
start = shape.parts[i]
end = shape.parts[i + 1]
part = shape.points[start:end]
points_transfer = list(
map(lambda p: transformer.transform(p[0], p[1]), part)
)
vertices_list.append(points_transfer)
else:
points_transfer = list(
map(lambda p: transformer.transform(p[0], p[1]), shape.points)
)
vertices_list.append(points_transfer)
return vertices_list
def create_dxf(vertices_list: list[tuple[float, float]]):
doc = ezdxf.new()
msp = doc.modelspace()
for vertices in vertices_list:
msp.add_lwpolyline(vertices)
doc.saveas("polylinestest.dxf")
def main():
shp_path = "N03-20240101_27.shp"
sf = shapefile.Reader(shp_path, encoding="utf-8")
transformer = create_transformer()
vertices_list = convert_points(sf, transformer)
create_dxf(vertices_list)
if __name__ == "__main__":
main()
