home about terms twitter

PythonとPlotlyで富士山の立体図をブラウザ上で表示する

date: 2020-02-20 | mod: 2020-02-20

Pythonを介したPlotlyの使用により、stlファイルとして得た富士山の立体データをxyz座標のみのデータに変換し、html出力してブラウザ上で表示させます。

  • Plotly …オープンソースのデータ分析・グラフ可視化ツール(公式サイト
  • stlファイル…法線ベクトルと三角形の座標を記録した三次元データファイル(Wikipedia

本ページ末に操作可能なサンプルのグラフを用意しています。


関連記事


index


コード全体

動作環境:

  • python 3.6.8
  • plotly 4.5.0
  • windows10 64bit
input output
dem.stl plotly_graph.html
import plotly.graph_objects as go
import numpy as np

data = open("dem.stl", "r")

with open('dem.txt', 'w') as f:
    lines = data.readlines()
    line_list = []
    for i in range(len(lines)):
        if 'outer loop' in lines[i - 1]:
            line = lines[i].replace('			 vertex ', '')
            line_list.append(line)

    for i in range(len(line_list)):
        if i % 50 == 0:
            f.write(line_list[i])

data.close()

pts = np.loadtxt('dem.txt')
x, y, z = pts.T

fig = go.Figure(data=[go.Mesh3d(x=x, y=z, z=y,
                intensity=y,
                colorscale='Inferno')]
                )

with open('plotly_graph.html', 'w') as f:
    f.write(fig.to_html(include_plotlyjs='cdn'))

後述しますが、stlファイルに記載されている座標を全て取得しておらず、座標データ数を150分の1に圧縮しています。ご理解いただいたうえでご使用ください。

重要な箇所を以下に示します。


データの取得

今回用いたstlファイルは、国土地理院が公開している地理院地図より取得しています(2020/02/20現在)。操作内容を簡単に示します。

上記URL先へ移動した後、まずは地図上の任意の座標へ移動します。
今回は富士山の山頂を選択しました。座標は[北緯: 35.363086°, 東経: 138.731060°]です。こちらへ移動すると、今回取得したものと同じ座標へ移動できます。

python-plotly-html-fujisan-1

任意の座標へ移動した後は、右上の【機能】->【3D】ボタンを移動していきます。大(2048 × 2048), 小(1024 × 1024)とありますが、今回はカスタムを選択しました。

python-plotly-html-fujisan-2

取得する領域を入力する画面が出てきます。今回は【400 ~ 400】と入力しました。

python-plotly-html-fujisan-3

そして【OK】を入力すると、以下のような画面が出てきます。

python-plotly-html-fujisan-4

立体的な富士山のイメージとともに、下部で各種ファイルが取得できるようになっています。
【STLファイル】、【VRMLファイル】、【WebGL用ファイル】とありますが、stlファイルを選択し、ダウンロードしました。

座標ファイルの変換

取得したstlファイルの内容は以下のようになっています。

solid 3d_data
	facet normal -0.918669216227036 0.3948009310723848 -0.013383082409230718
		outer loop
			 vertex 0 54.14 0
			 vertex 0 54.16 0.59
			 vertex 0.59 54.36 0
		endloop
	endfacet
	facet normal -0.9133536087789137 0.4034751157562689 -0.054708490272035296
		outer loop
			 vertex 0.59 54.36 0
			 vertex 0 54.16 0.59
			 vertex 0.59 54.44 0.59
		endloop
	endfacet
	facet normal -0.9232750122014152 0.38065627611484754 -0.051614410320656195
		outer loop
			 vertex 0.59 54.36 0
			 vertex 0.59 54.44 0.59
			 vertex 1.17 54.6 0
		endloop
	endfacet
    .
    .
    .

facet normalとして、三角形の法線ベクトルが記載されています。また、vertexとして、三角形の座標が示されています。
今回はこのvertexの座標のみを取得していきます。

Pythonで各行ごとにvertexが含まれるものを取得しても良いのですが、ファイルサイズが大きくなってしまうため、このプログラムではデータサイズを削減していくことを同時に行っていきます。具体的には「取得する座標の数を少なくしていく」ように設計していきます。

data = open("dem.stl", "r")

with open('dem.txt', 'w') as f:
    lines = data.readlines()
    line_list = []
    for i in range(len(lines)):
        if 'outer loop' in lines[i - 1]:
            line = lines[i].replace('			 vertex ', '')
            line_list.append(line)

    for i in range(len(line_list)):
        if i % 50 == 0:
            f.write(line_list[i])

data.close()

まず、ダウンロードしたdem.stlを変数dataとして読み込みます。
このdataについてreadlines()を使うことで、データの各行についてリスト(lines)として取得します。
得られたlinesについて、forループで処理を行います。

最初のループの中身としては、outer loopに一致したときに次の処理を行うようにしています。ここでのlines[i - 1]とは、目的とする座標データのひとつ前の行にouter loopが含まれていることを調べるために指定しています。三角形の三つの頂点の座標vertexのうち、outer loopが一つ前の行にある座標は一つしかないので、ここで3つ座標のうち1つだけを取得することで、取得する座標数を3分の1に削減しています。

そのような座標が見つかった後は、その座標の数値データのみが欲しいため、replace()を用いて vertex を消去しています。
最後に、得られた座標をリストline_listに追加しています。

最後のループでは、そのようにして得たline_listのうち、50の倍数の行にあるデータのみを取得します。この処理をすることで、取得する座標を50分の1に削減しています。先ほど3分の1にも削減したため、150分の1削減したことになります。

このようにして得られたデータは、data.txtとして出力します。このファイルの内容は次のようになります。

0 54.14 0
14.65 56.55 0
29.3 57.46 0
43.95 57.81 0
58.59 60.07 0
73.24 61.65 0
87.89 59.08 0
102.54 55.41 0
117.19 54.28 0
131.84 51.4 0
146.48 48.87 0
.
.
.

グラフパラメータの指定

pts = np.loadtxt('dem.txt')
x, y, z = pts.T

fig = go.Figure(data=[go.Mesh3d(x=x, y=z, z=y,
                intensity=y,
                colorscale='Inferno')]
                )

with open('plotly_graph.html', 'w') as f:
    f.write(fig.to_html(include_plotlyjs='cdn'))

得られたdem.txtをnumpyで読み込みます。読み込まれたデータは、1列目がx、2列目がy、3列目がzとして読み込まれます。

goとして読み込んだplotly.graph_objectsからFigure()を呼び出します。その中で各種パラメータを指定していきます。
go.Mesh3d()の引数を指定していきますが、今回得たstlファイルでは、富士山の高さの情報がy座標として与えられているようなので、y座標としてdem.txtのz座標、z座標としてdem.txtのy座標を指定しています。

また、intensity=y, colorscale='Inferno'として、yの数値によって色が変化するカラースケールを採用しています。

最後に、plotly_graph.htmlに出力させています。

出力

これにより得られる出力はhtmlですが、かなり長い(9万文字近い)ため詳細な内容は割愛させていただきます。
得られたhtmlの内容をブラウザ上で表示させると以下のようなデータとしてあらわされます。
<div>タグで囲まれた部分のみ貼り付け、<head>, <html>, <body>タグは省略しました。
あるいは出力したhtmlファイルを直接ブラウザで読み込ませることによっても見ることができます。

目的の立体データが出力されたことを確認できました。

関連記事