ETOPO1グローバル地形データセットで地球の全体地図をPythonで描いてみる

旧NGDC（米国地球物理データセンター）、現NCEI（国立環境情報センター）が提供している地形と海底地形を統合した無料のグローバル地形データセット。ETOPO1で地球全体の地図を描いてみました。

こんな感じ！

ETOPO1グローバル地形データセットで地球の全体地図を描いてみる

ETOPO1は、南極・北極の極域も含め全地球の海陸あわせた1分=約1.8kmの解像度の標高・水深のグリッドデータで、世界の海洋の体積を計算したり、地表のヒポソグラフ曲線を導出するために使用されているそうです。

で、提供されているデーターは南極とグリーンランドの氷床の上をなぞるものと「ETOPO1 Ice Surface」、氷床の下の地面をなぞるもの「ETOPO1 Bedrock」の２タイプがあるようで、それぞれ、メッシュのセルの中心に緯度と経度の線が通る「grid-registered」と、セルの端が緯度と経度の線に沿った「cell-registered」の２つの形式でデータが提供されています。

データの中身をのぞいてみると、メッシュの幅は地球全体で東西に21601セル、南北に10801セルの全体で約2億（233,312,401）セル、テレビやデジカメで言ったら地球全体で２０K相当のピクセル数のメッシュとして標高・水深データが格納されているようです。

ソースはこんな感じ

地球全体を一気に描くには233,312,401セルとデーターがデカイので、東西、南北ともメッシュを１０分の１に間引いて、２K、フルハイビジョン相当のピクセル数を使い地球全体を描いてみました。ETOPO1のことを調べはじめてから約２〜３時間。データの格納のされ方と、取り出し方がわかれば後はこっちのもの。データーを好きに加工して、やっつけのあんちょこソースですが描けるもんですね〜〜！　お気軽Pythonの本領発揮と言ったところでしょうか。

なお、以下ソースでは、X、Yの各ピクセルは小学校の算数でもよく使う方眼紙状のXY軸のメモリが等間隔の直行平面上に単純にマッピングしたNASAのBlue Marbleイメージと同じ正距円筒図法（plate carrée）となります。

NASA Visible Earth: December, Blue Marble Next Generation w/ Topography and Bathymetry

上の画像に前回描いた感じで、陰影起伏を重ねたらより立体的に表現できますし

memomemokun.hateblo.jp

Mayaviあたりを使い上の画像を球にテクスチャマッピングすれば、地球儀のように３次元でぐりぐり回せそうですね。

Mayavi: 3D scientific data visualization and plotting in Python — mayavi 4.6.2 documentation

なを、作図に使用したのは、南極とグリーンランドの氷床の下の地盤をなぞる「ETOPO1 Bedrock」、経緯線がメッシュのセルの中心（grid-registerednet）を通るnetCDF形式のものです。

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
    
    
from netCDF4 import Dataset
import numpy as np

import matplotlib.pyplot as plt
import matplotlib.colors as colors
from matplotlib import cm



etopo1 = Dataset('ETOPO1_Bed_g_gmt4.grd','r')

x=etopo1.variables['x'][:]
y=etopo1.variables['y'][:]
Z=etopo1.variables['z'][:]

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

# メッシュを間引く
X=X[::10,::10]
Y=Y[::10,::10]
Z=Z[::10,::10]


# 標高・水深地図を作成する
fig = plt.figure(figsize=(30, 20), dpi=72, facecolor='w', edgecolor='k')
plt.subplots_adjust(wspace=0.1, hspace=0.1)
ax = fig.add_subplot(111)


# 水深・標高用カラーマップ作成
cmap = 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 = [-5000, -4000, -2000, -500, 0, 100, 200, 500, 700, 1000, 1500, 2000, 3000, 4000, 5000, 6000, 7000]
norm = colors.BoundaryNorm(bounds, cmap.N)



# 標高・水深に応じて地形を塗りつぶす
pcolor = ax.pcolor(X, Y, Z, norm=norm, 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)
#cont = ax.contour(X, Y, Z, levels=elevation, cmap='seismic')
#cont = ax.contour(X, Y, Z, levels=elevation, cmap='coolwarm')
#cont.clabel(fmt='%1.1f', fontsize=14)

# 水深・標高バーをつける
cb = fig.colorbar(pcolor, shrink=0.6, aspect=20, extend='both', ticks=bounds, spacing='proportional')
cb.set_label('elevation')

plt.savefig('1.jpg', dpi=72)
plt.show()

上のソースを改良しまして

ETOPO1グローバル地形データセットで、地球全体を標高1m毎に別の色で塗り絵してみたのがこちら。

memomemokun.hateblo.jp