メッシュ水深データから海底地形図をmatplotlibで作成してみた
日本海洋データセンター（JODC）の
http://www.jodc.go.jp/jodcweb/index_j.html
500ｍメッシュ水深データで海底地形図を作成してみた。
J-EGG500のメッシュ水深データのデータの説明を読むと
『日本周辺の500mメッシュ海底地形データ （J-EGG500:JODC-Expert Grid data for Geography）は、海洋情報部をはじめとした 各種海洋調査機関によって得られた膨大な量の水深測量データを統合し、 多くの人が使用しやすいように等間隔で格子化した水深のデータセットです。』
『メッシュシステムは、各区域を２標準緯線によるランベルト正角円錐図法を使用してメートルスケールの平面座標に変換し、500 mの正方メッシュとして作成されています。経緯度に沿ったメッシュデータではありません。データは海域のみで、陸域のデータは含まれません。』
とあり、どうも、北極・南極を中心として緯線は同心円状に、経線は北極・南極から放射状に描くことで角度・形の歪みが少ないランベルト正角円錐図法上で500 m間隔の正方メッシュとして使用できるデータとして提供されているようなのですが、
Googleマップや国土地理院のタイル地図などで使われているタイル地図内のピクセル座標単位の経緯度に沿ったメッシュデータの方が、後々、国土地理院などの標高タイルデータと合体したりする際にも扱いやすいので、
J-EGG500の水深データを、タイル地図内のピクセル座標単位の格子状の水深データーに一度変換してからmatplotlibで海底地形図を描画してみました。
今回、作図に使ったのは、日本海洋データセンターの500mメッシュ水深データのダウンロードページ以下のように選択した、ピンク色の正方形の範囲の水深データ。
http://jdoss1.jodc.go.jp/vpage/depth500_file_j.html
フィリピン海プレートと太平洋プレートがぶつかり双方が北アメリカプレートの下に潜り込み、丁度、相模トラフ、日本海溝、伊豆・小笠原海溝の3つが出会う三重会合点と呼ばれる房総半島南東沖辺りの海底付近の海底地形平面図になります。
J-EGG500のデータフォーマットは以下の通り
・種別（0または1）、緯度（単位：度）、経度（単位：度)、水深（単位：ｍ）
・フォーマット[ I1、F10.5、F10.5、I6 ]
・種別 0： 計測水深または等深線から求めた水深、1： 補間処理により作成された水深 ・測地系は、世界測地系(WGS-84)を採用しています。
ソースはこんな感じ
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-20181115154754.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) X, Y = np.meshgrid(x, y) #海底地図を作成する fig = plt.figure(figsize=(14, 10), dpi=72, facecolor='w', edgecolor='k') plt.subplots_adjust(wspace=0.1, hspace=0.1) ax = fig.add_subplot(111) # ピクセル座標のY軸は北から南に値が大きくなるのでY軸は逆順にする ax.set_ylim([ymax, ymin]) #深度に応じて塗りつぶす #pcolor = ax.pcolor(X, Y, Z, cmap='Purples_r') pcolor = ax.pcolor(X, Y, Z, cmap='winter') #深度250m間隔で等高線を描く elevation = range(-10000,10,250) cont = ax.contour(X, Y, Z, levels=elevation, cmap='binary') cont.clabel(fmt='%1.1f', fontsize=14) ax.xaxis.get_major_formatter().set_useOffset(False) ax.yaxis.get_major_formatter().set_useOffset(False) ax.xaxis.set_major_locator(ticker.MaxNLocator(integer=True)) ax.yaxis.set_major_locator(ticker.MaxNLocator(integer=True)) #水深バーをつける cb = plt.colorbar(pcolor, shrink=0.5, aspect=10) #plt.savefig('map.jpg', dpi=72) plt.show()
上のソースに、以下コードを付け加え、10ピクセル単位に縦横に罫線を引いたのが以下で、各マス目に縦横10ピクセルで100個のピクセル座標が含まれている感じになります。
ax.xaxis.set_major_locator(ticker.MultipleLocator(100.0)) ax.xaxis.set_minor_locator(ticker.MultipleLocator(10.0)) ax.yaxis.set_major_locator(ticker.MultipleLocator(100.0)) ax.yaxis.set_minor_locator(ticker.MultipleLocator(10.0)) ax.grid(which='major', axis='x', linewidth=0.75, linestyle='-', color='0.75') ax.grid(which='minor', axis='x', linewidth=0.25, linestyle='-', color='0.75') ax.grid(which='major', axis='y', linewidth=0.75, linestyle='-', color='0.75') ax.grid(which='minor', axis='y', linewidth=0.25, linestyle='-', color='0.75')
さらに、ピクセル座標を緯度経度に直してX軸、Y軸を描くと、緯度経度も表示できます。
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 #ピクセル座標から緯度・経度を計算する 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-20181115160756.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, 10), dpi=72, facecolor='w', edgecolor='k') plt.subplots_adjust(wspace=0.1, hspace=0.1) ax = fig.add_subplot(111) #深度に応じて塗りつぶす #pcolor = ax.pcolor(X, Y, Z, cmap='Purples_r') pcolor = ax.pcolor(X, Y, Z, cmap='winter') #深度250m間隔で等高線を描く elevation = range(-10000,10,250) cont = ax.contour(X, Y, Z, levels=elevation, cmap='binary') cont.clabel(fmt='%1.1f', fontsize=14) ax.xaxis.get_major_formatter().set_useOffset(False) ax.yaxis.get_major_formatter().set_useOffset(False) ax.xaxis.set_major_locator(ticker.MaxNLocator(integer=True)) ax.yaxis.set_major_locator(ticker.MaxNLocator(integer=True)) ax.xaxis.set_major_locator(ticker.MultipleLocator(1.0)) ax.xaxis.set_minor_locator(ticker.MultipleLocator(0.1)) ax.yaxis.set_major_locator(ticker.MultipleLocator(1.0)) ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.1)) ax.grid(which='major', axis='x', linewidth=0.75, linestyle='-', color='0.75') ax.grid(which='minor', axis='x', linewidth=0.25, linestyle='-', color='0.75') ax.grid(which='major', axis='y', linewidth=0.75, linestyle='-', color='0.75') ax.grid(which='minor', axis='y', linewidth=0.25, linestyle='-', color='0.75') #水深バーをつける cb = plt.colorbar(pcolor, shrink=0.5, aspect=10) #plt.savefig('map.jpg', dpi=72) plt.show()
