IU Tips

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

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

2024.09.30

こんにちは。
IU BIM STUDIOの原田です。
今回も国土数値情報を使ってみたいと思います。

前回は大阪府の形状をDXF形式に変換し保存しました。

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

今回は大阪府の公園の形状をDXFに出力します。
データ形式は今回はCityGMLを使ってみたいと思います。

CityGMLとは

CityGMLは都市の情報を交換可能な状態で記述するためのデータ標準形式で、GML(Geography Markup Language)という地理情報を記述するデータ形式の一種となっています。GMLはXMLの一種でCityGMLもXMLをベースとしたデータ形式になっています。

データのダウンロード

下記のリンクから国土数値情報のページに飛びます。
https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-A55-2022.html

CityGMLを選び、今回は大阪市を指定して下部のボタンからダウンロードします。

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

Pythonには標準でxmlを扱うためのライブラリがあるのですが、いろいろ調べたところlxmlというライブラリを使った方が取り扱いが楽なようです。

というのも、ダウンロードしたCityGMLには名前空間が設定してあり、名前空間のあるXMLを取り扱うにはlxmlの方が適しているためです。XMLではタグというものを使ってデータの種別などを記述しますが、名前空間はデータが複雑になったときにタグの名称が重複しないように整理するためのものです

pip install lxml

前回と同じく ezdxfとpyprojを使います。必要に応じてインストールしてください。

ライブラリを読み込む

使用するライブラリを読み込んでおきます。

from lxml import etree
import ezdxf
import pyproj

CityGMLファイルの読み込み

CityGMLファイルはlxmlで読み込みます。

file_path = "27100_kouen.gml"
tree = etree.parse(file_path)
root = tree.getroot()

これで公園のデータを読み込むことができます。
Pythonデフォルトのxmlライブラリだと名前空間の設定が必要だったりするのですが、lxmlを使うと名前空間はうまいことやってくれているようです。

データ構造の確認

公園のデータを読み込むため元データの構造を見てみます。

<?xml version="1.0" encoding="UTF-8"?>
<core:CityModel xmlns:urf="https://www.geospatial.jp/iur/urf/2.0" xmlns:core="http://www.opengis.net/citygml/2.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:gml="http://www.opengis.net/gml" xmlns:xlink="http://www.w3.org/1999/xlink" xsi:schemaLocation="https://www.geospatial.jp/iur/urf/2.0 https://www.geospatial.jp/iur/schemas/urf/2.0/urbanFunction.xsd http://www.opengis.net/citygml/2.0 http://schemas.opengis.net/citygml/2.0/cityGMLBase.xsd http://www.opengis.net/gml http://schemas.opengis.net/gml/3.1.1/base/gml.xsd">
  <gml:boundedBy>
    <gml:Envelope srsName="http://www.opengis.net/def/crs/EPSG/0/6697" srsDimension="3">
      <gml:lowerCorner>34.58729966426313 135.41738406618532 0</gml:lowerCorner>
      <gml:upperCorner>34.76613152642622 135.5982481871067 0</gml:upperCorner>
    </gml:Envelope>
  </gml:boundedBy>
  <core:cityObjectMember>
    <urf:UrbanFacility gml:id="0243a6a2-c222-45bb-bdc9-c9287c5d7baa">
      <urf:function codeSpace="https://www.geospatial.jp/iur/codelists/2.0/Common_urbanFacilityType.xml">3021</urf:function>
      <urf:prefecture codeSpace="https://www.geospatial.jp/iur/codelists/2.0/Common_prefecture.xml">27</urf:prefecture>
      <urf:city codeSpace="https://www.geospatial.jp/iur/codelists/2.0/Common_localPublicAuthorities.xml">27100</urf:city>
      <urf:validFrom></urf:validFrom>
      <urf:validFromType codeSpace="https://www.geospatial.jp/iur/codelists/2.0/Common_validType.xml"></urf:validFromType>
      <urf:custodian></urf:custodian>
      <urf:notificationNumber></urf:notificationNumber>
      <urf:lod1MultiSurface>
        <gml:MultiSurface srsName="http://www.opengis.net/def/crs/EPSG/0/6697">
          <gml:surfaceMember>
            <gml:Polygon>
              <gml:exterior>
                <gml:LinearRing>
                  <gml:posList>34.764883424688776 135.5442919973323 0 34.76494279783606 135.5441268264767 0 34.76509040414205 135.5437343710176 0 34.76530774271833 135.5431738894819 0 34.76545853435434 135.54331324202607 0 34.76567278730101 135.5435201692168 0 34.76576256693442 135.54360927650146 0 34.76613152642622 135.54396590525533 0 34.766070310475314 135.54409955847262 0 34.76576068965223 135.54470062690004 0 34.76547369507133 135.5444946381494 0 34.765375258864815 135.54468907539203 0 34.76536850781795 135.54469184320428 0 34.76523291896809 135.5445942677456 0 34.76515383578045 135.5445400799259 0 34.76517172216408 135.54450174839917 0 34.764883424688776 135.5442919973323 0</gml:posList>
                </gml:LinearRing>
              </gml:exterior>
            </gml:Polygon>
          </gml:surfaceMember>
        </gml:MultiSurface>
      </urf:lod1MultiSurface>
    </urf:UrbanFacility>
  </core:cityObjectMember>
  <core:cityObjectMember>

...

見たところ<core:cityObjectMember>が公園一つ分のデータで、その奥深くにposListという緯度経度のデータが入っているようです。
これを取り出したいと思います。

まずは<core:cityObjectMember>を取り出してループ処理をしていきたいと思います。データを取り出す方法はいろいろありますが、今回はXpathを使います。

items = root.xpath("//gml:posList", namespaces=root.nsmap)

XpathはXMLやHTMLの要素を指定するための記述法で、意図したデータを効率的に取り出すために使います。上記の書き方はroot以下の全ての<gml:posList>を検索して取得するということになります。

緯度経度の取り出し

それでは取り出した公園データの緯度経度を取り出していきます。

def get_coord_list(items: list[etree._Element]) -> list[list[tuple[float, float]]]:
    coord_list = []
    for item in items:
        pos_list = item.text.split(" ")
        pos_float = list(map(float, pos_list))
        coords = [(pos_float[i + 1], pos_float[i]) for i in range(0, len(pos_float), 3)]
        coord_list.append(coords)
    return coord_list

coord_list = get_coord_list(items)

上記のCityGMLデータのように、posListは緯度、経度、高度がスペース区切りにされたものが連続した文字列になっています。それをスペースで分割、数値に変換、高度を除いた緯度経度のみを取り出します。それから公園ごとに緯度経度をリストにしてまとめます。

緯度経度から平面直交座標への変換

現在の緯度経度ではCAD上で正しい形状を表せないため、座標系を変換します。こちらは前回と同じなので、詳しい説明は省きます。

def create_transformer() -> pyproj.Transformer:
    src_proj = "EPSG:6697"
    dest_proj = "EPSG:6674"
    transformer = pyproj.Transformer.from_crs(src_proj, dest_proj, always_xy=True)
    return transformer

def convert_points(
    coord_list: list[list[tuple[float, float]]], transformer: pyproj.Transformer
) -> list[list[tuple[float, float]]]:
    polylines = []

    for coords in coord_list:
        points_transfer = list(map(lambda p: transformer.transform(p[0], p[1]), coords))
        polylines.append(points_transfer)

    return polylines

transformer = create_transformer()
polyline_vertices = convert_points(coord_list, transformer)

公園の範囲を単色の塗りつぶしで表示

前回は線のみの表示でしたが、今回は公園の範囲を緑色で塗りつぶしてみたいと思います。

def create_dxf(polyline_vertices: list[tuple[float, float]]):
    doc = ezdxf.new()
    msp = doc.modelspace()

    for vertices in polyline_vertices:
        msp.add_lwpolyline(vertices)
        hatch = msp.add_hatch(color=81)
        hatch.paths.add_polyline_path(vertices)

    doc.saveas("polylines.dxf")

create_dxf(polyline_vertices)

単色で塗りつぶすにはハッチングを作り、それに輪郭線を追加するという方法になります。
色の指定はRGBでも可能のようですが、今回はAutoCAD Color Indexを使いたいと思います。

AutoCAD Color Indexとは、あらかじめ決められた番号を指定することで256色の範囲でRGB値を入力することなく色を指定できるもののようです。


https://ezdxf.readthedocs.io/en/stable/concepts/aci.html#aci

今回使ったのは81番でリンク先の画像の通り、薄めの緑色です。
DXFを書き出して前回のデータと重ねて確認します。

このように公園が書き出せました。真ん中の大きいのが大阪城公園ですね。
無事、CityGMLの公園データをDXFにすることができました。

今回は以上です。
最後までお読みいただきありがとうございました。

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

from lxml import etree
import ezdxf
import pyproj


def get_xml_root(file_path: str) -> etree._Element:
    tree = etree.parse(file_path)
    return tree.getroot()


def get_coord_list(items: list[etree._Element]) -> list[list[tuple[float, float]]]:
    coord_list = []
    for item in items:
        pos_str = item.text.split(" ")
        pos_float = list(map(float, pos_str))
        coords = [(pos_float[i + 1], pos_float[i]) for i in range(0, len(pos_float), 3)]

        coord_list.append(coords)

    return coord_list


def create_transformer() -> pyproj.Transformer:
    src_proj = "EPSG:6697"
    dest_proj = "EPSG:6674"

    transformer = pyproj.Transformer.from_crs(src_proj, dest_proj, always_xy=True)

    return transformer


def convert_points(
    coord_list: list[list[tuple[float, float]]], transformer: pyproj.Transformer
) -> list[list[tuple[float, float]]]:
    polylines = []

    for coords in coord_list:
        points_transfer = list(map(lambda p: transformer.transform(p[0], p[1]), coords))
        polylines.append(points_transfer)

    return polylines


def create_dxf(polyline_vertices: list[tuple[float, float]]):
    doc = ezdxf.new()
    msp = doc.modelspace()

    for vertices in polyline_vertices:
        msp.add_lwpolyline(vertices)

        hatch = msp.add_hatch(color=81)
        hatch.paths.add_polyline_path(vertices)

    doc.saveas("polylines.dxf")


def main():
    file_path = "27100_kouen.gml"
    root = get_xml_root(file_path)

    items = root.xpath("//gml:posList", namespaces=root.nsmap)

    coord_list = get_coord_list(items)
    transformer = create_transformer()
    polyline_vertices = convert_points(coord_list, transformer)

    create_dxf(polyline_vertices)


if __name__ == "__main__":
    main()