ちょっと寄り道。
(1)から始めたのは、仮想空間内にではなくて現実空間での「立体模型」を作ろうかな、てのが動機なので、仮想空間での作業は「等高線を作る」だけでいいわけね。だからこそ「専門の GIS ソフトウェアではない」 NumPy + matplotlib でも十分に用を足せる、ということでもあるわけ。それどころか「専門の GIS ソフトウェア」には欠けていたり使いにくかったりするようなことを、自由にやりたいわけよね。スケーリングの自由さとかそんなのね。
そうなんだけど、NumPy + matplotlib でも標高データの色んな可視化が出来るので、要するに「仮想空間でも結構遊べる」。ちょっとやってみる。
まずは一番基礎的なヤツ。
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5
6 #lat, lon = 35.224167, 139.025373 # 箱根駒ケ岳
7 lat, lon = 35.603693, 139.646342 # 等々力渓谷
8 elevs, lc, uc = converted_loader.load(lat, lon)
9
10
11 #
12 import matplotlib
13 import matplotlib.pyplot as plt
14 import matplotlib.cm as cm
15
16 fig, ax = plt.subplots()
17
18 # ・単独の 2d array を可視化するだけなので X, Y は渡せない
19 # ・南西端が (0, 0) の array にしてあるので南北転地
20 ax.imshow(np.flipud(elevs), cmap=cm.rainbow)
21 # tick は、なので無意味なので、消しておく
22 ax.set_xticks([])
23 ax.set_yticks([])
24
25 plt.show()
converted_loader は(5)のと同じものね。こんな結果：
標高データだとそれそのもので水平面方向の形状がわかるのであんまり必要ないけれど、例えば気象データの可視化なんかの場合は、「まずは手っ取り早く」に imshow することもあるけれど、それよりは basemapを使った海岸線表示と組み合わせた可視化をするかな、ふつうは。basemap は matplotlib のサブプロジェクト扱いと思う。使用感は matplotlib にほぼ同じ。なんかの機会があったら紹介してみてもいいけど今回はやらない。
さて。この単純な imshow でも結構「きゃーステキ」と思うが、影(shade)を付けるともっと雰囲気が出る：
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5
6 #lat, lon = 35.224167, 139.025373 # 箱根駒ケ岳
7 lat, lon = 35.603693, 139.646342 # 等々力渓谷
8 elevs, lc, uc = converted_loader.load(lat, lon)
9
10 # NaN がいるとカラーマップの処理がおかしい。
11 # NaN はゼロではないのでイヤなんだけどやむなし。
12 elevs[np.isnan(elevs)] = 0.
13
14
15 #
16 import matplotlib
17 import matplotlib.pyplot as plt
18 import matplotlib.cm as cm
19 from matplotlib.colors import LightSource
20
21 fig, ax = plt.subplots()
22
23 # ・単独の 2d array を可視化するだけなので X, Y は渡せない
24 # ・南西端が (0, 0) の array にしてあるので南北転地
25 elevs = np.flipud(elevs)
26 ls = LightSource(azdeg=0, altdeg=65)
27 rgb = ls.shade(elevs, cm.rainbow)
28 ax.imshow(rgb)
29 # tick は、なので無意味なので、消しておく
30 ax.set_xticks([])
31 ax.set_yticks([])
32
33 plt.show()
NaN の問題には悩まされた。masked_array でもダメみたいで。ゼロではないんだよなぁ…。等々力渓谷の場合、多摩川水面部分が欠損データになってるの。当然全然ゼロじゃあない。これの結果はこんな：
なかなかステキね。
ワイヤーフレーム：
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5
6 #lat, lon = 35.224167, 139.025373 # 箱根駒ケ岳
7 lat, lon = 35.603693, 139.646342 # 等々力渓谷
8 elevs, lc, uc = converted_loader.load(lat, lon)
9
10 # 高解像度過ぎて画像としてキツいので、間引く。
11 elevs = elevs[::3, ::3]
12
13
14 #
15 import matplotlib
16 import matplotlib.pyplot as plt
17 import matplotlib.cm as cm
18 from matplotlib.ticker import FormatStrFormatter
19 from mpl_toolkits.mplot3d import axes3d
20
21 fig = plt.figure()
22 ax = fig.gca(projection='3d')
23
24 ax.xaxis.set_major_formatter(FormatStrFormatter('%.4f'))
25 ax.yaxis.set_major_formatter(FormatStrFormatter('%.4f'))
26
27 X = np.linspace(lc[1], uc[1], elevs.shape[1])
28 Y = np.linspace(lc[0], uc[0], elevs.shape[0])
29 XX, YY = np.meshgrid(X, Y)
30 ax.plot_wireframe(XX, YY, elevs)
31
32 plt.show()
plt.show() はインタラクティブな GUI になっているので、回転や拡大をしてこんな感じ：
さすがにただのワイヤーフレームだと「残念」感が否めない。わかりにくいわぁ。
3d なら plot_surface が一番「ぽい」かな：
1 # -*- coding: utf-8 -*-
2 import numpy as np
3 import converted_loader
4
5
6 #lat, lon = 35.224167, 139.025373 # 箱根駒ケ岳
7 lat, lon = 35.603693, 139.646342 # 等々力渓谷
8 elevs, lc, uc = converted_loader.load(lat, lon)
9
10 # NaN がいるとカラーマップの処理がおかしい。
11 # NaN はゼロではないのでイヤなんだけどやむなし。
12 elevs[np.isnan(elevs)] = 0.
13
14 # 高解像度過ぎて画像としてキツいので、間引く。
15 elevs = elevs[::3, ::3]
16
17
18 #
19 import matplotlib
20 import matplotlib.pyplot as plt
21 import matplotlib.cm as cm
22 from matplotlib.ticker import FormatStrFormatter
23 from mpl_toolkits.mplot3d import axes3d
24
25 fig = plt.figure()
26 ax = fig.gca(projection='3d')
27
28 ax.xaxis.set_major_formatter(FormatStrFormatter('%.4f'))
29 ax.yaxis.set_major_formatter(FormatStrFormatter('%.4f'))
30
31 X = np.linspace(lc[1], uc[1], elevs.shape[1])
32 Y = np.linspace(lc[0], uc[0], elevs.shape[0])
33 XX, YY = np.meshgrid(X, Y)
34 cmap = cm.rainbow
35 surf = ax.plot_surface(XX, YY, elevs, cmap=cmap,
36 rstride=1, cstride=1)
37 fig.colorbar(surf)
38
39 plt.show()
これの結果は、さっきと同じくインタラクティブな GUI で回転や拡大してこんな感じ：
まぁまぁでしょ。無論こんなのは「専門のソフト」にはかなわないし、これがしたいだけなら Google Earth とか 地理院地図3D使えばいいだけのことです、ハイ。
あたしがこれらを使わずにやろうとしてるのはね、普通なら「補助情報を一緒に出したい」でしょうが、避けたいのはまさにそこだから。純粋に標高データだけを使った地図が欲しいわけ。Google Earth でサーフィス画像を非表示にすることは出来ないし、地理院地図3Dで地勢(建築物表示など)非表示に出来ない。地形だけを楽しむにはちょいと邪魔なのね。
