三重会合点付近の立体(3D)海底地形図をmatplotlibで作成してみた

前回作成した、日本海洋データセンター（JODC）の500ｍメッシュ水深データを元にした海底地形図を立体(3D)地形図にしてみた。

memomemokun.hateblo.jp

今回も、房総半島南東沖のフィリピン海プレートと太平洋プレートがぶつかり双方が北アメリカプレートの下に潜り込む、三重会合点付近。

相模トラフや日本海溝、伊豆・小笠原海溝や房総沖の日本海溝内に陥没しつつある第１鹿島海山を始め、第２、 第３、第４鹿島海山 、拓洋第２海山、拓洋第３海山あたりが描画されてますね〜〜。

 

元データは日本海洋データセンター（JODC）の500ｍメッシュ水深データ

http://jdoss1.jodc.go.jp/vpage/depth500_file_j.html

の以下ピンクで選択した16メッシュ、点数60万点ほどの深度データーから

ソースは以下

from math import pi
from math import tanh
from math import sin
from math import asin
from numpy import arctanh

from io import StringIO
import string

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.ticker as ticker


#緯度・経度からピクセル座標を計算する
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

#地形データを読み込み、ピクセル座標を付加し変数 maptxt に格納
f = open('jodc-depth500mesh-20181115192653.txt', 'r')
jodc = pd.read_csv(f, delim_whitespace=True, header=None, 
                      names=('type', 'lat', 'lon', 'depth'), 
                      dtype={'col1':'S1', 'col2':'f8', 'col3':'f8', 'col4':'i1'})

lats   = jodc['lat'].values
lons   = jodc['lon'].values
depths = jodc['depth'].values
depths = depths * -1.0

f = StringIO()
for (lat, lon, depth) in zip(lats, lons, depths):
    pointX, pointY = fromLatLngToPoint(lat, lon, z)
    txt = u'{0:.0f},{1:.0f},{2:.12f},{3:.12f},{4:.0f}'.format(pointX, pointY,lat,lon,depth)
    f.write(txt + "\n")
    
maptxt = f.getvalue()
f.close()


#変数 maptxtから地形データを読み込む
jodc = pd.read_csv(StringIO(maptxt), header=None, 
                         names=('x', 'y', 'lat', 'lon', 'depth'), 
                         dtype={'col1':'i1', 'col2':'i1', 'col3':'f8', 'col3':'f8', 'col5':'i1'})


#ピクセル座標の最大、最小値などを求めておく
xmin = jodc['x'].min()
xmax = jodc['x'].max()
xdiff = jodc['x'].max()-jodc['x'].min()
ymin = jodc['y'].min()
ymax = jodc['y'].max()
ydiff = jodc['y'].max()-jodc['y'].min()
depth = jodc['depth'].min()

#深度データの点群データを格子データに変換する
Density = np.zeros( (ydiff+1, xdiff+1) ) #格子内の点の数
Intensity = np.zeros( (ydiff+1, xdiff+1) ) #格子ごと深度の総和
for index, row in jodc.iterrows():
    #print row.x-xmin, row.y-ymin, row.depth
    x = int(row.x-xmin)
    y = int(row.y-ymin)
    Density[y,x] = Density[y,x] + 1
    Intensity[y,x] = Intensity[y,x] + row.depth

#格子ごと深度の総和 / 格子内の点の数　として各格子単位の深度データを求める
Z=Intensity/Density

#欠損値は取り敢えず0にする
Z[np.isnan(Z)] = 0.0

#ピクセル座標を緯度経度に変換してX,Y軸のグリッドを生成
minlat, minlon = fromPointToLatLng(xmin,ymin,z)
maxlat, maxlon = fromPointToLatLng(xmax,ymax,z)
x = np.linspace(minlon,maxlon,xdiff+1)
y = np.linspace(minlat,maxlat,ydiff+1)
X, Y = np.meshgrid(x, y)

#海底地図を作成する
fig = plt.figure(figsize=(14, 8), 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=76, azim=-84)

# 横を1/0.5=2倍長く設定
ax.set_aspect(0.5, adjustable='box') 


#深度に応じて塗りつぶす
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, linewidth=0, alpha=.8, cmap='jet')


#深度250m間隔で等高線を描く
elevation = range(-10000,10,250)
cont = ax.contour(X, Y, Z, levels=elevation, cmap='binary')

plt.show()