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

[2025-11-14 おおはば改訂しました。]

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

格子点データのファイル形式として、気象庁では、WMOできめられた GRIB version 2 (GRIB2) という形式をつかっているが、研究者は、アメリカの Unidata できめられた NetCDF という形式をつかっていることが多い。

ここでは、気象庁の MSM (メソスケールモデル) の解析値の、京都大学生存圏研究所で NetCDFに変換された形のものを 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」から「2025」までの年 (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 (下層雲量)。水平格子間隔は 2006年3月以後は約 5 km (正確には、緯度 0.05 度、経度 0.0625度)、 2006年2月まではその2倍のあらさの約 10 km である。時間間隔は 1 時間だが、そのうち hour が 3の倍数の時刻は「初期値」、他の時刻は予報値がはいっている。降水量はすべて予報値である。観測値の情報がとりこまれるのは初期値がつくられるときである。

MSP-P データにはつぎの変数がふくまれる。 z (ジオポテンシャル高度)、w (鉛直p速度、気象学の教科書では ω [オメガ] と書かれることが多い)、u (風速の西風成分)、v (風速の南風成分)、temp (気温)、rh (相対湿度)。水平格子間隔は、2006年3月以後は約 10 km、2006年2月までは約 20 km (いずれも MSM-S の格子のひとつおき)である。気圧レベルは16個あり、1000, 975, 950, 925, 900, 850, 800, 700, 600, 500, 400, 300, 250, 200, 150, 100 hPa である。時間間隔は 3 時間であり、(2006年3月以降は) どの時刻の値も初期値である。

[注] ただし、w は、京大生存研のウェブページの注記 (2018-11-14) によれば、 2017-02-28 03 UTC 以後、すべて 0 がはいっている。鉛直流の初期値を観測値にもとづいてつくることはむずかしいので、初期値は 0 として、予報が進行すると値が変化することにしたのだ。
[注] 初期値・予報値の情報は http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/netcdf/README による。

ここにある MSM解析値の格子は、 3次元データの鉛直レベルは p 面であり、メソスケール数値予報モデルのレベルと同じではない。水平位置は緯度・経度それぞれ等間隔の格子であり、メソスケール数値予報モデルの格子とは同じではないらしい。したがって、配信された解析値がつくられるまでに、モデル格子からの空間内挿がおこなわれている。モデルの空間分解能は、2006年3月以後、ほぼ一定にたもたれているようだが、気象庁のモデルが改訂されると、その前後でデータの質がちがうかもしれないことに、注意が必要だろう。

研究室 Linux マシンでのデータの配置

hamana と niigata には、 2025-11-13 現在、2002-07-01 から 2025-11-11 までの日付のデータがある。 /raid/jma_msm/ の下に、つぎのようなディレクトリ編成でデータを配置している。

/raid/jma_msm/
- netcdf_p/ ... MSM-P 解析値データ
  - 2002/
    - 0701.nc
    - ...
    - 1231.nc
  - 2003/
    - 0101.nc
    - ...
    - 1231.nc
  - ...
  - 2025/
    - 0101.nc
    - ...
    - 1111.nc
- netcdf_s/ ... MSM-S 解析値データ
  - (上と同様)
- topo/ ... MSMの地形標高データ

補足説明

データ量は、(2006年3月以後の値のばあい) 1日ぶんの MSM-P が約 180 メガバイト、MSM-S が約 135 メガバイトである。とくに自分のパソコンで作業したいのでなければ、niigata または hamana で作業することをすすめる。
ファイル名ではMSM-S と MSM-P の区別や年の区別がつかないので、ダウンロードしてつかうばあいは、とりちがえないように注意してほしい。 (NetCDFなので、ファイルの中には変数名・格子の座標値・年月日などの情報はあるが。)
niigata での /raid/jma_msm は、/diskb/jma_msm に symbolic link でつけた別名である。

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

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

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

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

プログラム例

ここでは、榎本先生のプログラム例を参考にして、 xarray をつかって MSM-P と MSM-S を読んで、 cartopy で地図上に作図するプログラム例をつくってみた。

【2006年2月以前のデータには、それ以後のデータと仕様のちがいがある。増田がつかんでいるのは、ひとつは、格子間隔が、以前のデータは以後のデータの2倍であること、もうひとつは、気温の単位が、以前は ℃ (セルシウス度)、以後は K (ケルビン) であること、の2点である。つぎのブログ記事にややくわしく書いた。

[MSM 解析値 2006年2月以前は仕様がちがう] (2025-09-17) 】

2025年11月のプログラム例改訂で、この2つの問題点には対処し、プログラムを前後の期間を通じてつかえるようにした。しかし、2006年2月以前のデータについての動作テストは不じゅうぶんであり、未解決の問題があるかもしれない。】

上記のほかにもいくつかプログラム例をつくっているので、このページ、またはこのページからリンクされたページに追加していく予定。

【複数のプログラムで同じことをやっているところは、ひとつのファイルに書いて各プログラムからそれを共有するようにしたいのだが、それを実現する方針はまだたっていない。】

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

2025-11-13 現在、niigata と hamana の /home/masudako/ex/msm にある。動作確認は niigata でおこなった。プログラムのファイル名には「_」 (下線、アンダースコア) をふくんでおり、空白文字はふくんでいない。

MSM-P のひとつの気圧面 (例、850 hPa) の風ベクトルの矢印表示と、風速 (ベクトルの大きさ) の値によるぬりわけ
- [map_msm_p_wind_v20251103.py]
MSM-P のひとつの気圧面 (例、850 hPa) の風ベクトルの矢印表示と、気温によるぬりわけ
- [map_msm_p_wind_tmp_v20251103.py]
MSM-P のひとつの気圧面 (例、850 hPa) の風ベクトルの矢印表示と、相対湿度によるぬりわけ
- [map_msm_p_wind_rh_v20251103.py]
MSM-S の地上 (高さ 10 m) の風ベクトルの矢印表示と、風速 (ベクトルの大きさ) の値によるぬりわけ
- [map_msm_s_wind_v20251103.py]
MSM-S の地上 (高さ 10 m) の風ベクトルの矢印表示と、地上 (高さ 2 m) 気温によるぬりわけ
- [map_msm_s_wind_tmp_v20251103.py]

つぎの注は、ここにあげた例文に共通である。

プログラム (スクリプト) ファイル中では、日本語文字はコメントだけにふくまれ、その文字コードは UTF-8 になっている。ブラウザでは文字化けするかもしれないが、ブラウザに「テキストエンコーディングを修復」の機能があればそれをつかって読める。なければ、(右クリックから) ダウンロードしてパソコンで読むか、 niigata や hamana 上の emacs で読んでほしい。
気圧面のデータの作図では、気圧の値を hPa 単位で plev に代入する値としてあたえる。
対象日時を iyear, imon, idaystart, idayend に代入する値としてあたえる。 (1回の実行で処理できる期間はひとつの年のひとつの月のうちにかぎっている。) すると、idaystart から idayend までの毎日の世界時 (UTC) 0, 3, 6, 9, 12, 15, 18, 21時の図がつくられる。時刻をさらに限定したいときは、「for ihour in [...]」の角かっこ内を書きかえればよい。
緯度経度範囲は、lonmin, lonmax, latmin, latmax にあたえる値で指定する。いくつかの緯度経度範囲に略称をつけ、画像ファイル名にもそれをふくめることにした。仮の代表例として、北緯 34.5～38.5度、東経136.5～141.5度の領域 (関東地方のほか、中部地方の大部分と東北地方の南部をふくむ) を、 (うまい表現がなくてくるしまぎれだが) 「gkanto」 ("Greater Kanto", 「関東広域」) とした。
- 略称がつくられていない緯度経度範囲で作図したいときは、プログラム中で最後に指定している緯度経度範囲の値を書きかえたうえで、「region = ""」として変数 region に空文字列をあたえるのがよい。
- 【なお、「zenkoku」という緯度経度範囲略称は、解析雨量のデータが存在する領域をふくむように設定したものであり、 MSM解析値のデータが存在する領域よりもせまい。】
図は、変数 savefig の値が True ならば、imagefilepath で指定したところ (例のまま変更しなければ、"../image/") に png 画像ファイルで保存される。 False ならば X Window の画面表示となる。
画像ファイル名には、変数 savejst の値が True ならば日本標準時 (JST)、 False ならば世界時 (UTC) で日時がはいるようにしてある。
気温についての気象要素略号は NetCDF ファイル内では「temp」になっているが、増田のプログラム例や出力画像ファイル名では (増田の他の作業にあわせて) 「tmp」とした。
全部の格子点に矢印をかくとこみあってしまうので、東西方向、南北方向それぞれいくつおきの格子点にかくかを skip という変数で指定している。 skip の値を、MSM-P と MSM-S で、また、2006年2月までと2006年3月以後とで変えて、地球上での間隔がほぼ同じになるようにした。
風速と矢印の長さとの関係は、scale にあたえる数値できめている。同じ風速ならば、数値が小さいほど矢印が長くなる。凡例の矢印があらわす風速を、qk にあたえる数値で指定している。
風速の大きさや気温の数値の階級わけは配列 vallevs にあたえる数値で、色づけは、cmap にあたえるカラーマップ名 (文字列、matplotlib の Colormap reference 参照) で指定している。 (正と負の値をとる量のばあいは "coolwarm" で、負を青、正を赤にする。符号が一定の量のばあいは、"Blues" で青の濃淡表示にしてみた。)
- 各階級の色を直接にきめる方法の例は、 [地図をかく (3) 地形の標高を色で表示する] などを参照。
- [外] (matplotlib) [Colormap reference]
海岸線のほか、都道府県の境をいれている。その方法は [日本の都道府県境をかく / 中国の「省」の境をかく] に書いた。niigata・hamana 以外の計算機で実行するときは、あらかじめ都道府県境データを用意する必要がある。

等圧面・地表面のデータから気象変数を計算して地図上に作図するプログラム例

[注] 鉛直積算をするものは計算時間がかかるので、[プログラムをバックグラウンドで実行する] 方法をとるのがよいだろう。

水蒸気輸送に関するもの

MSM-P のひとつの気圧面 (例、850 hPa) の水蒸気輸送ベクトルの矢印表示と、比湿 q によるぬりわけ。比湿は空気の質量のうちの水蒸気の質量のわりあいである。気温から飽和水蒸気圧を Tetens の式 ([気候・気象学教材水の相変化、大気中の水蒸気] のページ参照) によって計算し、相対湿度とくみあわせて比湿を得る。ここでいう水蒸気輸送ベクトルは、比湿と風速をかけたものである。
- [map_msm_p_quv_v20251103.py]
MSM-P のひとつの気圧面 (例、850 hPa) の水蒸気輸送ベクトルの矢印表示と、水蒸気輸送ベクトルの絶対値によるぬりわけ。
- [map_msm_p_quvmag_v20251103.py]
MSM-P のひとつの気圧面 (例、850 hPa) の水蒸気輸送ベクトルの矢印表示と、水蒸気輸送ベクトルの水平発散によるぬりわけ。水平発散は、東西、南北それぞれ両どなりの格子点の値の差をとる「中央差分」で計算している。作図領域の端にあたる格子では水平発散を計算していない。
- [map_msm_p_divquv_v20251103.py]
地表面から MSM-P のひとつの気圧面 (例、500 hPa) まで鉛直積算した水蒸気輸送ベクトル (面積あたりの質量と風速の積) の矢印表示と、水蒸気量 (面積あたりの質量) によるぬりわけ。 MSM-P と MSM-S の両方をつかって計算する。
- [map_msm_quvcol_v20251103.py]
地表面から MSM-P のひとつの気圧面 (例、500 hPa) まで鉛直積算した水蒸気輸送ベクトル (面積あたりの質量と風速の積) の矢印表示と、その発散によるぬりわけ。水蒸気輸送ベクトルの発散 (値が負のばあいに「収束」) は、液体の水の標準の密度をつかって換算すると、降水量と同じ物理量の次元になるので、1時間降水量と同じ単位に換算してしめしている。
- [map_msm_divquvcol_v20251106.py]

風の発散と鉛直流に関するもの

地表面から MSM-P のひとつの気圧面 (例、500 hPa) まで風速を鉛直積算することによって求めた空気の質量輸送ベクトル (面積あたりの質量と風速の積) の発散によるぬりわけ。このベクトルの発散は、その気圧面の鉛直 p 速度 ω [Pa/s] の推定値になっている。【このプログラムは上記の divquvcol のプログラムを改変して作成作業中のものであり、変数名やコメントには不適切なところがある。】
- [map_msm_omega_v20251106.py]

モデルの地表面の標高とそれを作図するプログラム例

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) である。

hamana および niigata では、このデータを /raid/jma_msm/topo/TOPO.MSM_5K として置いた。

このデータを、Python で読み、地図上に標高段彩図を表示するプログラムを作成した。

2025-11-13 現在、つぎのものが niigata と hamana の /home/masudako/ex/msm にある。動作確認は niigata でおこなった。プログラムのファイル名には「_」 (下線、アンダースコア) をふくんでおり、空白文字はふくんでいない。

[map_msm_topo_v20250503.py]

緯度経度範囲は、例として、上記の「関東広域」としてある。

標高の値を階級わけして色と対応づける方法は、上記の風速や気温のばあいとちがい、 levels と colors の配列に値をあたえることによっている。 [気象庁の解析雨量を Pythonで読み、作図する] での降水量の階級わけと同じ方法である。

南北・鉛直断面や、東西・鉛直断面を作図するプログラム例

プログラム例はいずれも、年、月、日を指定し、世界時 (UTC) でその日のうち 3時間おきの 8時刻 (リスト変数 hours の内容を変更することによってしぼりこむことができる) について作図する。

鉛直の範囲は気圧 (hPa) で pmin, pmax の値で指定する。 pmin が上にくるように作図する。

等値線および色塗りの境界はリスト変数 vallevs の値で指定している。色はカラーマップで、正と負がありうる変数については coolwarm、符号が一定の変数については Blues を指定している。

savefig が True のときPNG画像ファイル出力、False のとき画面表示となる。

画像ファイルのファイル名には、savejst が True のとき日本標準時 (JST) で、 False のとき世界時 (UTC) で日時がはいる。

指定した経度での南北・鉛直断面

経度 (東経が正) を lonsel で、緯度 (北緯が正) の範囲を latmin, latmax で指定する。

風速の東西成分 (u) の南北・鉛直断面
- msm_u_lat_p_v20251114.py
風速の南北成分 (v) の南北・鉛直断面
- msm_v_lat_p_v20251114.py
気温 (tmp) の南北・鉛直断面
- msm_tmp_lat_p_v20251114.py
相対湿度 (rh) の南北・鉛直断面
- msm_rh_lat_p_v20251114.py
比湿 (q、気温と相対湿度から計算) の南北・鉛直断面
- msm_q_lat_p_v20251114.py
温位 (theta、気温と気圧から計算) の南北・鉛直断面
- msm_theta_lat_p_v20251114.py

指定した緯度での東西・鉛直断面

緯度 (北緯が正) を latsel で、経度 (東経が正) の範囲を lonmin, lonmax で指定する。

風速の南北成分 (v) の東西・鉛直断面
- msm_v_lon_p_v20251114.py

リンク集

数値予報の初期値のつくりかたについてのおおまかな話 (増田の教材ページ 2008年) [気象の数値予報でのモデルの使われかた] の「初期値の作りかた ... データ同化」のところ
気象庁 (www.jma.go.jp)
京大生存研 (www.rish.kyoto-u.ac.jp) の生存圏データベース (database.rish.kyoto-u.ac.jp)
- グローバル大気観測データ > 気象庁データ > [解析値を中心に再構成したNetCDFデータ (MSM, RSM)]
  - ページの表題は「気象庁データ」「NetCDF化した数値予報GPVデータ」となっている。
  - そのページのうちの「解析値を中心に再構成したNetCDFデータ (MSM, RSM)」という行 (リンクになっている) の (リンクさきではなく) 下につづく本文に、データの説明がある。
  - ただし、そこに書かれた「1時間データ(MSM-S/r1h, RSM-S/r1h) … 降水量，日射量」のファイルは、現在のファイル編成では地表面データに統合されており、r1h ディレクトリは存在しない。
  - データをダウンロードするための入り口は、[解析値を中心に再構成したNetCDFデータ (MSM, RSM)] の行からのリンクさきである。 (「OpenDAP入口」ではない。)
netCDFについて
- Unidata [NetCDF] (英語)
- Wikipedia 日本語版 [[NetCDF]]
- 地球流体電脳倶楽部 [NetCDF に関する情報]
GRIB2について
- Qiita (@e_toyoda さん) [WMO（世界気象機関）の格子データ形式GRIB2について] (2019年07月11日に投稿、2020年03月16日に更新)

(2023-01-16 ... 2023-08-22, 2025-05-03, 2025-05-04, 2025-06-14, 2025-09-18), 2025-11-14, 2025-12-25
増田耕一 (MASUDA Kooiti)

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