Pythonの Cartopy 球面プロットでいい感じに遠近感を出す方法

はじめに

球座標データを可視化する際によく使われるのが、Orthographic投影（正射図法）です。これは球を離れた場所から見たようなプロットになるので、せっかくなので影をつけてもう少し格好良い図を作りたいと思うのが人の性でしょう（下図参照）。２次元プロットなのですが、まるで３次元かのような遠近感が出ています。

Paraview などの洗練された３次元用描画ソフトを使えば、自動でいい感じに影をつけてくれるのですが、ここでは Pythonの Cartopy ライブラリを利用した手軽な方法（1例）を紹介します。

環境

テスト環境は以下の通り。既に Cartopy はインストール済みとします。インストール方法はこちらの記事等参照。

• Mac OS Monterey: 12.3
• Python: 3.9.13
• Numpy: 1.23.1
• Matplotlib: 3.5.2
• Cartopy: 0.20.3

Cartopy を使った Orthographic 投影プロット

まずは、球面上定義されたデータを cartopy のorthographic projection を用いて、可視化してみます。

prep.py

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import matplotlib.patches as mpatches
import cartopy.crs as ccrs

#grid
N_theta = 100
N_phi = 200

theta = np.linspace(-0.5*np.pi,0.5*np.pi,N_theta)
phi = np.linspace(-np.pi,np.pi,N_phi)

x,y = np.meshgrid(phi,theta)
rad = 180.0/np.pi

#spherical data
data = np.zeros([N_theta,N_phi])
for j in range(N_theta):
    for k in range(N_phi):
        data[j,k] = (np.cos(theta[j])**6)*np.sin(6*phi[k])

#--- figure ---#
fig = plt.figure(figsize=(6,6))

#view point
lon0=0.0
lat0=30.0
proj_ortho = ccrs.Orthographic(central_longitude=lon0,central_latitude=lat0)

ax = fig.add_subplot(111,projection=proj_ortho)
ax.set_global()

#plot spherical data (contour map)
im = ax.pcolormesh(x*rad,y*rad,data,cmap='RdBu',transform=ccrs.PlateCarree(),shading='gouraud',vmin=-1,vmax=1)

#latitude & longitude lines
gl=ax.gridlines(crs=ccrs.PlateCarree(),draw_labels=False,linewidth=1,color='k',linestyle='--')
gl.xlocator = mticker.FixedLocator(np.linspace(-180,180,13))
gl.ylocator = mticker.FixedLocator(np.linspace(-90,90,7))

#output
plt.show()

これで以下のような画像が出力されるはずです。

ここでは、簡単のためデータは $l=m=6$ の球面調和関数で与えています。
lon0 と lat0 は視点調整パラメータです。lat0=90 とすると北極から見ることに対応します。

ちなみに、コンター図なのに contourf ではなく pcolormesh を使っているのは、contourf を使うとプロットに要する時間が果てしなく長くなるからです。一説によると、同じデータのプロットに要する時間を Basemap と比較したところ、cartopy の方が約２０倍遅くなっているそうです（僕もやってみたのですが、グリッド数が大きいと大抵の場合計算が終わらずエラーとなります）。
これでは流石に使い物にならないので、現状 pcolormesh で shading='gouraud' とすることで凌いでいます。これは cartopy が今現在抱えている非常に重大なバグでこちらの掲示板等でも議論されています。

不透明度を変化させた影エフェクト分布

影効果は、グレースケールの分布関数を不透明にして重ねてプロットすることで実現しようと思います。そこで、まずは影分布関数を作っていきます。

shadow.py

#compute a distance between 2 points on a sphere 
def caldistance(theta1,phi1,theta2,phi2):
    x1 = np.cos(theta1)*np.cos(phi1)
    y1 = np.cos(theta1)*np.sin(phi1)
    z1 = np.sin(theta1)
    x2 = np.cos(theta2)*np.cos(phi2)
    y2 = np.cos(theta2)*np.sin(phi2)
    z2 = np.sin(theta2)
    d = np.arccos(x1*x2+y1*y2+z1*z2)
    return(d)

#construct shade and its transparency
shad = np.zeros([N_theta,N_phi])
alpha = np.zeros([N_theta,N_phi])
for j in range(N_theta):
    for k in range(N_phi):
        d = caldistance(0,0,theta[j],phi[k])
        r = np.abs(2*d/np.pi)
        shad[j,k] = r
        alpha[j,k] = np.minimum(abs(r**4),1)

ここで shad は実際にプロットされる値関数で、alpha は不透明度です。どちらも円の中心で０、端で１を取るのですが、中心は影の効果をなるべく無くしたいので不透明度を中心からの距離の４次関数で与えて、端付近で局所的に不透明度が大きくなるようにしてあります。

影分布関数を球面上にプロット

既に、Orthographic の中心視点は (lat0, lon0) で与えているため、そのまま shad をプロットすると中心の座標がズレてプロットがおかしくなってしまいます。そこで、cartopy の RotatedPole という機能を用いて、以下のように投影座標を変換してあげます。

plot_shadow.py

#rotate the axis of shadow map
pole_lat = 90-lat0
pole_lon = 180
ccrs_rot = ccrs.RotatedPole(pole_latitude=pole_lat,pole_longitude=pole_lon)

#plot the shadow with varying transparency
im2 = ax.imshow(shad,alpha=alpha,cmap='binary',interpolation='bicubic',transform=ccrs_rot,vmin=0,vmax=1) 

ここでは、imshow を用いて影分布をプロットしました。「先ほどと同じように pcolormesh を使えばいいのでは？」と思われるかもしれませんが、実は pcolormesh を使うと円の端に行くに従って無数の黒線が発生してしまうという問題があります。これは、そもそも pcolormesh では仕様上長方形の格子同士がわずかに重なっているため、不透明度を下げたプロット法を用いると、必ず格子の間に模様が浮かび上がってきてしまうのです。
しかし一方で、 imshow にも問題があります。下図に示した通り、円の端付近でどうしてもデータ欠損が発生してしまうのです。解像度を大幅に上げれば誤差の範囲内で無視できるようになるのですが、流石にコスパが悪いです。ここでは、球の大枠の線幅を太めにすることでデータ欠損を隠すことにしました。

重ねてプロット

最後に、data と shad を重ねてプロットします。陽に zorder を与えることで shad が上にくるようにします。

overplot.py

#plot the data first
im = ax.pcolormesh(x*rad,y*rad,data,cmap='RdBu',transform=ccrs.PlateCarree(),shading='gouraud',vmin=-1,vmax=1,zorder=0)

#overplot the shadow
im2 = ax.imshow(shad,alpha=alpha,cmap='binary',interpolation='bicubic',transform=ccrs_rot,vmin=0,vmax=1,zorder=1)

#thicken the spherical boundary line
clip_circle = mpatches.Circle(xy=[0,0],radius=6370000,facecolor='none',edgecolor='k',linewidth=3)

最後の行の radius=6370000 は地球半径の近似値です。最終的に下のようになって完成です。

正直こういう用途においては Basemap の方が描画が楽だったので、開発が終わってしまったのは残念です。まだ上記の contourf 問題があったり、cartopy が basemap の完全上位互換とは呼べない状況なのかなと思います。今後、さらに開発・アップデートが進みより便利な機能が追加されることに期待しましょう。