matplotlibとMayaviでNASAの画像を『地球儀』にしてみた。

前回、地球の全体地図をPythonで描いてましたら

Pythonで地球儀を描いてみた

むむ！　

地球の全体図を球面にマッピングすれば地球儀が描けるじゃん！と思い至り

memomemokun.hateblo.jp

NASAが公開している、経緯度の間隔が方眼紙と同じ等間隔の直行平面上に単純にマッピングされた正距円筒図法（plate carrée）「blue marble」のイメージを使い今回はPythonで『地球儀』を描いてみることに。

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

使用したデーターは、縦横5400x2700の「blue marble」のイメージ

f:id:memomemokun:20181227163426j:plain

それで、半径rの球面の球面座標 (r,θ,φ) から直交直線座標 (x,y,z) への変換は

f:id:memomemokun:20181227172503p:plain

以下で求まりますので、

$x = r\sin \theta \cos \phi$
$y = r\sin \theta \sin \phi$
$z = r\cos \theta$

球面座標 θ、φの範囲0 ≤ θ ≤ π、0 ≤ φ < 2πからx、y、zの取りうる範囲を求め、求めたx、y、zから3Dで球面を描き、その表面をNASAの「blue marble」イメージで塗りつぶして、『地球儀』をmatplotlibを使い描いてみると、こんな感んじ！

ソースは、以下ですが

import matplotlib.pyplot as plt 
import numpy as np 
from mpl_toolkits.mplot3d import Axes3D 

def sphere_plot(image_file):

    img = plt.imread(image_file)
    
    # 画像を少し間引く
    img=img[::10,::10]
    
    
    # 球の球面座標　r:動径、θ:極角(天頂角)、φ：方位角(偏角)
    R = 1
    theta = np.linspace(0, np.pi, img.shape[0])
    phi = np.linspace(0, 2*np.pi, img.shape[1])

    # 球の球面座標 (r,θ,φ) を直交直線座標 (x,y,z) に変換
    x = R * np.outer(np.cos(phi), np.sin(theta))
    y = R * np.outer(np.sin(phi), np.sin(theta))
    z = R * np.outer(np.ones(np.size(phi)), np.cos(theta))
    

    # 球の表面をイメージファイルで塗りつぶして3Dで描画する
    fig = plt.figure(figsize=(10, 10)) 
    ax = fig.add_subplot(111, projection='3d') 
    plt.subplots_adjust(left=0, right=1, bottom=0, top=1)

    ax.plot_surface(x.T, y.T, z.T, rstride=1, cstride=1, facecolors = img/256.) 

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


if __name__ == "__main__":

    # 5400x2700
    image_file = 'world.topo.bathy.200412.3x5400x2700.jpg'
    sphere_plot(image_file)

matplotlibだとめちゃんこ時間がかかります。

ということで、3Dグラフィックスのための強力なテクスチャマッピングとレンダリングエンジン ( VTK )を搭載した対話的3次元プロットパッケージMayaviで『地球儀』を描いてみるとこんな感じ！

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

速い！　速い！　

matplotlibの1000倍は早いかもと思えるような超高速の一瞬で『地球儀』を表示してきました。

ズーマップも思いのまま

ソースは、以下

    
from mayavi import mlab
from tvtk.api import tvtk

def sphere_plot(image_file):

    # 図形ウィンドウをジェネレート
    fig = mlab.figure(size=(800, 800))

    # マップ画像をテクスチャにロードする
    img = tvtk.JPEGReader()
    img.file_name = image_file
    texture = tvtk.Texture(input_connection=img.output_port, interpolate=1)


    
    R = 1 # 球の半径
    Nrad = 1000 # テクスチャの分割数

    # テクスチャをマッピングする半径Rの球面を生成
    sphere = tvtk.TexturedSphereSource(radius=R, theta_resolution=Nrad,
                                       phi_resolution=Nrad)

    # 図形ウィンドウにテクスチャのマッパーとアクターを紐づける
    sphere_mapper = tvtk.PolyDataMapper(input_connection=sphere.output_port)
    sphere_actor = tvtk.Actor(mapper=sphere_mapper, texture=texture)
    fig.scene.add_actor(sphere_actor)
    
    mlab.show()



if __name__ == "__main__":

    # 5400x2700
    image_file = 'world.topo.bathy.200412.3x5400x2700.jpg'
    sphere_plot(image_file)