気象庁のMSM解析値 (京大生存研のNetCDF版) を Python xarray で読む

気象データとしては、観測値のほかに、数値予報の初期値のためにつくられた格子点データがよくつかわれる。初期値をつくるためには、観測値と(前の時刻の初期値からの)予報値とをくみあわせる「データ同化」とよばれる技術がつかわれる。観測値をとりこんで得られた格子点値を「解析値」あるいは「客観解析値」という。それをそのまま初期値にするばあいと、初期値にするまでにもう一段の処理があるので解析値と初期値とが区別されるばあいがある。

格子点データのファイル形式として、気象庁では、WMOできめられた GRIB version 2 (GRIB2) という形式をつかっているが、研究者は、アメリカの Unidata できめられた NetCDF という形式をつかっていることが多い。ここでは、気象庁の MSM (メソスケールモデル) の解析値を、京都大学生存圏研究所で NetCDFに変換された形で、読むことにする。

[別のページ] では GrADS で読んだが、ここでは、Python の xarray パッケージで読むことにする。

京大生存研にある MSM データ

京大生存研の「生存圏データベース」http://database.rish.kyoto-u.ac.jp のサーバーで公開されている気象庁MSMの格子点データには、気象庁でつくられ、気象業務支援センターから GRIB2 形式で配信されたものと、それを京大生存研で編集して NetCDF 形式にしたものがある。 (なお、全球モデル GSM の格子点データは、GRIB2 形式のものだけが置かれている。)

ここでは、NetCDF形式になったものを読む。「生存圏データベース」→ 「グローバル大気観測データ」→ 「気象庁データ」→ 「解析値を中心に再構成したNetCDFデータ (MSM, RSM)」→ 「解析値を中心に再構成したNetCDFデータ (MSM, RSM)」[OpenDAPでないほう] とリンクをたどると、http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/netcdf/ に行く。

ここに、MSM-P と MSM-S のディレクトリがある。 MSM-P は気圧面 16層、MSM-S は地表面など 1層のデータをさしている。それぞれの下に、「2002」から「2023」までの年 (4けた数字) のディレクトリがある。たとえば「2002」の下をみると「0701.nc」から「1231.nc」までのファイルがある。これは、2002年7月1日から2002年12月31日まで、1日ごとのデータで、そこに、00 UTC から 23 UTC までの、MSM-S では 24時刻、MSM-P では 8時刻の値がふくまれている。

MSM-Sデータにはつぎの変数がふくまれる。 psea (海面気圧)、sp (モデル地表面気圧)、u (風速の西風成分)、v (風速の南風成分)、 temp (気温)、rh (相対湿度)、r1h (1時間降水量)、 ncld (雲量)、ncld_upper (上層雲量)、ncld_mid(中層雲量)、ncld_low (下層雲量)。水平格子間隔は実際は緯度経度だが約 5 km である。時間間隔は 1 時間だが、そのうち hour が 3の倍数の時刻は「初期値」、他の時刻は予報値がはいっている。降水量はすべて予報値である。観測値の情報がとりこまれるのは初期値がつくられるときである。

MSP-Pデータにはつぎの変数がふくまれる。 z (ジオポテンシャル高度)、w (鉛直p速度)、u (風速の西風成分)、v (風速の南風成分)、temp (気温)、rh (相対湿度)。【ただし、w は、増田が 2023-06-30 までに見た 2022年のいくつかの日時の netCDF ファイルに関するかぎり、すべて 0 がはいっている。】水平格子間隔は実際は緯度経度だが約 10 km (MSM-Sの格子のひとつおき)である。時間間隔は 3 時間で、(2006年3月以降は) すべて初期値である。

(初期値・予報値の情報は http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/netcdf/README による。)

モデルの地表面の標高の情報

MSMモデルでつかっている地表面の標高 (ジオポテンシャル高度) を格子点データの格子に内挿したものは、 http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/etc/ にある TOPO.MSM_5K である。 (同じディレクトリにある「README（支援センター）_201406.txt」、「README（支援センター）_201702.txt」の情報によった。) TOPO.MSM_5K のタイムスタンプは 2017年1月30日で、その他の時期も同じ値だったかどうかは未確認である。 (変化があっても数値としては小さなものらしいが)。このファイルは netCDF でも GRIB でもなく単純バイナリ (4バイト浮動小数点数、big endian) である。

このデータを、 [地形 (標高) データを Pythonで読み、作図する (試作版)] と同様に Python で読み、そのうち (下にのべる) 「関東」領域の値を地図上に表示するプログラムを作成した。【このプログラムは、Windows上の Anaconda Power Shell で実行し、 savefig() で png 画像ファイルをつくった。Spyder での動作は未確認。】【このプログラムは、contouf が cmap.set_under, cmap.set_over を認識して処理することを前提としている。もし set_under, set_over が無効ならばリスト levels の下端に下限の値、上端に上限の値、 colors の下端と上端に色を追加して対処すればよい。】

[TOPO.MSM_5K] 地形 (標高) データ [京大生存研から 2021-11-02 にダウンロードしたもの]
【[リンク訂正 2023-02-10] 】[map_msm_topo_kanto.py] 「関東」領域の地形標高段彩をするプログラム [結果画像例]

Python で NetCDF を読むためのパッケージ

Python で NetCDFを読み書きするためのパッケージとしては、netcdf4 というものがある。 (この 4 は NetCDF の version 4 をさしている。) そのほかに、xarray というものがつくられ、これでも NetCDF を読むことができる。 xarray のほうが、読んだ結果が numpy の ndarray と互換性が高いようなので、ここでは xarray をつかうことにする。

京大生存研で編集された気象庁 MSM の NetCDF を読んで作図するプログラム例として、京大防災研の榎本剛 [たけし] 先生によるものがある。これは MSM-P を読むものである。

(外) (京大防災研) [Python による気象データサイエンス: 数値予報データ]

ここでは、このプログラム例をもとに、MSM-P と MSM-S を読んで地図上に作図するプログラム例をつくってみた。

ローカルマシンでのディレクトリ編成

ディレクトリ編成をつぎのようにしておき、それにあわせてデータを配置する。【このデータセットのファイル名は、MSM-S と MSM-P の区別や年の区別がつかないので、とりちがえないように注意する。 NetCDFなので、ファイルの中に、変数名、格子数、対象日時などの情報はあるが。】

netcdf_p/
- 2022/
  - 1221.nc
  - ...
  - 1231.nc
- 2023/
  - 0101.nc
  - ...
  - 0105.nc
netcdf_s/
- 2022/ (上と同様)
- 2023/ (上と同様)
topo/ ... MSMの地形標高データ

hamana (Linuxマシン) のばあい (2023-06-30)

/raid/jma_msm/ の下に、上記の netcdf_p, netcdf_s, topo ディレクトリがある。 netcdf_p, netcdf_s には、2011-01-01 から 2023-06-28 までのデータを置いた (ときどき更新する予定)。

/home/masudako/ex/msm の下に python プログラム例がある。「cp -p /home/masudako/ex/msm .」として各自の作業ディレクトリをつくり、「cd msm」「python3 map_msm_p_wind_kanto_topo.py」のように実行することを想定している。

Windowsパソコンのばあい

2023年1月の実習用に、2022年12月21日から2023年1月5日までのデータをダウンロードしてみた。合計のデータ量は約 5ギガバイトである。 (1日ぶんの MSM-P が約 180 メガバイト、MSM-S が約 135 メガバイトである。) このうち12月21日から12月25日までの5日ぶん (データ量約 1.5 ギガバイト) を USB接続の外づけハードディスクで回覧し、各自コピーしてもらうことにした。

そのときの外づけハードディスクのルートディレクトリのトップの下に jma_msm があり、その下に上記の netcdf_p, netcdf_s, topo があり、またそれと並列に py があってそこにプログラム例を置いた。

気圧面や地表面での水平分布を地図上に作図するプログラム例

いずれも、日本語文字はコメントだけにふくまれ、その文字コードは UTF-8 になっている。

日付は「ds = 」の文のファイル名できまる。

ここでのプログラム例のならべかたは、プログラムをつけたしていったなりゆきで、

まず、1日8回の図を、横に4つ、縦に2行の形でならべて、1枚の画像ファイルにするもの
そのあとに、1時刻ごとに別々の画像ファイルにするもの

となっている。しかし、これからはじめる人への参考としては、下の 1ページに1時刻の画像だけを表示するプログラム例をさきに見たほうがよいと思う。

「関東」領域は、北緯 34.5～38.5度、東経 136～141度。実際には中部地方もふくむ。

MSM-P のひとつの気圧面 (例、850 hPa) の気温
- [map_msm_p_temp_japan.py] MSM 全領域、段彩 [結果画像例] (-30℃ 未満のところが白ぬきになっている。)
- [map_msm_p_temp_kanto.py] 「関東」領域、段彩と等値線をかさねる [結果画像例]
MSM-P のひとつの気圧面 (例、850 hPa) の風速
- [map_msm_p_wind_kanto.py] 「関東」領域、風速ベクトルの矢印と、風速の大きさの色分けと等値線をかさねる [結果画像例]
MSM-S の海面気圧
- [map_msm_s_psea_japan.py] MSM 全領域、段彩と等値線をかさねる [結果画像例]
- [map_msm_s_psea_kanto.py] 「関東」領域、段彩と等値線をかさねる [結果画像例]
MSM-S の地上 (2 m) 気温
- [map_msm_s_temp_japan.py] MSM 全領域、段彩 [結果画像例]
- [map_msm_s_temp_kanto.py] 「関東」領域、段彩と等値線をかさねる [結果画像例]
MSM-S の地上 (10 m) 風速。水平格子がこまかいので、skip を指定して矢印をかく格子点をまびいている。
- [map_msm_s_wind_kanto.py] 「関東」領域、風速ベクトルの矢印と、風速の大きさの色分けと等値線をかさねる [結果画像例]
ここまでのプログラム例の動作確認は Windows 上の Anaconda の Spyder でおこない、 plt.show() を実行して、右上の窓の plot で見られる画像を save した。
- [2023-06-08 補足] この一連の例文プログラムは、Spyder で実行することを想定してつくっている。シェルプロンプトで実行すると、図の大きさと文字の大きさのバランスがわるく、フォントサイズなどを変更する必要があるが、どう変更したらよいかは確認できていない。
- ファイル名は save した日時 (対象データの日時ではない) がはいったものになるが、あとに手作業でつけなおした。日時の数がふえると、これではめんどうだし、まちがえやすい。 plt.savefig() で画像ファイルをつくり、そのファイル名を日付によってきめることもできるが、まだプログラムにくみこんでいない。【たとえば [2023-01-10の教材ページ] のプログラム例 map_windarrow.py の最後のコメントにしてある部分を参考にしてほしい。】

未解決の問題点

風速と矢印の長さの関係は、変数 arrowscale であたえているが、その定義がよくわからず、試行錯誤できめている。【別の教材ページの例文では格子点ごとに quiver をよびだしたが、ここでは読みこんだ全領域の風ベクトルの u, v の配列を quiver にわたしている。】
quiverkey でかかれる凡例の矢印が見えにくい。
quiverkey でかかれる凡例の矢印を書く位置が、段彩の colorbar をかくときに確保された場所であり、等値線図をかかないと plt.show でしめされる図の外になってしまう。
contourf の階級しきい値と contour の線をひく数値は同じにしたが、いずれも仮にあたえたものであり、データ内容に適切か、検討できていない。 contourf のカラーマップは coolwarm または Orange をつかった。 contour の線の色は指定せず、coutour プログラムの機能にまかせた。これもくふうの余地がある。

1ページに1時刻の画像だけを表示するプログラム例 [2023-06-05 追加]

こちらは Spyder の外でもだいたい適切なフォントサイズになる。しかし plt.show() をつかう形のままである。

MSM-P のひとつの気圧面 (例、850 hPa) の風速
- [map_msm_p_wind_kanto_single.py] [結果画像例] 「関東」領域、風速ベクトルの矢印と、風速の大きさの色分けと等値線をかさねる
  - 「for hr in [0, 3, 6, 9, 12, 15, 18, 21]:」の角かっこ中のリスト内容を、作図したい時刻 (UTC の hour) だけにするとよい。
  - [2023-06-08 更新] 緯度・経度の数字ラベルをつけるようにした。
MSM-P のひとつの気圧面 (例、850 hPa) の風速と比湿 (specific humidity) [2023-08-23 追加]
- [map_msm_p_wind_sph_kanto_single.py] 「関東」領域、風速ベクトルの矢印と、比湿の色分けと等値線をかさねる。比湿は、相対湿度と気温から、気温から飽和水蒸気圧をもとめる Tetens の式 ([(気候・気象学) 水の相変化、大気中の水蒸気] 参照) をつかって計算している。
MSM-S の地上 (10 m) の風速 [2023-06-30 追加]
- [map_msm_s_wind_kanto_single.py] [結果画像例]
- MSM-P のばあいと同様。
- ただし、MSM-S データは格子間隔がこまかいので、矢印を表示する格子点を、東西、南北それぞれひとつおきにした。

風と地形標高をかさねる (1ページに1時刻の画像だけを表示する) [2023-06-18 追加]

MSM-P のひとつの気圧面 (例、850 hPa) の風速ベクトルと、MSMの地形の標高段彩をかさねる。風速の大きさの表示はしない。
- [map_msm_p_wind_kanto_topo.py] [結果画像例]

南北・鉛直断面を作図するプログラム例 [2023-02-10 追加]

MSM-P の 3次元格子データで、横に南北、縦に気圧をとった 2次元の等値線図をかき、色でぬりわける。

MSM-P の風速の東西成分 (西風が正) の、ひとつの経度での南北・鉛直分布 (横軸は緯度、縦軸は気圧 p [hPa])
- [u_lat_p.py] [結果画像例 (2022-12-23 09 UTC)]

未解決の問題点

縦軸に気圧をとるときは、対数尺度にしたほうが、実際の高さと比例するものに近くなる。プログラム例にコメントとして入れたように「ax.set_yscale('log')」をつかえば、対数尺度になる。しかし、そうすると、縦軸の目盛りにつける数値の「200」のようにしたいところが「2×10²」のような指数表示になってしまい読みにくい。その対策ができたら、対数尺度のプログラム例を追加したい。

応用

lon と lat のたちばをいれかえれば、ひとつの緯度での東西・鉛直断面を作図することができる。

リンク集

数値予報の初期値のつくりかたについてのおおまかな話 (増田の教材ページ 2008年) [気象の数値予報でのモデルの使われかた] の「初期値の作りかた ... データ同化」のところ
気象庁 (www.jma.go.jp)
京大生存研 (www.rish.kyoto-u.ac.jp)
- [気象庁データ]
- [NetCDF化した数値予報GPVデータ]
  - そのページにある「解析値を中心に再構成したNetCDFデータ (MSM, RSM)」の説明から引用【増田注: 現在のファイル編成では、降水量・日射量は地表面データに統合されたらしい。】
    - MSM, RSMのオリジナルデータを初期値(＝解析値)を中心に再構成した，1日1ファイルのNetCDFデータです．
      - 1日分のオリジナルデータの初期値を連ね，その間を直前の初期値からの予報値で埋める構成としています．このように再構成することで，過去の大気状態のベストな推定アーカイブとなっています．
      - 1時間降水量・日射量は初期値がありませんので，予報値のみを連ねた1時間データを別ファイルで作成しています．
    - データに含まれている物理量は以下のとおりです:
      - 気圧面データ(MSM-P, RSM-P) … ジオポテンシャル高度，東西風，南北風，気温，鉛直流(RSMと2006年2月28日までのMSMは300hPaまで)，相対湿度(300hPaまで)
      - 地表面データ(MSM-S, RSM-S) … 海面更正気圧，地上気圧(2006年3月1日以降のMSMのみ)，東西風(地上10m)，南北風(地上10m)，気温(地上1.5m)，相対湿度(地上1.5m)，上層雲量，中層雲量，下層雲量，全雲量(2006年3月1日以降のMSMのみ)
      - 1時間データ(MSM-S/r1h, RSM-S/r1h) … 降水量，日射量
    - 各データのディレクトリの下に yyyy/mmdd.nc (yyyy: 年, mmdd: 月日; UTC) のファイルがあります.
    - 各ファイルの具体的な時間構成については，READMEを参照してください．
- [解析値を中心に再構成したNetCDFデータ (MSM, RSM)] (データをダウンロードするための入口)
GRIB2について
- [WMO（世界気象機関）の格子データ形式GRIB2について] (@e_toyoda さん、2019年07月11日に投稿、2020年03月16日に更新)
netCDFについて
- [NetCDF] 本家 Unidata (英語)
- [[NetCDF]] Wikipedia日本語版
- [NetCDF に関する情報] 地球流体電脳倶楽部

2023-01-16, 2023-02-10, 2023-06-18, 2023-06-30, 2023-08-22
増田耕一 (MASUDA Kooiti)

増田の個人サイトでの [計算機を利用した気象分野の教材のページ] [立正大学の授業用の教材ページ] [教材目録ページ] [日本語入り口ページ] [ホームページ]へ