プログラミング系科目オール「C」の俺がPythonに一目惚れして局所天気図の描写に挑む話 ～第3話～

あらすじ(と謝罪)

大学の必修科目であるプログラミング(C言語)を機にプログラミングに嫌悪感を抱くようになったうえはーす。しかし別の科目で触れたPythonに一目惚れすると、局地天気図の描写に挑戦するほどのめり込む。一度は書けた天気図だが気に食わず、現在コード見直し中。いよいよ問題の等圧線に取り掛かる。

そして前回の投稿からかなり日が開いてしまったことをお詫びしたい。なぜなら、大学4年生であるがゆえに就活や卒研が本格化してしまって記事作成に割ける時間が少なくなってしまったのだ。今後もペースは確実に落ち込むだろうが、どうか温かい目で見守っていただけるとありがたい。

理想とする天気図の完成形

前回利用した地点名表示のコードを応用し、気圧を各プロット近傍に表示してやる。表示したい時間は熱的低気圧と雷雨性高気圧が同時に現れている様子がうかがえる17時にしたいので
df.iloc[i]['site']→df.iloc[i]['JST17']
に変える。

# 観測点気圧表示
for i in range(len(df)):
    ax.text(df.iloc[i]['lon'], df.iloc[i]['lat'], df.iloc[i]['JST17'], 
            ha='center', va='top', color='navy', fontweight='heavy', fontsize='small', zorder=7)

得られる図はこうなる。

熱的低気圧の顕現により松本付近の海面更生気圧が低い一方で、積乱雲のほぼ直下にあたる諏訪で雷雨性高気圧(メソハイ)が顕現し、海面更生気圧が高くなっている。

この図を基にすると、俺が今描かんとする理想の天気図はこんな感じである。マウスを使って描いているので線が汚いのは容赦してほしい。

1hPaごとに線を引いている

「それらしい線」で一筋縄ではいかなかった等圧線描写

等圧線の描写を描く上で数々の有志によるサイトを渡り歩いたが、簡単に言えば「あらかじめ連続値データが揃っている状態で、それをなぞる」方法が大半だった。

なので俺がやろうとしている「測定箇所の座標が散らばっている気圧データを補間して等圧線を描く」やり方を見つけるにはかなり苦戦した。

とりあえず下に示すリンク先の情報が一番手助けになったので紹介。そしてこの場を借りて感謝させていただきたい。

値を持つ座標群からヒートマップを表示 - Qiita

先に示したリンク先にて「scipy.interpolate.griddata」の扱い方がやさしく解説されているが、matplotlib.axes.Axes.textで引数の型を誤って与える前科(第2話参照)のある俺は本家のマニュアルを確認しなくてはならない。

scipy.interpolate.griddata — SciPy v1.12.0 Manual

座標にあたるpointsには浮動小数点(float)型の2次元配列を与えてやればいいそうだ。lonとlatの値をそのまま与えればいいだろう。

その座標における数値であるvaluesにはN次元配列(ndarray)、浮動小数点(float)型、複素数(complex)型、各次元の要素数を受け渡すshapeが対応しているようだ。ここはdf['JST17']渡しておけば問題ないか。

xiは補間対象とする座標(=2次元配列)を入れてやるそうだ。lonとlatをmeshgridで格子点に変換した座標を与えてみよう。

methodは補間の方法で、nearestで最も近いpointsのvaluesを返し、linearは線形補間、cubicは曲線補間。ちょっと何言ってるか分かんないが、とりあえずやってみれば何してるかわかるだろう。ちなみにデフォルトはlinearだそうだ。

あとは引き続き引数の型に気を付けながら「matplotlib.axes.Axes.contour」と組み合わせてこんな感じで等圧線表示を試みる。

import matplotlib.pyplot as plt
from scipy.interpolate import griddata

(中略)

# 時刻ごとの気圧データ抽出
JST17 = df['JST17']

# 等圧線表示
x, y = np.meshgrid(lon, lat)
z = griddata((lon, lat), JST17, (x, y), method='linear')
levels = [1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014] # 表示する等圧線の設定
cs1 = ax.contour(x, y, z, transform=ccrs.PlateCarree(), 
               levels=levels, linewidths=1.2, colors='k', zorder=9)
cs1.clabel(fmt="%.0f") # 等圧線のラベル表示

plt.show()

なんだこれは！(岡本太郎)

タタリ神のなりそこないみたいなものが出力されてしまった……

原因を探ってみたが、補間する座標(x, y)の設定が悪さをしていたようで次のように直してやればそれなりな図が出てきた。

# 等圧線表示_一応成功版
x, y = np.meshgrid(np.unique(lon), np.unique(lat)) # lon, latをnp.uniqueで抽出。　
z = griddata((lon, lat), JST17, (x, y), method='linear')
levels = [1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014]
cs1 = ax.contour(x, y, z, transform=ccrs.PlateCarree(), 
               levels=levels, linewidths=1.2, colors='k', zorder=9)
cs1.clabel(fmt="%.0f")

method='linear'

かなり直線的な等圧線ではあるが、概ね理想通りの等圧線である。補間をlinearからcubicに変えてみる。

method='cubic'

cubic補間だと線は滑らかになるぶん3次関数的に補間されてしまうためか、松本東方の領域に1007hPa以下の領域が現れている。なるほど、Prologueに示した天気図が納得いかなかった理由が読めてきた。あとはここからどう改良するかだが、その話はまた後日。

ということで今回はここまで。第4話は首を長くしてお待ちを。
(そろそろこのシリーズ以外のことも書きたいななどと思いながら)