IU Tips

# ALL# ARCHICAD# BIM# GIS# Plateau# Python# Revit# インタビュー# お知らせ# ワークショップ# 教育# 社員の声# 開発

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

2024.10.21

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

前回、前々回と国土交通省の国土数値情報という日本の地理データを活用する方法を書いてきました。今回はGeoJSON編です。GeoJSONのデータを使って、地理データと緯度経度から用途地域を調べるツールを作成したいと思います。

前回の記事はこちらです。

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

GeoJSONとは

GeoJSONとは地理データを記述するためのデータ記述形式で、地理情報を扱うためのデータ構造を定義したものになります。JSONというのは地理以外にも使われる一般的なデータの記述方法で、JavaScript Object Notationの略です。JavaScriptというプログラミング言語のオブジェクトというデータの表記法に倣った書き方ということを意味します。GeoJSONはJSONを元に作られた地理情報に特化したデータ形式ということになります。

データの取得

今回は用途地域のデータを使いたいと思いますが、用途地域は前回使用した都市計画決定情報データの中に入っています。こちらのデータは下記からダウンロードすることができます。

https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A55-2022.html

ライブラリのインストール

GeoJSONの読み込みはpythonのビルトインライブラリであるjsonでできるので、読み込みをするだけであればライブラリのインストールは必要ありません。後述しますが、matplotlibを使用しますのでインストールをしてください。matplotlibはグラフの描画などに使われるライブラリです。

pip install matplotlib

前回同様、ezdxfとpyprojは必要になりますので、未インストールの場合はインストールしてください。

GeoJSONファイルの読み込み

GeoJSONファイルの読み込みをしていきます。
先ほども書いたようにGeoJSONファイルの読み込みはjsonライブラリでできるので、通常のjsonと同じ方法です。読み込み後はpythonの辞書型として扱えます。

import json

def read_geojson(file_path: str) -> dict:
    with open(file_path, "r", encoding="utf-8") as f:
        data = json.load(f)

    return data

file_path = "27100_youto.geojson"
data = read_geojson(file_path)

トランスフォーマーの作成

緯度経度情報を平面直角座標に変換するのにpyprojのトランスフォーマーを作成します。前回も同様の作業をしているので詳しい説明は省きます。詳しくは第一回をご確認ください。

国土数値情報を使ってCADデータを作成する【Python】【DXF】【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

transformer = create_transformer()

指定地点の座標変換

用途地域の形状は後でまとめて処理するため、確認したい地点の座標を平面直角座標に変換しておきます。緯度経度をトランスフォーマーで平面直角座標に変換します。

lat = 34.691221
lon = 135.501892

p = transformer.transform(lon, lat)

図面データの準備

用途地域はいろいろあるため、レイヤーを分けて色分けして表示したいです。そのため中身を処理する前に用途地域に対応するレイヤーを作っておきます。

ezdxfの初期設定をしつつレイヤーを作成し、作成レイヤーの色を指定します。レイヤーの色は前回説明したAutoCAD Color Indexの番号です。

def create_layers(doc: ezdxf.document.Drawing) -> None:
    layers = [
        ("境界線", 0),
        ("第１種低層住居専用地域", 11),
        ("第２種低層住居専用地域", 31),
        ("第１種中高層住居専用地域", 51),
        ("第２種中高層住居専用地域", 71),
        ("第１種住居地域", 91),
        ("第２種住居地域", 111),
        ("準住居地域", 151),
        ("田園住居地域", 171),
        ("近隣商業地域", 211),
        ("商業地域", 231),
        ("準工業地域", 20),
        ("工業地域", 40),
        ("工業専用地域", 60),
    ]

    for name, color in layers:
        layer = doc.layers.add(name)
        layer.color = color
 
doc = ezdxf.new()
doc.styles.add("MSゴシック", font="msgothic.ttc")
msp = doc.modelspace()
create_layers(doc)

doc.stylesというのはテキストに適用する文字のスタイルです。今回は日本語を使用するので日本語用のフォントを指定しています。

ここでレイヤーを作っておけば、用途地域の処理時にレイヤーを指定することでレイヤーを分けることができます。

GeoJSONの内容確認

それでは、データの抽出をする前に、GeoJSONの中身を確認していきます。

データを見るとfeaturesの中に各用途の情報が入っていて、propertiesの中に容積率、建ぺい率のような情報、geometryの中に座標情報が入っているようです。

{
    "type": "FeatureCollection",
    "name": "27100_youto",
    "crs": {
        "type": "name",
        "properties": {
            "name": "urn:ogc:def:crs:EPSG::6668"
        }
    },
    "features": [
        {
            "type": "Feature",
            "properties": {
                "用途地域": "第２種中高層住居専用地域",
                "YoutoID": 4,
                "容積率": 200.0,
                "建ぺい率": 60.0,
                "Pref": "大阪府",
                "Citycode": "27100",
                "Cityname": "大阪市",
                "決定日": null,
                "決定区分": null,
                "決定者": null,
                "告示番号": null
            },
            "geometry": {
                "type": "Polygon",
                "coordinates": [
                    [
                        [
                            135.546525583,
                            34.766143877
                        ],
                        [
                            135.546495604,
                            34.766132497
                        ],
                        [
                            135.546473805,
                            34.766123401
                        ],
			　　　　　　　......
                    ]
                ]
            }
        },

これを取り出して利用していきます。

GeoJSONの処理

GeoJSON内の各用途の情報が入っているfeaturesをループ処理していきます。
まずは各要素から属性情報を情報を取り出します。
その後にポリゴンの情報を取り出し、
　・平面直角座標に変換
　・用途地域を作図
　・指定地点が用途地域の内部にあるか判定
　・情報をテキストで表示
という処理をしていきます。

これが全体の流れになります。

for item in data["features"]:
    use_district = item["properties"]["用途地域"]
    floor_area_ratio = item["properties"]["容積率"]
    building_coverage_ratio = item["properties"]["建ぺい率"]

    polygons = item["geometry"]["coordinates"]

    for polygon in polygons:
        coords = convert_points(polygon, transformer)

        draw_use_district(msp, coords, use_district)

        if is_inside_polyline(p, coords):
            draw_text(
                msp,
                lat,
                lon,
                p,
                use_district,
                floor_area_ratio,
                building_coverage_ratio,
            )

用途地域の作図

用途地域を作図する処理は下記の通りです。
まず境界線を作図、次に塗りつぶした形状を作図します。

def draw_use_district(
    msp: ezdxf.layouts.layout.Modelspace, coords: tuple[float, float], use_district: str
) -> None:
    msp.add_lwpolyline(coords, dxfattribs={"layer": "境界線"})

    hatch = msp.add_hatch(
        color=ezdxf.enums.ACI.BYLAYER, dxfattribs={"layer": use_district}
    )
    hatch.paths.add_polyline_path(coords, is_closed=True)

各データの頂点に対し、境界線は境界線レイヤーにハッチングは用途別に作成します。
塗りつぶしの色は個別ではなくレイヤーの色を使いたいため”BYLAYER”に設定します。

用途地域情報の表示

指定地点が用途地域にある場合その情報をテキストとして表示します。そのためにまず地点の座標がポリラインの中にあるかを判定する必要があります。

今回はmatplotlibライブラリを使います。Matplotlibはグラフなどの図形を描画するのによく使われますが、こんな機能もあるんですね。自力で判定しようとするとベクトルの外積がどうとかややこしいのですが、簡単に判定できてしまいます。

判定ができたら、その地点に○印と情報を表示するテキストを追加します。

from matplotlib.path import Path

def is_inside_polyline(p, coords) -> bool:
    path = Path(coords)
    return path.contains_point(p)

def draw_text(
    msp: ezdxf.layouts.layout.Modelspace,
    lat: float,
    lon: float,
    p: tuple[float, float],
    use_district: str,
    floor_area_ratio: str,
    building_coverage_ratio: str,
) -> None:
    msp.add_circle((p[0], p[1], 0), 3)

    text = msp.add_text(
        f"{lat},{lon} 用途地域:{use_district} 容積率:{floor_area_ratio} 建ぺい率:{building_coverage_ratio}",
        height=300,
        dxfattribs={"style": "MSゴシック"},
    )
    text.set_placement(p)

今回は情報として、緯度、経度、用途地域、容積率、建ぺい率をテキストに表示します。

ファイルの保存

処理が終わったらdxfファイルを保存します。

save_file_path = "youto.dxf"
doc.saveas(save_file_path)

DXFの確認

出来上がったDXFファイルを確認します。
Archicadで開いてみると、色分け表示された用途地域と情報を見ることができます。

小さいですが、丸の位置が緯度経度を指定した地点です。

まとめ

以上のような方法で用途地域の情報をDXFファイルに表示することができました。図面が必要なければ、情報だけ表示するなどの方法も可能ですね。

注意としては、このデータの問題なのか処理の問題なのかポリゴンが重なっている部分があったため、それに関しては使用する前に確認する必要があると思いました。

今回のコードでは穴あきの形状には対応できないので、そのあたりの改善などが必要になってくるかと思います。

また、敷地によっては用途地域をまたいでいる場合もありますが、今回のコードは1点での調査なので用途地域をまたぐケースは確認できません。国土数値情報のサイトにも書いてありますが、こちらのデータは現状あくまで参考のものです。正しい情報は自治体に確認するのが良いでしょう。

このような情報がWeb APIですぐに確認できるようになるといいですね。
今回は以上です。最後までお読みいただきありがとうございました。

コード全文を以下に載せておきます。

import json
import ezdxf
import pyproj
from matplotlib.path import Path


def read_geojson(file_path: str) -> dict:
    with open(file_path, "r", encoding="utf-8") as f:
        data = json.load(f)

    return data


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(
    polygon: list[tuple[float, float]], transformer: pyproj.Transformer
) -> list[tuple[float, float]]:
    coords = list(map(lambda p: transformer.transform(p[0], p[1]), polygon))

    return coords


def create_layers(doc: ezdxf.document.Drawing) -> None:
    layers = [
        ("境界線", 7),
        ("第１種低層住居専用地域", 11),
        ("第２種低層住居専用地域", 31),
        ("第１種中高層住居専用地域", 51),
        ("第２種中高層住居専用地域", 71),
        ("第１種住居地域", 91),
        ("第２種住居地域", 111),
        ("準住居地域", 151),
        ("田園住居地域", 171),
        ("近隣商業地域", 211),
        ("商業地域", 231),
        ("準工業地域", 20),
        ("工業地域", 40),
        ("工業専用地域", 60),
    ]

    for name, color in layers:
        layer = doc.layers.add(name)
        layer.color = color


def draw_use_district(
    msp: ezdxf.layouts.layout.Modelspace, coords: tuple[float, float], use_district: str
) -> None:
    msp.add_lwpolyline(coords, dxfattribs={"layer": "境界線"})

    hatch = msp.add_hatch(
        color=ezdxf.enums.ACI.BYLAYER, dxfattribs={"layer": use_district}
    )
    hatch.paths.add_polyline_path(coords, is_closed=True)


def is_inside_polyline(p, coords) -> bool:
    path = Path(coords)
    return path.contains_point(p)


def draw_text(
    msp: ezdxf.layouts.layout.Modelspace,
    lat: float,
    lon: float,
    p: tuple[float, float],
    use_district: str,
    floor_area_ratio: str,
    building_coverage_ratio: str,
) -> None:
    msp.add_circle((p[0], p[1], 0), 3)

    text = msp.add_text(
        f"{lat},{lon} 用途地域:{use_district} 容積率:{floor_area_ratio} 建ぺい率:{building_coverage_ratio}",
        height=300,
        dxfattribs={"style": "MSゴシック"},
    )
    text.set_placement(p)


def main():
    file_path = "27100_youto.geojson"
    data = read_geojson(file_path)

    transformer = create_transformer()

    doc = ezdxf.new()
    doc.styles.add("MSゴシック", font="msgothic.ttc")
    msp = doc.modelspace()
    create_layers(doc)

    lat = 34.691221
    lon = 135.501892

    p = transformer.transform(lon, lat)

    for item in data["features"]:
        use_district = item["properties"]["用途地域"]
        floor_area_ratio = item["properties"]["容積率"]
        building_coverage_ratio = item["properties"]["建ぺい率"]

        polygons = item["geometry"]["coordinates"]

        for polygon in polygons:
            coords = convert_points(polygon, transformer)

            draw_use_district(msp, coords, use_district)

            if is_inside_polyline(p, coords):
                draw_text(
                    msp,
                    lat,
                    lon,
                    p,
                    use_district,
                    floor_area_ratio,
                    building_coverage_ratio,
                )

    save_file_path = "youto.dxf"
    doc.saveas(save_file_path)


if __name__ == "__main__":
    main()