Pythonで3D地形図化 - Briswell Tech Blog

弊社のオフィスのある芝大門は、江戸時代より増上寺の門前町として多くの商店もあり賑わっています。以前、ご紹介しましたが江戸の大名屋敷もすぐそばにありました。戻れるものならこの目で見てみたいものです。

このあたりは、江戸時代に埋め立てがされ、それ以前は海の下（日比谷入江）だったようです。
国土地理院標高タイルの数値データより3D地形図を作成し、その事実を探っていきたいと思います。

matplotlibで国土地理院標高タイルから3D地形図を描いてみる - プログラム日記φ(．．)

こちらの記事を参考にさせていただきました。ありがとうございました。

では、行って参る。

1. モジュールのインポート

よく使われるモジュールですね。

import pandas as pd #csvファイル読み取り用
import requests #HTTP通信用
import numpy as np #数値計算用
import matplotlib.pyplot as plt #グラフ作成用
from io import StringIO #文字列をファイルのように扱える

2. 標高データの読込

こちらの国土地理院標高タイルについて、芝大門周辺の４つのタイルを取得します。

url_1 = 'http://cyberjapandata.gsi.go.jp/xyz/dem/13/7276/3226.txt' #芝大門を含むタイルNo.1
url_2 = 'http://cyberjapandata.gsi.go.jp/xyz/dem/13/7275/3226.txt' #西側のタイルNo.2
url_3 = 'http://cyberjapandata.gsi.go.jp/xyz/dem/13/7276/3227.txt' #南側のタイルNo.3
url_4 = 'http://cyberjapandata.gsi.go.jp/xyz/dem/13/7275/3227.txt' #南西側のタイルNo.4

urlList = [url_1, url_2, url_3, url_4]
z_result = []

for url in urlList:
  response = requests.get(url)
  maptxt = str.replace(response.text, u'e', u'-0.0')
  Z = pd.read_csv(StringIO(maptxt), header=None)
  z_result.append(Z)

3. 標高データを結合

読み込んだ４つのタイルの標高データを結合します。

z1 = np.concatenate((z_result[1],z_result[0]), axis = 1) #No.1とNo.2を結合
z2 = np.concatenate((z_result[3],z_result[2]), axis = 1) #No.3とNo.4を結合 

z_mix = np.concatenate((z1,z2), axis = 0) #No.1&2とNo.3&4を結合

4. 3Dグラフの設定

3Dグラフのデザイン部分です。

X, Y = np.meshgrid(np.linspace(0,255,512), np.linspace(255,0,512))
fig = plt.figure(figsize=(10, 8), dpi=80, facecolor='w', edgecolor='k') #w:白(White), k:黒 (Black)
ax = fig.gca(projection='3d') #3Dグラフ
ax.set_aspect('auto') #Axes3D currently only supports the aspect argument 'auto'.
ax.view_init(84, -67.5) #視点の設定

5. 3Dグラフを描画

標高値に応じて曲面を描きます。また、弊社のオフィス位置あたりに★を出力します。

surf = ax.plot_surface(X, Y, z_mix, cstride=1, rstride=1, cmap='rainbow', antialiased=False)
ax.text(140, 195, 2, "★", 'z', color='white', size='xx-large')
cb = plt.colorbar(surf, shrink=0.5, aspect=10)