🗺️

Plateau の 3D データから建物を個別にメッシュとして取り出す

2021/06/17に公開

Python

gis

PLATEAU

tech

Plateau のデータは CityGML と FBX の2つの形式で提供されています．本稿では，前者から OSS のみを使い，建物IDを指定して建物をメッシュとして取り出す方法について解説します．作業の流れは以下のとおりです．

CityGML ファイルを欲しい座標系に投影変換する．
得られた CityGML ファイルを CityJSON に変換する．
得られた CityJSON ファイルから建物データを取り出す．
面を三角形分割する．

今回使用した環境は以下のとおりです．

MacBook Pro (15-inch 2017)
macOS Big Sur 11.3.1
Docker Desktop 3.3.3
CityGML4j
Ubuntu:20.04 に以下のものを pip で入れたイメージ
- jupyterlab
- trimesh
- cjio

今回は MacOS 上で作業しましたが，コードはほとんど Docker 上の Ubuntu で走っています．そろそろ Linux マシンを買うべきかもしれない．

0. 注意点

テクスチャの扱いは今すぐにはできない

今回の方法ではテクスチャがついてこない（cjio に未実装のため）ので，テクスチャの必要なケースでは FBX とメタデータを空間結合などしたほうが早そうです．ほとんどの建物は緯度経度で適当な頂点と最近傍でマッチングすれば足りそうに見えます．

全部書き出すだけなら本稿の方法は不要

FBXを使うか，または cjio をコマンドラインから使って OBJ に直接変換することができます． OBJ にする場合，投影座標系をあらかじめ決めることができますが，テクスチャがついてこないという問題は残ります．CityJSON から OBJ に直接変換するには以下のコマンドを使います．

cjio in.json export --format obj out.obj

他のもので十分かもしれない

数日前に Google Maps に WebGL の API が追加されました．Web 上で動く何かを作るだけならこれで足りる地域やユースケースも多そうです． https://developers.google.com/maps/documentation/javascript/webgl

さらに CityGML は QGIS で直接読めるので GIS 用途ではほとんどの場合 QGIS を使ったほうがよいと思われます．

1. CityGML ファイルの投影変換

Plateau のデータは緯度経度に高さ方向をメートルで表した状態 (EPSG:6697) で入っているので，適当な投影座標系に直します．ここでは平面直角座標系 IX 系 (EPSG:6677) にします．変換ツールが提供されているのでありがたく使います．

docker run --rm -u 1000 \
-v データのあるディレクトリ:/data \
citygml4j/citygml-tools reproject --target-crs=6677 /data

完了すると reprojected と付いたファイルが同じ場所に生成されているので，移動します．

cd データのあるディレクトリ
mkdir -p ../bldg_reprojected
find *reprojected.gml | xargs -I{} mv {} ../bldg_reprojected/{}

2. CityGML to CityJSON

この変換プログラムはデフォルトで数値を小数点以下3桁で切り捨ててしまいますが，これでは 0.001° なので日本付近では約 100m 単位で丸められてしまいます．行政基本情報データ連携モデル標準ガイドライン群 ID:1015-4によれば，測量に使う電子国土基本図は小数点以下9桁の精度でデータを保持しているようなので，本稿でもそれに倣います．実際は6桁でも足りるでしょう．なお，再投影した場合はほとんどのケースでメートル単位でしょうから，桁数はデフォルトの3桁で足りるはずです．

docker run --rm -u 1000 \
-v 再投影されたデータのあるディレクトリ:/data \
citygml4j/citygml-tools to-cityjson --vertices-digits=9 --template-digits=9 /data

完了すると .json の拡張子をもつファイルが同じ場所に生成されているので，移動します．

cd 再投影されたデータのあるディレクトリ
mkdir -p ../bldg_cityjson
find *.json | xargs -I{} mv {} ../bldg_cityjson/{}

3. CityJSON ファイルから建物を取り出す

CityJSON を読み込みます．

import cjio
from cjio import cityjson

# Load CityJSON to dictionary
cj = cityjson.load('ファイルパス')

# get objects and building ids
cityobjects = cj.get_cityobjects(type='building')
building_ids = list(cityobjects.keys())

あとは cityobjects[building_ids[0]] のようにすれば CityObject を取り出せます．Unicode エスケープがそのまま見えているのでやや見苦しいですが，以下のような結果が返ってくるはずです．

CityObject の中身

{
  "id": "BLD_6f03fa39-e6bb-4970-b7f9-d367006f8d0c",
  "type": "Building",
  "attributes": {
    "measuredHeight": 7.7,
    "\u5efa\u7269ID": "13115-bldg-7748",
    "13_\u533a\u5e02\u753a\u6751\u30b3\u30fc\u30c9_\u5927\u5b57\u30fb\u753a\u30b3\u30fc\u30c9_\u753a\u30fb\u4e01\u76ee\u30b3\u30fc\u30c9": "13115008003",
    "\u753a\u30fb\u4e01\u76ee\u30b3\u30fc\u30c9": "3",
    "\u5927\u5b57\u30fb\u753a\u30b3\u30fc\u30c9": "8"
  },
  "children": [],
  "parents": [],
  "geometry_type": [
    "MultiSurface",
    "Solid"
  ],
  "geometry_lod": [
    0,
    1
  ],
  "semantic_surfaces": []
}

4. CityJSON から建物データを取り出す

本来は cjio に付いていてほしい機能ですが見つからなかったので，探し方が悪いか本当に存在しないかどちらかです．以下はとりあえず cjio のリポジトリからコピーしてきた三角形分割のコードに手を加えたものです．

三角形分割とメッシュ構築のコード

以下は調査中に残したメモ（推測含む）です．

vnp : v 頂点座標を np numpy.ndarray にしたもの，ということのようです
face: 頂点番号の配列の配列です．最初の配列が外側の境界，2番目以降が穴のようです
Javascript 版の earcut は3次元の三角形分割もできますが，これの C++ へのポートの Python バインディングである mapbox_earcut は2次元の分割のみ行えます．このため，一旦法線を求めてそれに垂直な面に投影してから三角形分割をしているようです

import numpy as np
import mapbox_earcut

# Taken from 
# A modification in to_2d was needed

def my_triangulate_face(face, vnp):
    if ((len(face) == 1) and (len(face[0]) == 3)):
        # print ("Already a triangle")
        return face
    sf = np.array([], dtype=np.int32)
    for ring in face:
        sf = np.hstack((sf, np.array(ring)))
    sfv = vnp[sf]
    rings = np.zeros(len(face), dtype=np.int32)
    total = 0
    for i in range(len(face)):
        total += len(face[i])
        rings[i] = total
        # 1. normal with Newell's method
    n = cjio.geom_help.get_normal_newell(sfv)
    sfv2d = np.zeros((sfv.shape[0], 2))
    for i, p in enumerate(sfv):
        xy = cjio.geom_help.to_2d(p, n[0]) # FIXED # TODO CityJSON to OBJ の変換では問題なく動いているようだが未調査
        sfv2d[i][0] = xy[0]
        sfv2d[i][1] = xy[1]
    result = mapbox_earcut.triangulate_float64(sfv2d, rings) # FIXED float64のほうを使わないと緯度経度でデータを持っているときに時々失敗する

    for i, each in enumerate(result):
        result[i] = int(sf[each])

    #    print (type(result.reshape(-1, 3).tolist()[0][0]))
    return result.reshape(-1, 3).tolist()

def cityobject_to_trimesh(city_object, lod=1):
    geom_lod_filtered = list(filter(lambda x : x.lod == lod, city_object.geometry))
    
    result = []
    
    for geometry in geom_lod_filtered:
    
        face_by_indices = geometry.build_index()
        boundaries = face_by_indices[0]
        vertices = np.array(list(face_by_indices[1].keys()))
        
        # Triangulate
        trifaces = np.empty((0, 3), dtype=np.int64)
        for shell in boundaries:
            for face in shell:
                tris = np.array(my_triangulate_face(face, vertices))
                trifaces = np.vstack((trifaces, tris))

        mesh = trimesh.Trimesh(vertices=vertices, faces=trifaces)
        result.append(mesh)

    return result

あとは，この関数を使って以下のように CityObject をメッシュに変換して STL などに書き出すことができます．

cityobject = cityobjects[building_ids[0]]
mesh = cityobject_to_trimesh(cityobject)[0]
mesh.export('mesh.stl');

5. 結果

先頭から 100 個くらいを書き出してみました．