数十年前に山登りの講習を受けたとき，自分で短時間に天気図を描けるように練習しました。天気図の情報はNHKラジオ第2で決まった時間に放送されており，天気図を描くための地図が登山用品店で売られていました。ラジオからは，アジア地区の各観測所の天気・気圧・風向・風力や高気圧・低気圧の位置情報などが提供されていました（いまでも，放送はしているみたいです）。その放送から「木浦（モッポ）」「厦門（アモイ）」など，あまりなじみのない地名と場所を覚えたものです。放送のあと，等圧線や前線を描いて仕上げるのですが，きれいに仕上げるにはセンスが必要でした。
しかし，いまではコンピュータが解析したデータに一般人でもアクセスでき，研究用であれば無料で使用することができます。そこで，プログラム言語であるpythonを使って天気図を作成できるようにしました。なお，私はプログラミングの初心者です。そのため，パソコンでpythonが動くようにしたり，いろいろなやり方をネットやchat-GPTで調べながら四苦八苦して取り組みました。
Pythonがパソコンで動くようにする環境を整える
いちばん参考になったのは下記リンクのサイトでした。
参考）Windowsの場合｜Pythonでやってみる気象データの可視化と解析
これに従い，「Miniconda」をインストールし，レポジトリを初期状態から「conda-forge」に変更しました。指示に従い仮想環境も構築しました。
Pythonのコードを書いたり，実行・デバッグをするツールとして「PyCharm」をインストールしました。PyCharmの使い方はYouTubeにより独学で勉強しました。最初は慣れない中で始めたのでなかなかの苦痛でしたが，最低限の使い方を覚えてプログラミングを開始しました。
Pythonのプログラミングについて，私はほぼ初心者です。大昔にFortranを少し勉強したことがある程度です。Pythonの本を購入したり，ネットで検索したりして変数の扱いや基本的な構文を理解しました。当然，可読性やメンテナンス性に優れるオブジェクト指向のプログラミングができるわけもなく，寄せ集めの力業で書いたスパゲティコードで進めています。
天気図作成のステップ
Pythonが使えるようになったら，天気図を作成するために何が必要か考えます。大きく分けて以下のステップがあります。
- 作画用データの入手
気象庁がコンピュータで計算しているデータを研究用に公開しているサイトがあるので，そこにアクセスしてダウンロードする - 入手したデータを使えるようにする
データに何が含まれるか理解し，必要なデータを取り出す - 天気図を描くための地図を用意する
- 天気図を描画する
今回は地上天気図だけでなく高層天気図も作成するので，天気図に何の物理量をどのように描くか決めて作図する - 天気図の保存
作図用データの入手
気象庁が出しているデータは商用で使う場合は有料ですが，研究用に無料で公開しています。例えば下記サイトから入手できます。
参考）https://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/
予報の気象データをダウンロードする場合，各時間のファイルを読み込む必要があるので手動でダウンロードするのは大変です。Pythonでダウンロード用のプログラムを作成しました。参考にしたサイトは以下のものです。
参考）【Python】Pythonでネット上にあるファイルをダウンロードする方法 #Python – Qiita
今回は気象庁のGSMのデータをダウンロードしました。データは計算の初期値となる日時のものと，そこから先の6時間ごとのデータになります。ファイル名には命名規則があるので，目的とするファイル名をプログラム内で指定し，順番に自動的にダウンロードできるようにしました。
GSMデータから必要な物理量の取り出し
GSMのデータの形式は「grib2」で，Pythonでこのファイルから必要なデータを取り出すものとして「pygrib」というライブラリがあります。また，気象データの処理や計算を行うライブラリ「MetPy」も使います。
処理の流れとしては
- GSMのデータをpygribを使ってPythonで開く
- 何の物理量が欲しいか指定して，pygribを使ってデータを取り出す
- 取り出したデータは物理量，緯度，経度のndarrayが入ったタプルになっている
- 取り出したデータをアンパックして物理量，緯度，経度の変数に入れる
- 物理量には必要に応じて単位を付ける（MetPyを使用）
- 渦度・相当温位・湿数など，GSMデータにないものは計算する（MetPyを使用）
- ガウシアンフィルタを使って平準化する
- 物理量・緯度・経度のデータから描画する
- 描画した天気図を保存する
参考）GRIB2をPythonで読み込んで気象データを地図に描く方法 – Looop TechBlog
また，2.でpygribでデータを取り出すとき，どの物理量にするか指定する必要がありますが，その方法が3通りあります。”name”，”shortName”，”parametereName”です。一覧を下図に示します。ここを理解するのが大変でした。
高層天気図を描くときの各気圧面には下表のような物理量が提供されています。
天気図上の地図
地図を描くライブラリとして「Cartopy」を使います。下記サイトが詳しいです。
参考）Windowsの場合｜Pythonでやってみる気象データの可視化と解析
PythonでGMTを使わずにGSMのGPVデータから天気図を描く｜天文と気象（とその他いろいろ）
物理量への単位の付与
単位の付与は下図のように行いました。変数の中身を見ると，単位情報が付加されていることがわかります。ここは，バージョンによって書き方がいろいろあるみたいで，ネット検索で理解するのに苦労したところです。
渦度・相当温位・湿数など，GSMデータにないものの計算
ライブラリMetPyが活躍しました。相対渦度の計算については，ネット上で検索してもなかなか理解できる資料がなく，chat-GPTに聞いたらサクッとコードを書いてくれました。
MetPyを使って下の式で計算するのですが，dx，dyの計算式がわからずchat-GPTの提案式（ここには記載せず）を使いました。検証が必要ですので，今後やるつもりです。
import metpy.calc as mpcalc
vort = mpcalc.vorticity(gpv_values_u, gpv_values_v, dx=dx, dy=dy)
湿数については，露点温度が気温と湿度から計算できますので，湿数＝気温-露点温度で求められます。このとき，使用する単位にはお気を付けください。
相当温位は気圧pressure[hPa]，気温gpv_tmp_values[℃]，露点温度[℃]からMetPyを使って次式で求められます。
theta_e=mpcalc.equivalent_potential_temperature(pressure,gpv_tmp_values.to(units('degC')),dewpoint)
まとめ
最初に掲載したのは海面更正気圧と850hPaにおける気温でした。最後に，500hPaにおける高度と気温を示します。気圧面と描画したい物理量を指定すれば，いろいろな天気図を描くことができます。今後は，より見やすい天気図を目指して改良していきます。また，いまはスパゲティコードですが，オブジェクト指向のプログラムにして見やすく改良しやすいコードにしていきたいと思います。皆さんにお知らせできるようなトピックができたらまた記事を掲載するつもりです。
