Rhinoceros x Pythonで3Dデータの分析をしよう(実践編)

# はじめに

https://qiita.com/pensuke_san/items/092db0188d2901538389

こちらの前回の記事で、RhinocerosをPythonを使って操作する環境のセットアップができました。本記事では実際に3Dデータを使用したRhinocerosにおける3Dデータの分析例をご紹介したいと思います。ここではオープンな3Dの都市データが提供されているPlateauのデータを使用します。

# Plateauデータセットの準備

## Plateauとは

Plateauとは国土交通省日本の都市の3Dモデル化プロジェクトで、3D都市モデルを整備し、様々な活用の取り組みがなされています。誰でもデータを利用できるように整備された3D都市モデルがオープンなデータとして提供されています。

## データの準備

PlateauのWebサイトからデータをダウンロードします。正確にはG空間情報センターからProject PLATEAUのデータをダウンロードすることができます。

https://www.geospatial.jp/ckan/dataset/plateau

本記事では2020年東京23区のデータを使用します。

## データの仕様

### データ形式
　
PlateauのデータセットはCityGMLという形式で配布されています。CityGMLは3D都市モデルを扱う標準的なフォーマットであり、すべての情報が含まれています。ただし他の3Dモデリングソフトも同様ですが、Rhinocerosはそのままの形式では対応していないためデータ形式を変換する必要があります。データの詳細やデータ変換の方法については公式HPに詳細な記載がありますので、こちらを見ていただけると良いと思います。

https://www.mlit.go.jp/plateau/learning/tpc03-1/

簡素化されたデータではあるものの、ポータルサイトに.OBJ形式にすでに変換されているデータがアップロードされているのでそちらを使います。


### 各種データの説明
ダウンロードしたデータはフォルダごとに以下の種類に分かれています。

- bldg:建築
- brid:橋
- tran:道路
- dem:地形

またLOD（Level of Details）ごとに建物の詳細度合いが異なります。はじめから変換されているOBJファイルは一部の都心の地域のみLOD2がありテクスチャーなどリッチな情報を持っています。今回はbldgのlod1のみ使用します。
そしてそれぞれの種類のフォルダごとには以下の命名規則でデータが入っています。

```
# bldgの場合
./bldg/lod1/(メッシュコード)_bldg_6677_obj/(メッシュコード)_bldg_6677.obj
```

Plateauではすべての都市モデルが一つのモデルになっているのではなく、メッシュコードごとに分割されて保存されています。メッシュコードと実際の地域の関係はデータの中に`13100_indexmap_op.pdf`というファイルがあるのでそれを参照してください。
また6677とは東京都の一部が属するEPSGコードと呼ばれるもので平面直角座標系と呼ばれる、本来3次元の球体の地球を2次元で表現するときにある緯度経度の基準点を原点としてモデルを平面に変換しています。そのため実際の地図の緯度経度からモデル内の位置を割り出すにはこちらの基準点から相対的な位置座標を計算する必要があります。

それぞれのデータの単位はm単位で、一つのメッシュのモデルは1km x 1kmのサイズとなっています。モデルや読み込むソフトによっては縮尺スケールを気にする必要があるようなのですが、（例えばUnityで都市モデルの中を人や車が動くときには他のオブジェクトのスケールと対応させる必要がある）OBJファイルの場合はおよそメッシュのサイズが1kmに一致していたので、さほど気にしなくて問題なさそうです。

# 分析チュートリアル
それではPlateauの3DモデルをRhinocerosとPythonを使って分析してみます。
今回の分析としては指定のメッシュのエリアにおける建物のそれぞれの高さを取得しどのような高さの分布になっているか可視化するというのをゴールとしたいと思います。

## 1. 3Dデータを取り込み建物ごとに分解する

3Dモデルを取り込むには次の関数を使用します。
```python
import rhinoscriptsyntax as rs

def model_import(path):
    command = "-_Import {}".format(path)
    command += "_Enter "
    rs.Command(command, True)
```
`rs.Command()`はRhinocerosのGUIのコマンドをそのまま使用することができます。また文字列でOptionを指定することができます。
特徴的なエリアの例としてMESH_ID=53394611の東京駅のエリアを取り込みます。（以下の画像は取り込んだモデルをPerspectiveビューを選択し、RenderedモードでZoom extentsコマンドで視点調整したものになります。

![image-8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/178986/e96742b2-cf6d-7397-111e-829422c97f19.png)

ちなみにlod2のテクスチャー用モデルを取り込むと以下のようになります。

![image-9.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/178986/70b94cec-9cb5-de5c-a338-7aafc377d30d.png)

分布はオフィスビルなどの高い建物と一般的な商業施設や住居などの建物に分かれそうですね。早速それぞれの建物の高さを取得していきたいところですが、3Dモデルはエリアのメッシュで1つのオブジェクトとなっているので、建物それぞれに分解する必要があります。分解するためには次の処理を行います。

```python
import rhinoscriptsyntax as rs

objs = rs.AllObjects()
for obj in objs:
  exploded_objs = rs.ExplodeMeshes(obj, delete=True)
```

`rs.AllObjects()`はよく使う関数でエディタに取り込んでいるすべてのオブジェクトのGUID（uuid形式の一意なID）を取得することができます。それらを`rs.ExplodeMeshes(obj)`に渡すことで結合したMeshを建物一つ一つに分解できます。下の画像のように建物を分解し、単体で選択できるようになりました。

![image-10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/178986/553499d7-3d33-8b14-1b8a-e73b688f855d.png)

## 2. 分解した建物データのそれぞれの高さを取得する

次に分解した建物のそれぞれの高さを取得し、建物ごとに高さを記録したデータを作成してpandasなどで扱いやすいようにjsonファイルに保存します。一つ前のセクションで結合したメッシュを分解し建物一つ一つに扱えるようになりました。
早速建物の高さを測っていきますが、3Dオブジェクトは形状が複雑なので大まかな形状を取得するときにはバウンディングボックスを使います。

```python
import rhinoscriptsyntax as rs

def calc_bbox_coord(obj):
    points = rs.BoundingBox(obj)
    x_list = []
    y_list = []
    z_list = []
    for point in points:
        x_list.append(point.X)
        y_list.append(point.Y)
        z_list.append(point.Z)
    return {'x_min': min(x_list), 'x_max': max(x_list),
            'y_min': min(y_list), 'y_max': max(y_list),
            'z_min': min(z_list), 'z_max': max(z_list)}
```

`rs.BoundingBox`はオブジェクトのboundingboxの8点が返ってきます。今回は高さだけ分かれば良いので、各軸のmin・maxに直しています。この関数を使用して各建物のGUIDをキーとして、建物情報を格納したjsonファイルを作ります。

```python
import json
import rhinoscriptsyntax as rs

objs = rs.AllObjects()
data = {}
for obj in objs:
    bbox = calc_bbox_coord(obj)
    data[obj.ToString()] = bbox

filepath = "/path/to/result.json"
with open(filepath, "w") as f:
  json.dump(data, f, indent=4)
```

RhinocerosのオブジェクトはC#のObjectクラスを継承しているので、ToString()メソッドを使うことができます。GUIDごとにバウンディングボックスの情報を格納することができました。
```json
{
    "fd92cd12-dd84-4b66-989e-3f6fdeaf71c1": {
        "x_min": -5554.2038,
        "x_max": -5539.8295,
        "y_min": -35263.2091,
        "y_max": -35254.3327,
        "z_min": 3.915,
        "z_max": 13.918
    },
}
```

## 3. データを分析化して可視化する

ここまで来たら後はpandasやmatplotlibで各種前処理・可視化をするだけです。建物の高さの分布を可視化してみましょう。

```python
import pandas as pd
import matplotlib.pyplot as plt

df = pd.read_json("./result.json",orient='index')
df['height'] = df['z_max'] - df['z_min']

plt.hist(df['height'],bins=40,edgecolor='black')
plt.title('Height Distribution')
plt.xlabel('Height [m]')
plt.ylabel('The Number of Buildings')
plt.show()
```

![image-11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/178986/4d90dc02-81b9-5962-7aa0-f925bb43a5fc.png)

ほとんどが50m以下の建物ですが、やはり東京の中心地であり100m以上ビルや中には200m以上の超高層ビルも見られます。また違ったエリアでは異なる特徴が見られると思いますので、色々なエリアや建物の数や密度など異なる特徴を比較してみるのも面白いと思います。

## 4. まとめ

本記事ではPlateauのオープンなデータとRhinocerosを使って3Dデータの分析を行ってみました。Rhinocerosは比較的Pythonから扱いやすく、モデリングを自動化するだけでなく、大量の3Dオブジェクトの集計値を求めることにも活用することができます。この記事が読者の方の3Dデータを活用した開発や分析に興味を持つきっかけとなれば幸いです。