ジェネリックPLATEAU VIEWを作りたい -CityGMLからメッシュモデルの作成-
CityGMLというデータ形式の勉強をするついでに、このデータをOpen3Dのメッシュデータ形式(TriangleMesh)に変換し描画する処理を書いたので記録としてまとめます。
ちょっと綺麗になった
— sakamo (@ppyhr) April 12, 2024
高さ30m以上の建物をマゼンタで塗ってます https://t.co/X77wHi355P pic.twitter.com/7Rx9sdsQLH
目次
はじめに
みなさんは3Dモデルがどのように表現されているか知っていますか？
実は、よく見る3Dモデル(CAD, CGなど)は中身の詰まっていない、頂点・面・エッジ(線)で構成されるメッシュモデルとして表現・描画されていることが多いです*。
以下の画像のイメージです。
(https://sketchfab.com/3d-models/stanford-bunny-43f266d6cd6e4c6888b9943557528c0f)
* 中身が詰まっているモデルもあり、これらは「ソリッド」メッシュといいます。メッシュモデルと呼ぶ際には通常「サーフェス」モデルを指します。
メッシュモデルについてはこちらが参考になります。
https://help.autodesk.com/view/CIV3D/2024/JPN/index.html?guid=GUID-A6232957-5039-4AB7-8B1D-8FD0AD98F77B
しかし、PLATEAUが公開しているCityGML形式のデータはメッシュモデルの形式で書かれていません。したがって、PLATEAUデータから都市モデルを作成するためには、gml形式で書かれたデータをメッシュモデルに変換する必要が生じるわけです。
そこで、今回はPLATEAUのCityGMLデータをOpen3Dのメッシュデータに変換し、作成したメッシュモデルを描画する処理を書いてみました。
CityGMLについてはこちら
メッシュモデルの原理が分かると何が嬉しい？
PLATEAUのデータを自分で書き換えることが可能になります。つまり追加・削除が可能となるため、自分で任意の形状のオブジェクトを追加したり、自分で決定した閾値でオブジェクトを色分けしたりできるようになります！
形式への理解があれば、作成したモデルをPLATEAUと全く関係のない分野に持ち込み扱うことも容易になります。
処理の流れ
1. データの読み込み
文京区のCityGMLからGMLファイルを読み込みます。
読み込みにはgeopandasライブラリを使っています。
ここで、PLATEAUの3D都市モデルには
・建物のフットプリント(の頂点)情報
・高さ情報
が含まれており、gmlファイルを読み込んだDataFrame(表形式のデータ構造)を実際に見てみると、geometry列にフットプリントの頂点情報が、measuredHeight列に建物の高さ情報が含まれていることがわかります。
先頭10要素のmeasuredHeightとgeometry(フットプリントの頂点情報)
measuredHeight geometry
0 12.0 MULTIPOLYGON Z (((139.73105 35.70079 29.13326,...
1 8.1 MULTIPOLYGON Z (((139.72776 35.70382 17.10876,...
2 17.2 MULTIPOLYGON Z (((139.73242 35.70671 4.49996, ...
3 25.1 MULTIPOLYGON Z (((139.73661 35.70784 4.78996, ...
4 6.8 MULTIPOLYGON Z (((139.72815 35.70801 5.40966, ...
5 10.1 MULTIPOLYGON Z (((139.72832 35.70525 9.67996, ...
6 8.5 MULTIPOLYGON Z (((139.73103 35.70249 26.17996,...
7 8.0 MULTIPOLYGON Z (((139.73092 35.70164 27.67116,...
8 7.4 MULTIPOLYGON Z (((139.73038 35.70137 28.38996,...
9 9.7 MULTIPOLYGON Z (((139.73504 35.70701 5.01006, ...この1列1列が、個々の建物のフットプリントと高さ情報を表しています。
ただのPolygonではなくMultiPolygonとなっているのは、建物が複数のメッシュモデルで構成されていることを意味しています。
2. フットプリントの頂点座標を平面直角座標系に変換
PLATEAUのデータに入っているフットプリントの頂点情報は地理座標系(緯度経度)で表されています。一方で高さ情報はメートル単位で表されており、メッシュを作成する際には同じ単位に合わせる必要があります。
今回はpyprojライブラリを用いてフットプリントの頂点座標を平面直角座標系(メートル単位)に変換することで、単位の統一を行います。
(例) 緯度, 経度の値を平面直角座標系の座標x, yに変換
import pyproj
EPSG4612 = pyproj.Proj("+init=EPSG:4612") # 地理座標系
EPSG2451 = pyproj.Proj("+init=EPSG:2451") # 平面直角座標系
# 地理座標系から平面直角座標系への変換
x, y = pyproj.transform(EPSG4612, EPSG2451, <Latitude>, <Longitude>)3. 建物メッシュを構成する天井面と床面の頂点verticesを生成
Open3Dのメッシュモデルは三角形の面(3つの頂点と各頂点を結ぶエッジ)の集合として表現され、メッシュモデルを作成する際には
・vertices(メッシュ頂点のリスト)
・triangles(「三角形を形成する3頂点の番号リスト」のリスト)
を設定する必要があります。
頂点のリストは、床面を構成する頂点のリストと、天井面を構成する頂点のリストを結合して作成しました。
4. 建物メッシュを構成する三角形trianglesを作成
建物のメッシュモデルを作成するため、3で作成した頂点のリストについて、どの頂点で三角形を形成するかを示すリストtrianglesを作成します。
ここがこの処理のキモです！シンプルな処理によって、建物を三角形のみで表現することが可能です。
3, 4の処理をまとめるとこんなイメージ
trianglesの作成では天井面、床面、側面1, 側面2という風に頂点番号を指定する処理を分け、それらを統合することで建物全体を表現できる三角形triangleのリストtrianglesを作っています。
# 側面のtrianglesを作成
triangles = []
for i in range(num_base_vertices - 1):
# 側面1 (図の緑の三角形)
triangles.append([i, (i + 1) % num_base_vertices, i + num_base_vertices])
# 側面2 (図の赤の三角形)
triangles.append(
[
(i + 1) % num_base_vertices,
(i + 1) % num_base_vertices + num_base_vertices,
i + num_base_vertices,
]
)
# 底面と天井面のtrianglesを作成
for i in range(2, num_base_vertices):
# 底面
triangles.append([0, i - 1, i])
# 天井面 (図の青の三角形)
triangles.append(
[num_base_vertices, num_base_vertices + i, num_base_vertices + i - 1]
)5. 建物メッシュモデルの作成
4, 5で作成したvertices, trianglesの情報から、メッシュモデルを作成します。
mesh = o3d.geometry.TriangleMesh()
mesh.vertices = o3d.utility.Vector3dVector(vertices)
mesh.triangles = o3d.utility.Vector3iVector(triangles)
6. 建物メッシュの集合(都市メッシュ)の作成
複数の建物メッシュを統合します。これらは最終的に一つのシーンで可視化されます。
7. 建物メッシュの色付け&可視化
Open3Dの paint_uniform_color メソッドを利用し、高さ閾値を超える建物とそうでない建物に対して異なる色を適用してみます。
今回はmeasuredHeightが30[m]以上の建物をマゼンタで、その他の建物を灰色で描画してみました。
可視化にはOpen3Dのvisualizationメソッドを使いました。
建物のモデルがきちんと作成されていることがわかります。
まとめ
今回はPLATEAUのCityGMLデータを読み取りOpen3Dのメッシュモデルに変換して描画する処理について解説しました。
「処理の流れ」では概要しか説明できていないので、詳しく知りたい方はぜひ実装を見てください。
実装
from typing import List, Tuple
import geopandas as gpd
import open3d as o3d
import pyproj
import tqdm
from shapely.geometry import MultiPolygon, Polygon
class CreateCityMesh:
def __init__(self, building_gdf: gpd.GeoDataFrame):
self.footprint: Tuple[Polygon, MultiPolygon] = building_gdf["geometry"]
self.height: float = building_gdf["measuredHeight"]
self.EPSG4612 = pyproj.Proj("+init=EPSG:4612")
self.EPSG2451 = pyproj.Proj("+init=EPSG:2451")
def _create_mesh_triangles(self, num_base_vertices: int) -> List[List[int]]:
"""
建物の3Dメッシュのtrianglesを作成する。
Args:
num_base_vertices(int): 建物の底面の頂点数
Returns:
List[List[int]]: 建物の3Dメッシュのtriangles
"""
# 側面のtrianglesを作成
triangles = []
for i in range(num_base_vertices - 1):
triangles.append([i, (i + 1) % num_base_vertices, i + num_base_vertices])
triangles.append(
[
(i + 1) % num_base_vertices,
(i + 1) % num_base_vertices + num_base_vertices,
i + num_base_vertices,
]
)
# 底面と天井面のtrianglesを作成
for i in range(2, num_base_vertices):
# 底面
triangles.append([0, i - 1, i])
# 天井面
triangles.append(
[num_base_vertices, num_base_vertices + i, num_base_vertices + i - 1]
)
return triangles
def _create_building_mesh(
self, bldg_footprint: Polygon
) -> Tuple[o3d.geometry.TriangleMesh, float]:
"""
建物のフットプリントの情報から3Dメッシュを作成する。
Args:
bldg_footprint(Polygon): 建物のフットプリント
Returns:
o3d.geometry.TriangleMesh: 建物の3Dメッシュ
"""
# フットプリントの外周座標・高さ情報を取得
exterior_info = np.array(bldg_footprint.exterior.coords)
exterior = exterior_info[:, :2]
# 緯度経度から平面直角座標系に変換
for i in range(len(exterior)):
x, y = pyproj.transform(
self.EPSG4612, self.EPSG2451, exterior[i, 0], exterior[i, 1]
)
exterior[i] = [x, y]
# 建物の底面と天井面の頂点を作成
bottom_vertices = np.hstack([exterior[:-1], np.zeros((len(exterior) - 1, 1))])
top_vertices = np.hstack(
[exterior[:-1], np.full((len(exterior) - 1, 1), self.height)]
)
# verticesのリストを作成
vertices = np.vstack([bottom_vertices, top_vertices])
num_base_vertices = len(bottom_vertices)
# trianglesを作成
triangles = self._create_mesh_triangles(num_base_vertices)
# open3dのTriangleMeshオブジェクトを作成
mesh = o3d.geometry.TriangleMesh()
mesh.vertices = o3d.utility.Vector3dVector(vertices)
mesh.triangles = o3d.utility.Vector3iVector(triangles)
mesh.compute_vertex_normals()
mesh.compute_triangle_normals()
return mesh
def _paint_mesh_by_height(
self,
mesh: o3d.geometry.TriangleMesh,
height_threshold: float = 30,
) -> o3d.geometry.TriangleMesh:
"""
heightがheight_thresholdより大きい場合はマゼンタ、
height_threshold以下の場合は灰色でメッシュを塗りつぶす。
Args:
mesh(o3d.geometry.TriangleMesh): 3Dメッシュ
height_threshold(float): 3Dメッシュを塗りつぶす高さの閾値
Returns:
o3d.geometry.TriangleMesh: 塗りつぶされた3Dメッシュ
"""
if self.height > height_threshold:
mesh.paint_uniform_color([1, 0, 1])
else:
mesh.paint_uniform_color([0.5, 0.5, 0.5])
return mesh
def create_city_mesh(
self,
meshes: List[o3d.geometry.TriangleMesh],
) -> List[o3d.geometry.TriangleMesh]:
"""
建物メッシュを作成し、meshes(建物の3Dメッシュの集合)に追加して返す。
Args:
meshes (List[o3d.geometry.TriangleMesh]): 建物の3Dメッシュの集合
Returns:
List[o3d.geometry.TriangleMesh]: 建物の3Dメッシュの集合
"""
if isinstance(self.footprint, Polygon):
mesh = self._create_building_mesh(self.footprint)
mesh = self._paint_mesh_by_height(mesh, 30)
meshes.append(mesh)
elif isinstance(self.footprint, MultiPolygon):
polygons = [polygon for polygon in self.footprint.geoms]
mesh_bldg = o3d.geometry.TriangleMesh()
for i, polygon in enumerate(polygons, start=1):
mesh = self._create_building_mesh(polygon)
mesh_bldg += mesh
mesh_bldg = self._paint_mesh_by_height(mesh_bldg, 30)
meshes.append(mesh_bldg)
return meshes
if __name__ == "__main__":
gdf = gpd.read_file("<gmlファイルのパス>")
num_buildings = len(gdf.index)
print(f"Number of buildings: {num_buildings}")
# measuredHeightのデータが無い場合は終了
if "measuredHeight" not in gdf.columns:
print("measuredHeight is not found")
exit(1)
meshes: List[o3d.geometry.TriangleMesh] = []
for i in tqdm.tqdm(range(num_buildings)):
building_gdf = gdf.iloc[i]
meshes = CreateCityMesh(building_gdf).create_city_mesh(meshes)
o3d.visualization.draw_geometries(meshes, window_name="Building Mesh")この記事が気に入ったらサポートをしてみませんか？