伊豆・小笠原沖の3D海底地形図をplotlyでブログに表示してみました

先日、Python matplotlibで日本海洋データセンター（JODC）の500ｍメッシュ水深データで海底地形図を描いてみましたが、

memomemokun.hateblo.jp

立体図形をWeb上でも表示できるplotlyで海底地形図を描いてみました。

plot.ly

こんな感じ！

上で描いてみたのは以下のピンク色で囲まれた海域の水深データになります。

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

plotlyについて調べはじめたばかりで、使い方がよくわけってませんが、ほとんど前回コードを直すことなく実行出来ます。

ソースはこんな感じ。

mport plotly.plotly as py
from plotly.graph_objs import Surface


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 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


# ズームレベル
z = 7

#地形データを読み込み、ピクセル座標を付加し変数 maptxt に格納
f = open('jodc-depth500mesh-20181229134551.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軸のグリッドを生成
x = np.linspace(xmin,xmax,xdiff+1)
y = np.linspace(ymin,ymax,ydiff+1)

# Y軸は北の方が値が小さいので逆順にする
y = y[::-1]

X, Y = np.meshgrid(x, y)

# グリッド端に欠損値があるようなので周りの1ピクセルを省いちゃう
X = X[1:X.shape[0]-1,1:X.shape[1]-1]
Y = Y[1:Y.shape[0]-1,1:Y.shape[1]-1]
Z = Z[1:X.shape[0]-1,1:X.shape[1]-1]

surface = Surface(x=X, y=Y, z=Z, colorscale='Jet') 

data=[surface]
py.iplot(data,filename='Ocean bottom topography of Izu-Ogasawara')