NOAA/GOES の見慣れたX-ray flux 3days plotをPythonで自力で作る
注: 使用するデータなど、かなり専門的な内容が含まれています。物理や用語の説明等は然るべき文献に譲ります。検索等から見に来てくださった方、一部モジュールの解説は4. コード各部の解説から行っていますので、そちらをご覧ください。
同業者の方々、車輪の再発明は承知していますが、もしもっと画期的な方法があればぜひご教授ください。ただ、IDLからPythonへの脱却は1つの大きな前提とさせていただきます。
1. 背景
1.1. 何をどうしてわざわざ自分で作るのか
太陽関連の研究をしていればどこかで見たことがあるはずの、以下の図を自分で作ります。作った図は5. 実際に描画するにありますので、そこだけでも是非見ていってください。
最近~~(最近見てなかったのでいつからか分からない)~~、GOES Xray fluxのページがインタラクティブなグラフに変わりました。これが問題だという訳ではないのですが、
- 確かに値は見やすいが、何時何分にフラックスがどうだったかという情報をリアルタイムで必要とする機会は正直ない
- 長い期間このフォーマットのX線タイムプロファイルが使われてきたので、このフォーマットの画像を見ると安心する人が多い
- 過去の3days plotは上記のページあるいはSolar Monitorから取得できるが、スライド等に直接貼り付けると線が細かったりで怒られる
- GOES-16の稼働に伴ってGOES-15,14の運用が終了する(らしい)ので、今後同じフォーマットの図が更新されるのか分からない
など、色々思うところがあるので、それならmatplotlibの勉強がてら、いっそのこと自分で作れるようになろう、というのがこの投稿の趣旨です。
1.2. 前提知識
この記事では、こちらの記事中の
にて登場するaxesという言葉を多用します。以下の内容及びコードを少しでも理解することを目的とするならば、FigureとAxesを簡単に把握した上で読み進めることをおすすめします。
また、序文で記載の通り、専門用語等の解説はしていません。あらかじめご了承ください。モジュール等の解説は可能な限り専門用語を排して行うつもりです。
2. プログラムの作成
2.1. モジュールのimport
今回、見た目にこだわるために色々と初めて使うモジュールをimportしました。
import numpy as np
import datetime,calendar,requests
import pandas as pd
from io import StringIO
import matplotlib
%matplotlib inline ## jupyter notebook上での描画用
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import matplotlib.transforms as transforms
from matplotlib.dates import DateFormatter
from matplotlib.dates import date2num
from matplotlib.offsetbox import AnchoredOffsetbox, TextArea, HPacker, VPacker
2.2. データの置き場所
基本的にGOES 15までのデータは全て公開されています。例えば、2017年9月、GOES-15のX線1分平均データはcsv形式でここに置いてあります。今回は、ローカルにファイルをダウンロードすることなくグラフを作成できるよう進めていきます。
また現時点では、グラフを作成したい時期のGOES primaryとGOES secondaryの番号が分かっている前提で進めます。確かIDLのコマンドにはprimaryの情報は入っていたと思うので、可能そうであればいつか追記します。
3days plotを作成するため、まずは元データを参照します。GOESの1分平均X線データは1ヶ月ごとにまとまってcsv形式で公開されているので、月替わりに気をつけてデータを参照します。
def read_csv_skiprows(url):
resp = requests.get(url)
data_text = resp.text
is_skip, skip_row = False,1
for line in data_text.split('\n'):## 毎行分割して'data:'を検索
if 'data:' in line:
is_skip = True
break
skip_row += 1
if not is_skip:
skip_row = 0
return pd.read_csv(StringIO(data_text),skiprows=skip_row)
def get_goes_csv(goesnum,date,trange=3):
## プロットの最初の日付:tstart
if type(date)==list:
ts = datetime.date(date[0],date[1],date[2])
elif (type(date)==datetime.date)or(type(date)==datetime.datetime):
ts = date
## プロットの最後の日付:tend
te = ts+datetime.timedelta(days=trange)
url = []
if ts.month!=te.month: ## グラフが月をまたぐ場合
for i,t in enumerate([ts,te]): ## enumerateでindexも取得
## 各月の最終日を取得
last_day = calendar.monthrange(t.year,t.month)[1]
url.append('https://satdat.ngdc.noaa.gov/sem/goes/data/avg/{:04}/{:02}/goes{:02}/csv/g{:02}_xrs_1m_{:04}{:02}{:02}_{:02}{:02}{:02}.csv'.format(t.year,t.month,goesnum,goesnum,t.year,t.month,1,t.year,t.month,last_day))
if i==0:
d = read_csv_skiprows(url[i])
else: ## 今回は、2回以上月をまたぐことはないものとします
df = pd.concat([d,read_csv_skiprows(url[i])])
else:
last_day = calendar.monthrange(ts.year,ts.month)[1]
url.append('https://satdat.ngdc.noaa.gov/sem/goes/data/avg/{:04}/{:02}/goes{:02}/csv/g{:02}_xrs_1m_{:04}{:02}{:02}_{:02}{:02}{:02}.csv'\
.format(ts.year,ts.month,goesnum,goesnum,ts.year,ts.month,1,ts.year,ts.month,last_day))
df = read_csv_skiprows(url[0])
## 日付のカラムがstr型なので、時間を処理できる型に変更
df.time_tag = pd.to_datetime(df.time_tag)
## tsとteの間のデータを返す
df_obj = df[(df.time_tag>=ts)&(df.time_tag<=te)]
return df_obj
csv形式とは言ったものの少しめんどくさいことがあって、ここを見れば分かるように、欲しいデータはかなり下の方から始まっていて、最初の方はデータの説明等が記述されています。
pandasのread_csvはこういう記述があった場合に、skiprowsとかheaderの引数を指定すれば目的の行まで飛べるのですが、何が面倒かというと、ファイル(というか観測時期)によってデータに辿り着くまでの行数がバラバラです。
ただし、データの直前には'data:'という記述があるので、これを目標にスキップする行数を求めます。この方法は以前、teratailで教えていただきました。
以上の処理で、GOESのX線データを取得することができます。とりあえず自分でGOESのX線プロットを作りたければ、Sunpy1にもGOESデータを取得するコマンドがあるのでそちらを使うこともできます。
ただし、1つの時期に1つの衛星からのデータしか得られないのが欠点です(大半の場合問題にならないので欠点とは言い過ぎかもしれませんが)。例えば2017年9月6日のフレアでは、8時57分に発生したフレアのデータがちょうどGOES-15で落ちていて、GOES-13のデータで補完してやる必要があります。こういう場合のことを考えると、NOAAに取りに行った方がいいかもしれません。
3. 完成したコード
結論から書くと、以下のようなコードで出力ができました。非常に長いです。4. コード各部の解説で部分的に解説しています。とりあえずグラフが欲しいという人は5. 実際に描画するまで飛んでください。
plt.rcParams['font.family'] = 'Simplex'
plt.rcParams['font.size'] = 10
plt.rcParams['figure.dpi'] = 100
plt.rcParams['axes.linewidth'] = 1
plt.rcParams['lines.antialiased'] = False
plt.rcParams['text.antialiased'] = False
plot_antialiased = False
## 全てpixel単位
wdict = dict(
grid = 1, ## グリッドの線幅
spine = 1, ## 枠の線幅
tick = 1, ## 目盛りの線幅
line = 0.01, ## グラフの線幅
)
## pixelをpointに変換
def px2pt(px):
dpi = plt.rcParams['figure.dpi']
return px*72/dpi
def plot_goes(t, trange ,goes ,savename, wdict, width=6.4, ylim=(10**-9,10**-2)):
'''
t : プロット開始の日付 (type=datetime.date or list)
(e.g. datetime.date(2017,9,4) or [2017,9,4])
trange : プロットする期間 日単位 (type=int) (e.g. trange=3 > 3 days plot)
goes : 取得するGOES衛星の番号 list形式 (e.g. [15,13])
index = 0 がprimary (1-8A>red, 0.5-4A>blue)
index = 1 がsecondary (1-8A>yellow, 0.5-4A>purple)
width : プロットの横幅 default=6.4
if it is 'extend', extend the width with trange
ylim : y軸の範囲 default=(10**-9,10**-2)
おそらく変更の必要はないが、変えても大丈夫なように設計
'''
if type(t)==list:
t = datetime.date(t[0],t[1],t[2])
elif (type(t)==datetime.date)or(type(t)==datetime.datetime):
pass
else:
print('t should be a list or a datetime object.')
fs_def = plt.rcParams['font.size']
dpi = plt.rcParams['figure.dpi']
y0_fig = plt.rcParams['figure.subplot.bottom']
y1_fig = plt.rcParams['figure.subplot.top']
## 幅の設定
if width=='extend':
width = (dpi*6.4 + (dpi*6.4 - fs_def*10)/3*(trange-3))/dpi
elif type(width)==float:
width = width
else:
print('width is \'extend\' or float')
## 左右の余白はfontsize * 10 (figureに対する割合)
margin = fs_def/(100*width)*10
## 1日あたりの幅 (figureに対する割合)
w_day = (1-2*margin)/trange
ax = []
## 枠(axes)の作成
fig = plt.figure(figsize=(width,4.8))
for i in range(trange):
t_s = t + datetime.timedelta(days=i)
t_e = t_s + datetime.timedelta(days=1)
if i==0: ## 1日目
ax.append(fig.add_axes([margin+i*w_day,y0_fig,w_day,(y1_fig-y0_fig)],\
yscale='log', xlim=[t_s,t_e], ylim=ylim, zorder=-i))
new_xticks = date2num([t_s,t_e])
else: ## 2日目以降
ax.append(fig.add_axes([margin+i*w_day,y0_fig,w_day,(y1_fig-y0_fig)],\
yscale='log', xlim=[t_s,t_e], ylim=ylim, yticklabels=[], zorder=-i))
new_xticks = date2num([t_e])
ax[i].xaxis.set_major_locator(ticker.FixedLocator(new_xticks))
ax[i].xaxis.set_major_formatter(DateFormatter('%b %d'))
ax[i].xaxis.set_minor_locator(matplotlib.dates.HourLocator(interval=3))
ax[i].yaxis.grid(lw=px2pt(wdict['grid']),color='black')
ax[i].tick_params(axis='x', which='major', bottom=True, top=True, direction='in', length=6, width=px2pt(wdict['tick']), pad=fs_def)
ax[i].tick_params(axis='x', which='minor', bottom=True, top=True, direction='in', length=3, width=px2pt(wdict['tick']))
for j in ax[i].spines.keys(): ## 枠の線幅の調整
ax[i].spines[j].set_linewidth(px2pt(wdict['spine']))
## 目盛りの調整
if i==0: ## 1日目
ax[i].tick_params(axis='y', which='both', right=True, left=True, direction='in', length=6, width=px2pt(wdict['tick']))
ax[i].spines['right'].set_visible(False)
elif i==trange-1: ## 最終日
ax[i].tick_params(axis='y', which='both', right=True, left=False, direction='in', length=6, width=px2pt(wdict['tick']))
ax[i].spines['left'].set_visible(False)
else: ## その他
ax[i].tick_params(axis='y', which='both', right=True, left=False, direction='in', length=6, width=px2pt(wdict['tick']))
ax[i].spines['left'].set_visible(False)
ax[i].spines['right'].set_visible(False)
## Flare class text
f_pos = ax[-1].get_position()
trans = transforms.blended_transform_factory(fig.transFigure,transforms.IdentityTransform())
X10_y_display = ax[-1].transData.transform((date2num(t_e),10**-3))[1]
X_y_display = ax[-1].transData.transform((date2num(t_e),10**-4))[1]
w = X10_y_display - X_y_display
## 図全体に対するフォントサイズの相対的な大きさ
## クラスのテキストを置く際に、外枠に張り付かないよう余白を持たせるために使用
cls_text_size = plt.rcParams['font.size']/(plt.rcParams['figure.figsize'][0])/72
for i,cls in enumerate(['X', 'M', 'C', 'B', 'A']):
pos = [f_pos.x1,X10_y_display - w/2 - i*w]
if (pos[1]>=ax[-1].transData.transform((date2num(t_e),ax[-1].get_ylim()[0]))[1])&\
(pos[1]<=ax[-1].transData.transform((date2num(t_e),ax[-1].get_ylim()[1]))[1]):
ax[-1].text(pos[0]+cls_text_size/5,pos[1],cls,horizontalalignment='left',verticalalignment='center',transform=trans)
## GOES number and wavelength
ybox_b = TextArea('GOES{} 0.5-4.0 A'.format(goes[0]), textprops=dict(color='#001DFF', size='large' ,rotation=90,ha='left',va='bottom')) #青
ybox_r = TextArea('GOES{} 1.0-8.0 A'.format(goes[0]), textprops=dict(color='#FF0000', size='large' ,rotation=90,ha='left',va='bottom')) #赤
ybox_primary = VPacker(children=[ybox_r, ybox_b],align='bottom', pad=0, sep=2*fs_def)
if len(goes)==2:
ybox_p = TextArea('GOES{} 0.5-4.0 A'.format(goes[1]), textprops=dict(color='#5400A8', size='large' ,rotation=90,ha='left',va='bottom')) #紫
ybox_y = TextArea('GOES{} 1.0-8.0 A'.format(goes[1]), textprops=dict(color='#FFA500', size='large' ,rotation=90,ha='left',va='bottom')) #黄
ybox_secondary = VPacker(children=[ybox_y, ybox_p],align='bottom', pad=0, sep=2*fs_def)
ybox = HPacker(children=[ybox_primary,ybox_secondary],align='bottom', pad=0, sep=fs_def)
anchored_ybox = AnchoredOffsetbox(loc=5, child=ybox, pad=0., frameon=False, bbox_to_anchor=(1-margin*0.25, 0.5),
bbox_transform=fig.transFigure, borderpad=0.)
else:
ybox = ybox_primary
anchored_ybox = AnchoredOffsetbox(loc=5, child=ybox, pad=0., frameon=False, bbox_to_anchor=(1-margin*0.5, 0.5),
bbox_transform=fig.transFigure, borderpad=0.)
ax[-1].add_artist(anchored_ybox)
## 軸ラベル作成 (dummy axes)
ax_dummy = fig.add_axes([margin,y0_fig,1-margin*2,y1_fig-y0_fig],frameon=False)
ax_dummy.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
ax_dummy.set_xlabel('Universal Time',size='large')
ax_dummy.set_ylabel('Watts m$^{-2}$',size='large')
ax_dummy.set_title('GOES Xray Flux (1-minute data)',size='x-large')
ax_dummy.text(1, 1-cls_text_size/5, 'Begin: {} {} {} 0000 UTC'.format(t.year,t.strftime('%b'),t.day),
horizontalalignment='right',
verticalalignment='top',
transform=fig.transFigure)
## X-ray data plot
color = iter(['#FF0000','#001DFF','#FFA500','#5400A8']) ## red, blue, yellow, purple
zo = iter([4,3,2,1])
for g in goes:
wavelength = iter(['1.0-8.0 A','0.5-4.0 A'])
data = get_goes_csv(g,t,trange)
xray_long = np.ma.masked_where((data.B_AVG<=ylim[0])|(data.B_AVG>=ylim[1]),data.B_AVG)
xray_short = np.ma.masked_where((data.A_AVG<=ylim[0])|(data.A_AVG>=ylim[1]),data.A_AVG)
ax[0].plot(data.time_tag,xray_long, lw=px2pt(wdict['line']),c=next(color),label='GOES{} 1.0-8.0 A'.format(g),aa=plot_antialiased,clip_on=False,zorder=next(zo))
ax[0].plot(data.time_tag,xray_short,lw=px2pt(wdict['line']),c=next(color),label='GOES{} 0.5-4.0 A'.format(g),aa=plot_antialiased,clip_on=False,zorder=next(zo))
## save graph
fig.savefig(savename)
4. コード各部の解説
各部分の紹介をしていきます。とりあえず使えるようになりたい、という方は読み飛ばしてください。
4.1. フォント
IDLで使用されているフォントは、こちらによればSimplex Romanというフォントだそうです。今回は見た目にこだわったため、フリーのフォントサイトからダウンロードして使いました。Mac OSには標準では搭載されていなかったため、こだわりたい方はダウンロードしてください。
必要なければ、
plt.rcParams['font.family'] = 'Arial'
などとしておけば見やすいグラフになると思います。
4.2. エイリアシング
plt.rcParams['lines.antialiased'] = False
plt.rcParams['text.antialiased'] = False
plot_antialiased = False
matplotlibで、pngやjpgなどのビットマップ形式で出力する場合、自動で補完をしてギザギザがなくなるよう、アンチエイリアシングを行ってくれます。
普通に図を書くならここを気にする必要はないのですが、今回はこだわりました。上記のパラメータを全てTrueにすればアンチエイリアシングを行いますし、その方が無難だと思います。なお、tickのアンチエイリアシングを解除する方法は、どうやらないようです。
4.3. pixelとpointの変換
def px2pt(px):
dpi = plt.rcParams['figure.dpi']
return px*72/dpi
これも無駄なこだわりです。元画像では、あらゆる線の幅が1 pixelに指定されていたので、それのために作成しました。というのも、matplotlibで幅などを指定する際の単位は全てpointです。1 pointは1 inchの1/72で、inchとpixelの関係はdpiによって変化するので、それを補正しています。
4.4. 幅の設定
if width=='extend':
width = (dpi*6.4 + (dpi*6.4 - fs_def*10)/3*(trange-3))/dpi
elif type(width)==float:
width = width
else:
print('width is \'extend\' or float')
元画像は縦が480pixel、横が640pixelだったので、基本はそれに合わせています(dpiは100)。ただし、3日よりも長いプロットを作成する場合は、width=='extend'と指定すれば、1日あたりの幅を3days plotの時と同じにしたまま、横のサイズを伸ばせるようにしました。
また、普通にfloat型で幅を指定すればそのサイズになるように設定していますが、動作保証はできません。
なお、よくよく考えたら3日より短いplotを作成した場合はエラーを吐きますので、後で記述します(6. シンプルに描画する)。
4.5. 余白の設定
## 左右の余白はfontsize * 10 (figureに対する割合)
margin = fs_def/(100*width)*10
## 1日あたりの幅 (figureに対する割合)
w_day = (1-2*margin)/trange
左右の余白を設定します。コメントアウトにある通り、フォントサイズのおよそ10倍に設定しています。
4.6. 枠の設定
## 枠(axes)の作成
fig = plt.figure(figsize=(width,4.8))
for i in range(trange):
t_s = t + datetime.timedelta(days=i)
t_e = t_s + datetime.timedelta(days=1)
if i==0: ## 1日目
ax.append(fig.add_axes([margin+i*w_day,y0_fig,w_day,(y1_fig-y0_fig)],\
yscale='log', xlim=[t_s,t_e], ylim=ylim, zorder=-i))
new_xticks = date2num([t_s,t_e])
else: ## 2日目以降
ax.append(fig.add_axes([margin+i*w_day,y0_fig,w_day,(y1_fig-y0_fig)],\
yscale='log', xlim=[t_s,t_e], ylim=ylim, yticklabels=[], zorder=-i))
new_xticks = date2num([t_e])
ax[i].xaxis.set_major_locator(ticker.FixedLocator(new_xticks))
ax[i].xaxis.set_major_formatter(DateFormatter('%b %d'))
ax[i].xaxis.set_minor_locator(matplotlib.dates.HourLocator(interval=3))
ax[i].yaxis.grid(lw=px2pt(grid_width),color='black')
ax[i].tick_params(axis='x', which='major', bottom=True, top=True, direction='in', length=6, width=px2pt(tick_width), pad=fs_def)
ax[i].tick_params(axis='x', which='minor', bottom=True, top=True, direction='in', length=3, width=px2pt(tick_width))
for j in ax[i].spines.keys(): ## 枠の線幅の調整
ax[i].spines[j].set_linewidth(px2pt(spine_width))
## 目盛りの調整
if i==0: ## 1日目
ax[i].tick_params(axis='y', which='both', right=True, left=True, direction='in', length=6, width=px2pt(tick_width))
ax[i].spines['right'].set_visible(False)
elif i==trange-1: ## 最終日
ax[i].tick_params(axis='y', which='both', right=True, left=False, direction='in', length=6, width=px2pt(tick_width))
ax[i].spines['left'].set_visible(False)
else: ## その他
ax[i].tick_params(axis='y', which='both', right=True, left=False, direction='in', length=6, width=px2pt(tick_width))
ax[i].spines['left'].set_visible(False)
ax[i].spines['right'].set_visible(False)
多分一番細かいところに一番時間をかけてこだわった部分です。というのも、元のグラフの内部にあるメモリ、これを描画する手段はどうやらmatplotlibにはないようです。
ax.hline()などで範囲を区切った横棒を無理やり書くことも考えましたが、あまり力技はやりたくありませんでした。それで、IDLを扱える同期に聞いたところ、IDLには任意の場所に軸を追加する機能があるそうです。何それmatplotlibにもくれ。
仕方がないので、axesを分割することで対処することにしました。具体的には、add_axesで座標を指定することで以下の図のように領域を分割し、枠線(matplotlibではspineと呼ぶそうです)の表示を変えることで、無事に枠の作成が完了しました。
4.6. フレアクラスのテキスト
## Flare class text
## 左端のaxesの座標を取得
f_pos = ax[-1].get_position()
trans = transforms.blended_transform_factory(fig.transFigure,transforms.IdentityTransform())
X10_y_display = ax[-1].transData.transform((date2num(t_e),10**-3))[1]
X_y_display = ax[-1].transData.transform((date2num(t_e),10**-4))[1]
w = X10_y_display - X_y_display
## 図全体に対するフォントサイズの相対的な大きさ
## クラスのテキストを置く際に、外枠に張り付かないよう余白を持たせるために使用
cls_text_size = plt.rcParams['font.size']/(plt.rcParams['figure.figsize'][0])/72
for i,cls in enumerate(['X', 'M', 'C', 'B', 'A']):
pos = [f_pos.x1,X10_y_display - w/2 - i*w]
if (pos[1]>=ax[-1].transData.transform((date2num(t_e),ax[-1].get_ylim()[0]))[1])&\
(pos[1]<=ax[-1].transData.transform((date2num(t_e),ax[-1].get_ylim()[1]))[1]):
ax[-1].text(pos[0]+cls_text_size/5,pos[1],cls,horizontalalignment='left',verticalalignment='center',transform=trans)
このコードの中で一番可読性の低い部分です。すみません。
フレアのクラスは縦軸の値によって決まるので、正しい値のところにaxes.textで文字を置いてやる必要があります。
axes.textを置く時の座標の設定の仕方はいくつかあって、
- 図全体のピクセルに対応する座標
- 図の左下を(0,0)、右上を(1,1)に定めた相対座標
- 図の中にあるaxes(枠)の中で、左下を(0,0)、右上を(1,1)に定めた相対座標
- 図が持っている値に対応する座標
なのですが、この図は横軸が時刻なので、(4)を諦めました(date2numを使えば可能？)。そのため、(1)〜(3)を混ぜて扱えるblended_transform_factoryというものを使いながら、枠線の外側の正しい位置に文字を置けるように設計しました。具体的には、
trans = transforms.blended_transform_factory(fig.transFigure,transforms.IdentityTransform())
とすることでx座標は(2; figure.transFigure)を、y座標は(4; transforms.IdentityTransform())を参照して、textを置けるようにしています。
ここはめちゃくちゃ改善点ありそうです...
4.7. レジェンドの作成
## GOES number and wavelength
## 各色のテキストを90度回転させて設定
ybox_b = TextArea('GOES{} 0.5-4.0 A'.format(goes[0]), textprops=dict(color='#001DFF', size='large' ,rotation=90,ha='left',va='bottom')) #青
ybox_r = TextArea('GOES{} 1.0-8.0 A'.format(goes[0]), textprops=dict(color='#FF0000', size='large' ,rotation=90,ha='left',va='bottom')) #赤
## VPackerで縦方向に連結
ybox_primary = VPacker(children=[ybox_r, ybox_b],align='bottom', pad=0, sep=2*fs_def)
## primaryとsecondaryの両方をプロットする場合は、
## secondaryも同様に作成して HPackerで横方向に連結
if len(goes)==2:
ybox_p = TextArea('GOES{} 0.5-4.0 A'.format(goes[1]), textprops=dict(color='#5400A8', size='large' ,rotation=90,ha='left',va='bottom')) #紫
ybox_y = TextArea('GOES{} 1.0-8.0 A'.format(goes[1]), textprops=dict(color='#FFA500', size='large' ,rotation=90,ha='left',va='bottom')) #黄
ybox_secondary = VPacker(children=[ybox_y, ybox_p],align='bottom', pad=0, sep=2*fs_def)
ybox = HPacker(children=[ybox_primary,ybox_secondary],align='bottom', pad=0, sep=fs_def)
## AnchoredOffsetboxで設置 locの番号ごとに設置する場所が決められている
## loc=5 は'center right' つまり図の右中央
anchored_ybox = AnchoredOffsetbox(loc=5, child=ybox, pad=0., frameon=False, bbox_to_anchor=(1-margin*0.25, 0.5), bbox_transform=fig.transFigure, borderpad=0.)
else: ## primaryだけの場合はprimaryだけで
ybox = ybox_primary
anchored_ybox = AnchoredOffsetbox(loc=5, child=ybox, pad=0., frameon=False, bbox_to_anchor=(1-margin*0.5, 0.5), bbox_transform=fig.transFigure, borderpad=0.)
ax[-1].add_artist(anchored_ybox)
ここも可読性低めです。
元画像ではレジェンドが枠の右外側に90度回転して置いてあります。matplotlibでは、レジェンドを枠外に置くことはできても、回転させる手段はどうやらないようです。さらに、axes.text中で色を変えることもできないようなので、色の異なるtextを別々に配置する必要があります。これは、単なる座標指定を用いるだけでは非常に困難です。そのため、相対位置で配置が可能なmatplotlib.offsetboxのVPackerやHPackerを用いて、図の外に回転したレジェンドを再現しました。かなり力技の手法です。
まず、TextAreaで各々のtextを生成し、これらをVPackerやHPackerを用いて連結していきます。そうして完成した擬似レジェンドを、AnchoredOffsetboxを用いて図右端に設置しました。
locの数字ごとに設置する場所が決まっており、loc=5はcenter rightで右端に設置することを意味します。
4.8. 軸ラベルの作成
## 軸ラベル作成 (dummy axes)
ax_dummy = fig.add_axes([margin,y0_fig,1-margin*2,y1_fig-y0_fig],frameon=False)
ax_dummy.tick_params(labelcolor='none', top='off', bottom='off', left='off', right='off')
ax_dummy.set_xlabel('Universal Time',size='large')
ax_dummy.set_ylabel('Watts m$^{-2}$',size='large')
ax_dummy.set_title('GOES Xray Flux (1-minute data)',size='x-large')
ax_dummy.text(1, 1-cls_text_size/5, 'Begin: {} {} {} 0000 UTC'.format(t.year,t.strftime('%b'),t.day),
horizontalalignment='right',
verticalalignment='top',
transform=fig.transFigure)
図の内側に目盛を表示するためにaxesを分けて表示するという方法~~(荒技)~~を使いましたが、そのおかげで横軸のラベル設定位置が定まらなくなってしまいました。そこで、何も表示されない大きなダミーのaxesを作り、それに対して軸ラベルを設定しています。
また、元画像の上下左右に表示されている文字群の詳細な再現は(特に下部の文字は必要がないため)諦めました。ただ、右上にあるデータの開始時刻がないと、データの年の情報が分からなくなるため、この部分は急遽表示することにしました。
4.9. X-rayデータのプロット
## X-ray data plot
color = iter(['#FF0000','#001DFF','#FFA500','#5400A8']) ## red, blue, yellow, purple
zo = iter([4,3,2,1])
for g in goes:
wavelength = iter(['1.0-8.0 A','0.5-4.0 A'])
data = get_goes_csv(g,t,trange)
xray_long = np.ma.masked_where((data.B_AVG<=ylim[0])|(data.B_AVG>=ylim[1]),data.B_AVG)
xray_short = np.ma.masked_where((data.A_AVG<=ylim[0])|(data.A_AVG>=ylim[1]),data.A_AVG)
ax[0].plot(data.time_tag,xray_long, lw=px2pt(line_width),c=next(color),label='GOES{} 1.0-8.0 A'.format(g),aa=plot_antialiased,clip_on=False,zorder=next(zo))
ax[0].plot(data.time_tag,xray_short,lw=px2pt(line_width),c=next(color),label='GOES{} 0.5-4.0 A'.format(g),aa=plot_antialiased,clip_on=False,zorder=next(zo))
ようやくX線データのプロットまで来ました。
GOESの1分平均データのcsvファイルのカラムのうち、X線フラックスをプロットする場合に必要なのは、'time_tag'と'A_AVG'と'B_AVG'です。'time_tag'はそのまま毎分の時刻が、'A_AVG'と'B_AVG'はそれぞれ、0.05 - 0.4 nmと0.1-0.8 nmに対応するX線の平均データが格納されています。これらを、pandasのDataFrameの文法に基づいて取り出しています。
zorderは描画の順番です。この値が大きいものほど前面で表示されます。元画像ではGOES primaryの値が前面に表示されていたため、そうなるように数字を設定しました。
その他に行なっていることとして、プロット時にclip_on=Falseを指定しています。今回は日付分のaxesを用意しているので、本来ならその数だけaxes.plotを呼び出す必要がありますが、clip_on=Falseを指定すると用意されたaxesからはみ出してプロットができるので、1日目のaxesで全日分の時間を指定したデータをプロットしています。更にこのclip_on=Falseを行うことで、枠線の上にもプロットを乗せることができています。
加えて、np.ma.masked_whereを用いた値のマスクをしています。今回使用するデータは、データがロストしているときは-99999が格納されており、これをそのままログスケールでプロットしてしまうと、井戸型ポテンシャルみたいな図が出来上がってしまいます。
さらに、デフォルトで設定したylim=(10**-9,10**-2)の範囲で描画するなら基本的に問題はないですが、仮にここを変更した場合、clip_on=Falseを設定している影響で上下方向にもはみ出してプロットされる可能性があります。そのため、
xray_long = np.ma.masked_where((data.B_AVG<=ylim[0])|(data.B_AVG>=ylim[1]),data.B_AVG)
のようにして、-99999を含め、設定したy軸の範囲から出るデータを全てマスクするようにしています。こうすることで、マスクされた部分のプロットは繋がることなく、空白になります。
通常、np.nanを格納することでもプロットは空白になるため、あまりnumpy.maを使う例は出てこないように思いますが、全く別のグラフを描画する際にnp.nanの格納ではなくマスクをかけた方が都合がいいことがありました。そのため、ここでは紹介のつもりでnumpy.maを使用しています。
5. 実際に描画する
注:Simplexのフォントはデフォルトで各PCに用意されていない可能性があります。
上記の関数を使って描画していきます。基本的に、必要なパラメータは最初に宣言していますが、
plt.rcParams['font.family'] = 'Simplex'
plt.rcParams['font.size'] = 10
plt.rcParams['figure.dpi'] = 100
plt.rcParams['axes.linewidth'] = 1
plt.rcParams['lines.antialiased'] = False
plt.rcParams['text.antialiased'] = False
plot_antialiased = False
## 全てpixel単位
wdict = dict(
grid = 1, ## グリッドの線幅
spine = 1, ## 枠の線幅
tick = 1, ## 目盛りの線幅
line = 0.01, ## グラフの線幅
)
plot_goes([2017,9,4], 3 ,[15,13] ,'goes.png', wdict, width=6.4, ylim=(10**-9,10**-2))
のようにすればプロットできます。
並べて比較するとこのような感じ。
結構再現できていませんか？？ダウンロードしたSimplexフォントがかなり細いため、エイリアシングの関係で文字の一部分が表示されないことと、上下左右の文字情報以外はほぼ再現できたのではないかと思います。~~枠のサイズとかはもうあと手作業でチューニングするしかなさそうだしパス。~~これにてコピー完了です。
各入力を少し変えても問題なく描画されるはずです。
plot_goes([2017,9,4], 5 ,[15] ,'goes_extend.png', wdict, width='extend', ylim=(10**-9,10**-2))
6. シンプルに描画する
5. 実際に描画するまでで、本家の再現はほぼ完了しました。ここからは、もう少し実用的な(シンプルな)図にします。
上記のgoes_plot()のコードが煩雑になっている要素のトップ2は、
- 図中に目盛を表示しようとしたこと
- レジェンドを90度回転させて置こうとしたこと
だと思います。ポスターや発表スライドでそこまでこだわる必要はないと思いますので、これらを排除してコードの簡素化を目指します。加えて、上記の関数goes_plot()はこれを実行するたびにサーバーにデータを取りに行くので、時間がかかるしアクセス過多になりかねません。そこで、一度ダウンロードしたデータを使いまわして、図の細部を調整できるように一部改造します。(当然ですが、データファイルをローカルにダウンロードして、それを読み込むのが一番です。)
6.1. コード
plt.rcParams['font.family'] = 'Arial'
plt.rcParams['font.size'] = 11
plt.rcParams['figure.dpi'] = 100
plt.rcParams['axes.linewidth'] = 1
plt.rcParams['lines.antialiased'] = True
plt.rcParams['text.antialiased'] = True
plot_antialiased = True
## 全てpixel単位
wdict = dict(
grid = 1, ## グリッドの線幅
spine = 1, ## 枠の線幅
tick = 1, ## 目盛りの線幅
line = 2, ## グラフの線幅
)
## pixelをpointに変換
def px2pt(px):
dpi = plt.rcParams['figure.dpi']
return px*72/dpi
def plot_goes_simple(t, trange ,goes ,savename, wdict, width=6.4, ylim=(10**-9,10**-2)):
'''
t : プロット開始の日付 (type=datetime.date or list)
(e.g. datetime.date(2017,9,4) or [2017,9,4])
trange : プロットする期間 日単位 (type=int) (e.g. trange=3 > 3 days plot)
goes : 取得するGOES衛星の番号 list形式 (e.g. [15,13])
index = 0 がprimary (1-8A>red, 0.5-4A>blue)
index = 1 がsecondary (1-8A>yellow, 0.5-4A>purple)
width : プロットの横幅 default=6.4
if it is 'extend', extend the width with trange
ylim : y軸の範囲 default=(10**-9,10**-2)
おそらく変更の必要はないが、変えても大丈夫なように設計
'''
if type(t)==list:
t = datetime.date(t[0],t[1],t[2])
elif (type(t)==datetime.date)or(type(t)==datetime.datetime):
pass
else:
print('t should be a list or a datetime object.')
fs_def = plt.rcParams['font.size']
dpi = plt.rcParams['figure.dpi']
y0_fig = plt.rcParams['figure.subplot.bottom']
y1_fig = plt.rcParams['figure.subplot.top']
## 幅の設定
if width=='extend':
width = (dpi*6.4 + (dpi*6.4 - fs_def*10)/3*(trange-3))/dpi
elif type(width)==float:
width = width
else:
print('width is \'extend\' or float')
## 左右の余白はfontsize * 10 (figureに対する割合)
margin = fs_def/(100*width)*10
## 1日あたりの幅 (figureに対する割合)
w_day = (1-2*margin)/trange
ax = []
## 枠(axes)の作成
t_s = t
t_e = t_s + datetime.timedelta(days=trange)
fig = plt.figure(figsize=(width,4.8))
## 枠内に通常のレジェンドを置くなら、無理にadd_axesで指定することなく、
## デフォルトのfig.add_subplotでいいと思います
#ax = fig.add_axes([margin,y0_fig,1-2*margin,y1_fig-y0_fig],yscale='log', xlim=[t_s,t_e], ylim=ylim)
ax = fig.add_subplot(111,yscale='log', xlim=[t_s,t_e], ylim=ylim)
new_xticks = date2num([t_s+datetime.timedelta(days=i) for i in range(trange+1)])
ax.xaxis.set_major_locator(ticker.FixedLocator(new_xticks))
ax.xaxis.set_major_formatter(DateFormatter('%b %d'))
ax.xaxis.set_minor_locator(matplotlib.dates.HourLocator(interval=3))
ax.yaxis.grid(lw=px2pt(wdict['grid']),color='black')
ax.tick_params(axis='x', which='major', bottom=True, top=True, direction='in', length=6, width=px2pt(wdict['tick']), pad=fs_def)
ax.tick_params(axis='x', which='minor', bottom=True, top=True, direction='in', length=3, width=px2pt(wdict['tick']))
ax.tick_params(axis='y', which='both', right=True, left=True, direction='in', length=6, width=px2pt(wdict['tick']))
for j in ax.spines.keys(): ## 枠の線幅の調整
ax.spines[j].set_linewidth(px2pt(wdict['spine']))
## Flare class text
f_pos = ax.get_position()
trans = transforms.blended_transform_factory(fig.transFigure,transforms.IdentityTransform())
X10_y_display = ax.transData.transform((date2num(t_e),10**-3))[1]
X_y_display = ax.transData.transform((date2num(t_e),10**-4))[1]
w = X10_y_display - X_y_display
## 図全体に対するフォントサイズの相対的な大きさ
## クラスのテキストを置く際に、外枠に張り付かないよう余白を持たせるために使用
cls_text_size = plt.rcParams['font.size']/(plt.rcParams['figure.figsize'][0])/72
for i,cls in enumerate(['X', 'M', 'C', 'B', 'A']):
pos = [f_pos.x1,X10_y_display - w/2 - i*w]
if (pos[1]>=ax.transData.transform((date2num(t_e),ax.get_ylim()[0]))[1])&\
(pos[1]<=ax.transData.transform((date2num(t_e),ax.get_ylim()[1]))[1]):
ax.text(pos[0]+cls_text_size/2,pos[1],cls,horizontalalignment='center',verticalalignment='center',transform=trans)
## ここのコメントアウトを外せば、枠外にレジェンドを置けます
'''
## GOES number and wavelength
ybox_b = TextArea('GOES{} 0.5-4.0 A'.format(goes[0]), textprops=dict(color='#001DFF', size='large' ,rotation=90,ha='left',va='bottom')) #青
ybox_r = TextArea('GOES{} 1.0-8.0 A'.format(goes[0]), textprops=dict(color='#FF0000', size='large' ,rotation=90,ha='left',va='bottom')) #赤
ybox_primary = VPacker(children=[ybox_r, ybox_b],align='bottom', pad=0, sep=2*fs_def)
if len(goes)==2:
ybox_p = TextArea('GOES{} 0.5-4.0 A'.format(goes[1]), textprops=dict(color='#5400A8', size='large' ,rotation=90,ha='left',va='bottom')) #紫
ybox_y = TextArea('GOES{} 1.0-8.0 A'.format(goes[1]), textprops=dict(color='#FFA500', size='large' ,rotation=90,ha='left',va='bottom')) #黄
ybox_secondary = VPacker(children=[ybox_y, ybox_p],align='bottom', pad=0, sep=2*fs_def)
ybox = HPacker(children=[ybox_primary,ybox_secondary],align='bottom', pad=0, sep=fs_def)
anchored_ybox = AnchoredOffsetbox(loc=5, child=ybox, pad=0., frameon=False, bbox_to_anchor=(1-margin*0.25, 0.5),
bbox_transform=fig.transFigure, borderpad=0.)
else:
ybox = ybox_primary
anchored_ybox = AnchoredOffsetbox(loc=5, child=ybox, pad=0., frameon=False, bbox_to_anchor=(1-margin*0.5, 0.5),
bbox_transform=fig.transFigure, borderpad=0.)
ax.add_artist(anchored_ybox)
'''
## 軸ラベル作成
ax.set_xlabel('Universal Time',size='large')
ax.set_ylabel('Watts m$^{-2}$',size='large')
ax.set_title('GOES Xray Flux (1-minute data)',size='x-large')
ax.text(1, 1-cls_text_size/5, 'Begin: {} {} {} 0000 UTC'.format(t.year,t.strftime('%b'),t.day),
horizontalalignment='right',
verticalalignment='top',
transform=fig.transFigure)
## X-ray data plot
color = iter(['#FF0000','#001DFF','#FFA500','#5400A8']) ## red, blue, yellow, purple
zo = iter([4,3,2,1])
for i,g in enumerate(goes):
data = data_temp[i]
xray_long = np.ma.masked_where((data.B_AVG<=ylim[0])|(data.B_AVG>=ylim[1]),data.B_AVG)
xray_short = np.ma.masked_where((data.A_AVG<=ylim[0])|(data.A_AVG>=ylim[1]),data.A_AVG)
ax.plot(data.time_tag,xray_long, lw=px2pt(wdict['line']),c=next(color),label='GOES{} 1.0-8.0 A'.format(g),aa=plot_antialiased,clip_on=False,zorder=next(zo))
ax.plot(data.time_tag,xray_short,lw=px2pt(wdict['line']),c=next(color),label='GOES{} 0.5-4.0 A'.format(g),aa=plot_antialiased,clip_on=False,zorder=next(zo))
## simple legend
leg = ax.legend(handlelength=0, handletextpad=0, ncol=2)
for line, text in zip(leg.get_lines(), leg.get_texts()):
text.set_color(line.get_color())
## save graph
fig.savefig(savename)
axesを複数置く必要がなくなったので、変数axまわりの表記が若干変わっています。アレンジする場合はplot_goes()との違いに注意してください。
## simple legend
leg = ax.legend(handlelength=0, handletextpad=0, ncol=2)
for line, text in zip(leg.get_lines(), leg.get_texts()):
text.set_color(line.get_color())
で、レジェンドのマーカーを見えなくし、レジェンドテキストの色を各プロットに対応した色に変更しています。
6.2. 描画
いくつか試しましたが、以下のパラメータがちょうど良さそうに思います。
plt.rcParams['font.family'] = 'Arial'
plt.rcParams['font.size'] = 11
plt.rcParams['figure.dpi'] = 100
plt.rcParams['axes.linewidth'] = 1
plt.rcParams['lines.antialiased'] = True
plt.rcParams['text.antialiased'] = True
plot_antialiased = True
## 全てpixel単位
wdict = dict(
grid = 1, ## グリッドの線幅
spine = 1, ## 枠の線幅
tick = 1, ## 目盛りの線幅
line = 2, ## グラフの線幅
)
以下のようにしてdata_tempにX線データを格納し、プロットします。
t = [2017,9,4]
trange = 3
goes = [15,13]
data_temp = []
for g in goes:
data_temp.append(get_goes_csv(g,t,trange))
plot_goes_simple(t, trange, goes, 'goes_simple.png', wdict, width=6.4, ylim=(10**-9,10**-2))
こんな感じなら、ポスターに載せても比較的見やすいと思います。
7. おわりに
こだわりだしたら止まらなくなってしまった。axes周りのことはめちゃくちゃ勉強になりました。もし機会があれば、論文や発表スライド、ポスターなどに活用してください。万が一役に立ったら、どこかでお会いした時にでも「役に立ったよ」と言ってくださると、「やったぜ」って気持ちになって大変嬉しいです。
質問や改善点などあればぜひご連絡ください。GOES-16はjson形式でリアルタイムデータが公開されていますが、直近1週間のデータしか置いていないんですよね。過去のデータはそのうち公開されるのだろうか。~~されないとそのうち困るんですが。~~またフォーマットが定まれば、対応できるように作ってみたいと思います。
参考文献
matplotlib全般
matplotlib.figure.Figure — Matplotlib 3.1.2 documentation
早く知っておきたかったmatplotlibの基礎知識、あるいは見た目の調整が捗るArtistの話 - Qiita
matplotlibのめっちゃまとめ - Qiita
matplotlib逆引きメモ（簡単なグラフ） - Qiita
matplotlib逆引きメモ（フレーム編） - Qiita
tick, spineまわり
matplotlib.datesで時系列データのグラフの軸目盛の設定をする | 分析トレイン
Matplotlib 目盛と目盛ラベル、目盛線
matplotlibで枠線を消したグラフをつくる - Qiita
text配置と座標指定(transform)
Text properties and layout — Matplotlib 3.1.2 documentation
Transformations Tutorial — Matplotlib 3.1.2 documentation
複数subplotに対する軸ラベルの設置(ダミーaxesの設定)
python - pyplot axes labels for subplots - Stack Overflow
legend markerとlegend text colorの調整
python - How do you just show the text label in plot legend? (e.g. remove a label's line in the legend) - Stack Overflow
Option to set the text color in legend to be same as the line · Issue #10720 · matplotlib/matplotlib · GitHub
90度回転したlegendの作成(matplotlib.offsetbox)
python - Matplotlib: y-axis label with multiple colors - Stack Overflow
matplotlib.offsetbox — Matplotlib 3.1.2 documentation
numpy.ma
The numpy.ma module — NumPy v1.17 Manual
-
この部分、あえて古いバージョンのdocumentを参照しています。というのも、最新版のdocumentにはなぜか図が表示されていなくて直感的にわかりづらいからです。 ↩