プログラム日記φ(..)

おもにPython関係のプログラムメモ

matplotlibで富士山近辺の立体地形図を描いてみた(駿河湾、相模湾の海底地形図付き)

前回まで書いてきた立体地図の作成結果を元に、海底地形図もついた富士山周辺の立体地形図をmatplotlib作成してみた。

 

f:id:memomemokun:20181118175700j:plain

 

これまでの、記事を読んでもらえば作るのは簡単で、 前回までの知識を元に、今回は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を使っています。

 

f:id:memomemokun:20181119053542j:plain

 

標高 0m, 10m, 20m, 70m, 100m, 200mの区間は色分けが細かいので、拡大するとこんな感じ。私、配色センスがないので色分けはイマイチですが、色の指定はRGBの16進でも指定できますので、自分の好みに応じて自由に色分けを変えられます。

 

f:id:memomemokun:20181119060109j:plain

 

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の地点だけ色付けしたいなんてのも簡単にできちゃいます。

 

f:id:memomemokun:20181119054830j:plain

 

国土地理のデータについては著作権の問題もあるので、著作権フリーの地形データ、スペースシャトルの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()