matplotlib+cartopyで描いた全球プロット画像をダジックアースを利用してデジタル地球儀化して遊ぶ。
はじめに
雲の分布などの全球プロットを正距円筒図法などで描くと、極域が引き延ばされたりして印象がずいぶんと変わってしまう。一方で正射図法などで描くと全球を見るためには視点を変えた複数の画像を用意する必要があり、面倒である。何か良い方法がないか考えた所、正距円筒図法を正射図法に変換し、地球儀のようにインタラクティブに操作することができるソフトウェア「Dagik Earth」の画像フォルダを差し替えてしまえば、自前の画像を地球儀のように表示・操作することができると考えた。
以下ではmatplotlib+cartopyでダジックアース用に調整した全球画像を出力し、ダジックアースの画像フォルダと差し替えて表示する方法をまとめていく。
実行環境
本記事で紹介するプログラムの実行環境は前半Jupyter Notebook,後半は後述する理由により、シェルあるいはAnaconda Promptを想定しています。
OSはwindows8.1とCnetOSで動作確認をしました。いずれもAnaconda3を導入しています。
cartopyとxarrayはanacondaの標準パッケージに含まれていないので、Documentを参照しながらインストールしておきます。
ダジックアースのダウンロード
このページからコンテンツのダウンロードページを開き、その他の項目にある、コンテンツ作成用テンプレートをダウンロードする。ダウンロードしたファイルを展開してDagik_sample/data/imagesを開くと、map,screenフォルダが存在する。mapには正距円筒図法で書かれた地球画像が、screenにはmapの画像を説明する画像が保存されています。
このimagesフォルダ以下を自前の画像で置き換えてしまえば自前のコンテンツを作成できるというのが本記事の目標です。
画像の作成
データの用意
コンターを描くための元データがほしいのでNOAAなどから適当にnetcdfファイルをダウンロードしてくる(本記事では自前のnetcdfファイルを用いている)。netcdfファイルはxarrayで読み取ることにします。
ファイルの読み込み部分はデータの形式によって書き換える必要があるかもしれない。(例えばfortranバイナリの場合はnp.fromfileが使える)
追記
NOAAのサイトがアメリカの政府閉鎖の影響でアクセスできなくなっているようです。
一応こちら(京都大学さんのミラーサイト)が利用できそうです。
NOAAの閉鎖が解除されました。(1/30追記)
#In[1]
#必要なモジュールのインポート
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr
import numpy as np
import pandas as pd
import cartopy.util as cutil
from PIL import Image
import subprocess
#import os windows環境で利用
#In[2]
#netcdfファイルの読み込み
ds = xr.open_dataset('test.nc')
ds
#Out[2]
#<xarray.Dataset>
#Dimensions: (lat: 64, lon: 128, time: 366)
#Coordinates:
# * lon (lon) float32 0.0 2.8125 5.625 8.4375 ... 351.5625 354.375,357.1875
# * lat (lat) float32 87.8638 85.09653 82.31291 ... -85.09653 -87.8638
# * time (time) datetime64[ns] 2016-01-01 2016-01-02 ... 2016-12-31
#Data variables:
# pm25 (time, lat, lon) float64 ...
#In[3]
#全てnumpy.ndarrayに変換
da=ds['pm25'].values
dtime=pd.to_datetime(ds['time'].values,'D') #ここだけpandasのDatetimeIndexを使用
lon=ds['lon'].values
lat=ds['lat'].values
#経度方向にサイクリックになるようにデータを追加
data,cyclic_lon=cutil.add_cyclic_point(da,coord=lon)
xx,yy=np.meshgrid(cyclic_lon,lat)
データをxarray.dataarrayではなくnumpy配列にするのは数百枚の画像を出力するので、少しでも描画を高速にするためです。
またcutil.add_cyclic_point()は経度360°のデータを経度0°のデータと同一にして元データに追加してやらないと、本初子午線上のデータがないことになり、コンター図に白線ができて残念なことになってしまうからである(下図参照)。
ちなみにadd_cyclic_point()はxarray.dataarrayに対してうまく働いてくれないみたいなので、その点でもnumpy配列にしてしまった方が都合が良いです。
map画像とscreen画像の作成
データを読み込んだところで次はmatplotlibで全球コンターを枠線をすべてoffに設定してからプロットする。
なおmatplotlibでのプロットはオブジェクト指向インターフェイスを基本に記述しています。そのあたりについてはこの記事に非常に詳しく解説されています。
#In[4]
def map_plot(timestep=0):
fig = plt.figure(figsize=(20.48,10.24),dpi=150,frameon=False) #figsize
ax = fig.add_subplot(1,1,1
,projection=ccrs.PlateCarree(central_longitude=180.0))
ax.set_global() #プロットする範囲を全球に設定
ax.coastlines() #海岸線を描く
ax.stock_img() #cartopy標準の地球画像を描く(なくてもいいかも)
cf = ax.contourf(xx,yy,data[timestep,:,:]
,levels=levels #コンター間隔を指定
,cmap='viridis' #カラーマップを指定
,extend='max'
,transform=ccrs.PlateCarree())
fig.savefig(outdir+'/temp/temp_'+str(timestep)+'.jpg'
,bbox_inches='tight'
,edgecolor='none'
,pad_inches=0.0)
plt.close()
return cf
def screen_plot(cf,timestep=0,):
fig = plt.figure(figsize=(5.12,10.24),dpi=100)
ax=fig.add_subplot(1,1,1,facecolor='k')
cax= fig.add_axes([0.2,0.3,0.1,0.5]) #カラーバー用のaxesを作成
ax.text(0.05,0.95,title+'\n'+str(dtime[timestep]),color='w',fontsize=20)
cbar=plt.colorbar(cf,cax=cax,extend='max')
cbar.set_ticks(cticks)
cbar.set_label(clabel,fontsize=16,color='w')
cbar.ax.tick_params(colors='w',labelsize=18)#目盛の色とフォントサイズ
# pic = plt.axes([0.1,-0.09,0.4,0.4]) #Dagik Earthロゴ用
# pic.imshow(img)
fig.savefig(outdir+'/images/screen_'+str(timestep)+'.jpg'
,facecolor='k')
plt.close()
return 0
map_plot()では正距円筒図法で全球コンター図を描き、ax.contourfの返り値cfをカラーバーを描くために返却しています。
screen_plot()ではこのcfを利用してカラーバーをプロットしています。
map画像の調整
map画像の外枠を調整しても、画像に出力された際に細い外枠線が残ってしまい、また画像サイズも変わってしまうので、pillowを使って、トリミングと画像サイズの調整を行う。
#In[5]
def map_trim(timestep=0):
factx1,facty1=0.002e0,0.0025e0
factx2,facty2=0.998e0,0.9975e0
factxy=1 #解像度を上げるとき用(なくてもいいかも)
img=Image.open(outdir+'/temp/temp_'+str(timestep)+'.jpg')
figx,figy=img.size #map画像のサイズを取得
im_trim=img.crop((figx*factx1,figy*facty1,
figx*factx2,figy*facty2))
new_im_trim=im_trim.resize((int(2048*factxy),int(1024*factxy)))
new_im_trim.save(outdir+'/images/map/map_'+str(timestep)+'.jpg')
return 0
#上記の関数を実行するためのメイン関数
def main(timestep):
cf=map_plot(timestep=timestep)
screen_plot(cf,timestep=timestep)
map_trim(timestep=timestep)
return 0
グローバル変数の設定と作業用ディレクトリの作成
各プロット間で共通の変数は関数の引数にするよりもグローバル変数として扱ったほうが楽なので、ここではその定義を行う。また出力先のディレクトリもここで作成するようにする。
#In[6]
levels=np.arange(0,50,0.5) #コンター間隔の設定
cticks=np.arange(0,50.1,10)#カラーバーのticks
clabel='pm25[ng/m3]' #カラーバーのlabel
outdir='~/work' #出力先ディレクトリ
subprocess.check_call(['mkdir','-p',outdir+'/{temp,images/{map,screen}}'])#ディレクトリの作成
title='my title' #screen画像内にタイトルを入れる
#screen画像にDagik Earth のロゴを入れる
#img = np.array(Image.open('dagik.jpg'))#テンプレートから切り出して用意する
Windows環境ではsubprocessからディレクトリを作成できなかったのでos.makedirsで代用します。
subprocess.check_call(['mkdir','-p',outdir+'/{temp,images/{map,screen}}'])の行を下記と置き換えてください。
os.makedirs(outdir+'/temp',exist_ok=True)
os.makedirs(outdir+'/images/map',exist_ok=True)
os.makedirs(outdir+'/images/screen',exist_ok=True)
exist_ok=Trueはディレクトリが存在しても例外を返さないという意味です。
プログラムを実行してみる
いきなりすべての画像を生成してしまうと、やり直すのに時間がかかってしまうのでまずは1枚だけ吐き出してみる。
#In[7]
main(0)
大丈夫そうならループで実行。
#In[8]
for i in range(366):
main(i)
詳しい原因はよくわからないが、jupyter notebook上で上記のプログラムを実行すると100枚あたりからPCのメモリをほとんど消費してしまう現象が起きた。ループを回すときは.py形式にしてシェルで実行する方がよいかも。
multiprocessingで効率的に描画する。
上記のプログラム普通に実行すると数百枚の画像を吐き出すのにそれなりに時間がかかってしまう。ここではpythonの並列化の機能を用いて、より効率的に画像を生成できるようにする。変更はそれほど多くありません。
まずプログラムの最初にインポートするモジュールにmultiprocessingを追加します。
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr
import numpy as np
import pandas as pd
import cartopy.util as cutil
from PIL import Image
import subprocess
import multiprocessing as mp #追加
#plt.switch_backend('agg') #環境によってはmatplotlibのバックエンドを指定する必要があるかもしれません
#<以下省略>
次にループでmain関数を実行していたところを下記のように書き換えます。
#<省略>
if __name__=='__main__':
pool=mp.Pool(4) #プロセスを4個生成
pool.map(main,range(366)) #366枚の画像を4プロセスで並列に生成する。
pool.close() #画像の生成が終了したらプロセスを終了させる
特にWindows環境ではif__name__=='__main__'がないと並列化できません。詳しくはこちらを参照。
使用できるプロセスの最大数はmp.cpu_count()で取得できます。
実行はシェルから行います。
$ python3 dagik_multi.py
> python dagik_multi.py
コード
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr
import numpy as np
import pandas as pd
import cartopy.util as cutil
from PIL import Image
import subprocess
import multiprocessing as mp
#plt.switch_backend('agg')
levels=np.arange(0,50,0.5) #コンター間隔の設定
cticks=np.arange(0,50.1,10)#カラーバーのticks
clabel='pm25[ng/m3]' #カラーバーのlabel
outdir='~/work'
subprocess.check_call(['mkdir','-p',outdir+'/{temp,images/{map,screen}}'])#ディレクトリの作成
title='my title'
#img = np.array(Image.open('dagik.jpg'))
#データの準備
ds = xr.open_dataset('test.nc')
da=ds['pm25'].values
dtime=pd.to_datetime(ds['time'].values,'D') #ここだけpandasのDatetimeIndexを使用
lon=ds['lon'].values
lat=ds['lat'].values
data,cyclic_lon=cutil.add_cyclic_point(da,coord=lon)
xx,yy=np.meshgrid(cyclic_lon,lat)
def map_plot(timestep=0):
fig = plt.figure(figsize=(20.48,10.24),dpi=150,frameon=False) #figsize
ax = fig.add_subplot(1,1,1
,projection=ccrs.PlateCarree(central_longitude=180.0))
ax.set_global() #プロットする範囲を全球に設定
ax.coastlines() #海岸線を描く
ax.stock_img() #cartopy標準の地球画像を描く(なくてもいいかも)
cf = ax.contourf(xx,yy,data[timestep,:,:]
,levels=levels
,cmap='viridis'
,extend='max'
,transform=ccrs.PlateCarree())
fig.savefig(outdir+'/temp/temp_'+str(timestep)+'.jpg'
,bbox_inches='tight'
,edgecolor='none'
,pad_inches=0.0)
plt.close()
return cf
def screen_plot(cf,timestep=0,):
fig = plt.figure(figsize=(5.12,10.24),dpi=100)
ax=fig.add_subplot(1,1,1,facecolor='k')
cax= fig.add_axes([0.2,0.3,0.1,0.5]) #カラーバー用のaxesを作成
ax.text(0.05,0.90,title+'\n'+str(dtime[timestep]),color='w',fontsize=20)
cbar=plt.colorbar(cf,cax=cax,extend='max')
cbar.set_ticks(cticks)
cbar.set_label(clabel,fontsize=16,color='w')
cbar.ax.tick_params(colors='w',labelsize=18)#目盛の色とフォントサイズ
# pic = plt.axes([0.1,-0.09,0.4,0.4]) #Dagik Earthロゴ用
# pic.imshow(img)
fig.savefig(outdir+'/images//screen/screen_'+str(timestep)+'.jpg'
,facecolor='k')
plt.close()
return 0
def map_trim(timestep=0):
factx1,facty1=0.002e0,0.0025e0
factx2,facty2=0.998e0,0.9975e0
factxy=1
img=Image.open(outdir+'/temp/temp_'+str(timestep)+'.jpg')
figx,figy=img.size
im_trim=img.crop((figx*factx1,figy*facty1,
figx*factx2,figy*facty2))
new_im_trim=im_trim.resize((int(2048*factxy),int(1024*factxy)))
new_im_trim.save(outdir+'/images/map/map_'+str(timestep)+'.jpg')
return 0
def main(timestep):
cf=map_plot(timestep=timestep)
screen_plot(cf,timestep=timestep)
map_trim(timestep=timestep)
return 0
if __name__=='__main__':
pool=mp.Pool(4)
pool.map(main,range(366))
pool.close()
ディレクトリ操作を含んでいるので、実行は慎重に行ってください。
描画例
画像の生成が終わったらimagesをダジックアースのimagesフォルダと置換します。置換が終了したら、ダジックアースフォルダのDagik_Earth_64bit.exeを実行します。map画像の枚数が500枚以上だと実行直後に「応答なし」と表示されることがありますが、画像の読み込みは行われているのでタスクマネージャーのメモリ使用率でも眺めながら気長に待ちます。900枚程度の画像がある場合はメモリが8GBほどのPCではメモリ不足でクラッシュする可能性があります。その場合はメモリの増設を検討しましょう。
以下実行例です。
ずっと眺めていられます。
まとめ
画像の差し替えを行うことで高度な専門知識がなくてもGoogleEarthっぽいことができました。ダジックアース様様です。
ミニプロジェクターで球面に投影したりすると結構楽しいです。
matplotlibはいろいろなことができて便利なのですが、描画速度が遅く、数百枚ともなると並列化しても10分～の時間がかかってしまいます。
また環境が変わるとmatplotlibがしばしばエラーを出すのも悩みどころ...
参考URL
2015年7月版ダジック・アース・マニュアル
Python公式ドキュメント
https://qiita.com/skotaro/items/08dc0b8c5704c94eafb9
https://qiita.com/yukiB/items/203a6248c2d466b80d38