matplotlibで富士山近辺の立体地形図を描いてみた(駿河湾、相模湾の海底地形図付き)
前回まで書いてきた立体地図の作成結果を元に、海底地形図もついた富士山周辺の立体地形図をmatplotlib作成してみた。
これまでの、記事を読んでもらえば作るのは簡単で、 前回までの知識を元に、今回はMySQLの以下の標高格納用TABLEに国土地理と日本海洋データセンターの標高・水深データを、とりあえず地図タイルズームレベル7の解像度で、x軸方向のピクセル座標、y軸方向のピクセル座標ごとに一旦突っ込み、
あとはそこから指定の緯度経度の標高データを読み込み、メッシュしてsurface、contourするだけ。pythonなので超お気楽モード。
今回の範囲の地図を作成するのにズームレベル7解像度の246行201列分の国土地理と日本海洋データセンターの標高・深度データを合体して使用しています。
富士山は当然のことながら、駿河湾、相模湾あたりの海底地形や、伊豆半島、三浦半島、御前崎を始め、塩見岳、大無間山、甲武信ケ岳、蓼科山、御嶽山、槍ケ岳、丹沢山、箱根山、天城山の山々や、伊豆大島、利島、新島も描かれてますね。
なお、標高・水深を以下の感じで -3000m, -2000m, -1000m, -500m, 0m, 10m, 20m, 70m, 100m, 200m, 400m, 800m, 1200m, 1600m, 2000m, 3000m, 4000m ごとに色分けするカラーマップの作成にListedColormapとBoundaryNormを使っています。
標高 0m, 10m, 20m, 70m, 100m, 200mの区間は色分けが細かいので、拡大するとこんな感じ。私、配色センスがないので色分けはイマイチですが、色の指定はRGBの16進でも指定できますので、自分の好みに応じて自由に色分けを変えられます。
cmap = mpl.colors.ListedColormap(['blue', 'royalblue', 'cornflowerblue', 'dodgerblue', 'limegreen', 'g', 'darkkhaki', 'papayawhip', 'moccasin', 'wheat', 'lightsalmon', 'tomato', 'orangered', 'red', 'brown', 'snow']) cmap.set_over('white') cmap.set_under('darkblue') bounds = [-3000, -2000, -1000, -500, 0, 10, 20, 70, 100, 200, 400, 800, 1200, 1600, 2000, 3000, 4000] norm = mpl.colors.BoundaryNorm(bounds, cmap.N)
また、カラーマップの区間指定をちょっといじると、標高100mから200mの地点だけ色付けしたいなんてのも簡単にできちゃいます。
国土地理のデータについては著作権の問題もあるので、著作権フリーの地形データ、スペースシャトルの1秒メッシュ(約30m)地形データSRTM-1が使えるか次は試してみよう。
それと、標高単位で色分けしてるので、湖とか描かれてませんけどね。それは、要課題ですね。地形データをどっからかもってきてポリゴンで描けばいいのかな。
標高格納用TABLE
CREATE TABLE IF NOT EXISTS `dem_07` ( `x` int(12) NOT NULL, `y` int(12) NOT NULL, `depth` float NOT NULL ) ENGINE=MyISAM DEFAULT CHARSET=utf8; ALTER TABLE `dem_07` ADD PRIMARY KEY (`x`,`y`), ADD KEY `idx_x` (`x`), ADD KEY `idx_y` (`y`);
ソースはこんな感じ
from math import pi from math import tanh from math import sin from math import asin from numpy import arctanh import pandas as pd import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import matplotlib.ticker as ticker import matplotlib.colors as colors from matplotlib.mlab import bivariate_normal import MySQLdb e_db_host = "127.0.0.1" e_db_port = 3306 e_db = "earthquake" e_db_user = "ユーザー名" e_db_passwd = "パスワード" #緯度・経度からピクセル座標を計算する def fromLatLngToPoint(lat,lon, z): L = 85.05112878 pointX = int( (2**(z+7)) * ((lon/180) + 1) ) pointY = int( (2**(z+7)/pi) * (-1 * arctanh(sin(lat * pi/180)) + arctanh(sin(L * pi/180))) ) return pointX, pointY #ピクセル座標から緯度・経度を計算する def fromPointToLatLng(pixelX,pixelLonY, z): L = 85.05112878 lon = 180 * ((pixelX / 2.0**(z+7) ) - 1) lat = 180/pi * (asin(tanh(-pi/2**(z+7)*pixelLonY + arctanh(sin(pi/180*L))))) return lat, lon # ズームレベル z = 7 # 作図する地図の中心点の緯度経度 # 今回は富士山を中心に作図 look_Lat = 35.3625 look_Lon = 138.7306 # 中心点から緯度経度が1.1度離れた左下右上の座標を計算 gap = 1.1 minLat = look_Lat - gap minLon = look_Lon - gap maxLat = look_Lat + gap maxLon = look_Lon + gap point1X, point1Y = fromLatLngToPoint(minLat, minLon, z) point2X, point2Y = fromLatLngToPoint(maxLat, maxLon, z) # 標高メッシュデータを読み込む conn = MySQLdb.connect(db=e_db, host=e_db_host, port=e_db_port, user=e_db_user, passwd=e_db_passwd) cursor = conn.cursor() sql = "select * from dem_07 where x>=%d and x<=%d and y>=%d and y<=%d" % (point1X, point2X, point2Y, point1Y) cnt = cursor.execute( sql ) rec = cursor.fetchall() cursor.close() conn.close() # X座標、Y座標、標高地をそれぞれ配列にセーブ points = [] for point in rec: points.append(point) Xs = [ int(e[0]) for e in points ] Ys = [ int(e[1]) for e in points ] Elvs = [ float(e[2]) for e in points ] # ピクセル座標の最大、最小値などを求めておく xmin = min(Xs) xmax = max(Xs) xdiff = xmax-xmin ymin = min(Ys) ymax = max(Ys) ydiff = ymax-ymin # 標高データを格子データに格納する Z = np.zeros( (ydiff+1, xdiff+1) ) print Z.shape for (x, y, elv) in zip(Xs, Ys, Elvs): Z[y-ymin, x-xmin] = elv # ピクセル座標を緯度経度に変換してX,Y軸のグリッドを生成
# なを、標高メッシュの投影法が高緯度ほど間隔が伸びるメルカトル図法のためこれで生成される緯度は正確じゃないです。ここは要改造ポイントだな。 x = np.linspace(minLon,maxLon,xdiff+1) y = np.linspace(maxLat,minLat,ydiff+1) X, Y = np.meshgrid(x, y) # 立体地図を作成する fig = plt.figure(figsize=(24, 16), dpi=72, facecolor='w', edgecolor='k') plt.subplots_adjust(left=0, right=1, bottom=0, top=1) ax = fig.add_subplot(111, projection='3d') ax.view_init(elev=77, azim=-80) # 横を1/0.5=2倍長く設定 ax.set_aspect(0.5, adjustable='box') # 水深・標高用カラーマップ作成 cmap = mpl.colors.ListedColormap(['blue', 'royalblue', 'cornflowerblue', 'dodgerblue', 'limegreen', 'g', 'darkkhaki', 'papayawhip', 'moccasin', 'wheat', 'lightsalmon', 'tomato', 'orangered', 'red', 'brown', 'snow']) cmap.set_over('white') cmap.set_under('darkblue') bounds = [-3000, -2000, -1000, -500, 0, 10, 20, 70, 100, 200, 400, 800, 1200, 1600, 2000, 3000, 4000] norm = mpl.colors.BoundaryNorm(bounds, cmap.N) # 標高に応じて塗りつぶす surf = ax.plot_surface(X, Y, Z, norm=norm, rstride=1, cstride=1, linewidth=0, alpha=.7, cmap=cmap) # 標高100間隔で等高線を描く elevation = range(-10000,10000,100) cont = ax.contour(X, Y, Z, levels=elevation, edgecolor='none', cmap='binary', alpha=1.0) # 水深・標高バーをつける cb = fig.colorbar(surf, shrink=0.4, aspect=20, extend='both', ticks=bounds, spacing='proportional') cb.set_label('elevation') #plt.savefig('map.jpg', dpi=72) plt.show()