国土数値情報を使ってCADデータを作成する【Python】【DXF】【Shapefile】

2024.09.03

こんにちは。
IU BIM STUDIOの原田です。

今回は国土交通省が提供する日本の地形等に関するデータ、国土数値情報を使ってみたいと思います。

国土数値情報とは

国土数値情報とは国土交通省が行政区域、河川などの情報をGISデータで整備しているものです。
無料でダウンロードして使うことができます。
ダウンロードは下記のリンクからすることができます。

https://nlftp.mlit.go.jp

データのダウンロード

上記のサイトにアクセスしたら、画像の赤枠の部分、国土数値情報をクリックします。

データのダウンロードページに遷移しますので、今回は行政区域を選択します。

ダウンロードしたい都道府県を選びます。
今回は大阪府を選びます。

一覧に飛びますのでダウンロードボタンを押します。

ダウンロードできたら、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()