WXBC 2020年度テクノロジー研修
「メッシュ気象データ分析チャレンジ！」講習テキスト

1. GRIB2ファイルからの気象データの取り出し¶

目次：¶

  • はじめに
  • 1.1 GRIB2ファイル処理ツールwgrib2
    • 1.1.1 コマンドラインプログラム
    • 1.1.2 コマンドラインプログラムのPythonからの実行
  • 1.2 wgrib2によるGRIB2ファイルの操作
    • 1.2.1 GRIB2ファイル内のデータのインベントリー
    • 1.2.2 GRIB2ファイルからのデータの取り出し
  • 1.3 WXBCオリジナルライブラリ wxbcgribX
    • 1.3.1 Pythonパッケージ Xarray
    • 1.3.2 Datasetオブジェクト
    • 1.3.3 DataArrayオブジェクト
    • 1.3.4 データ読み出し関数 getgpv
  • おわりに
  • ANNEX
  • 著作権について

はじめに¶

　気象庁は、アメダスのような地点ごとのデータとは別に、面的/立体的、時間的な広がりを持つ気象データも配信しています。例えば、気象庁 メソ数値予報モデルGPV (通称MSM-GPV) は、海上を含む日本周辺域の気象予測データで、東西方向に約5km ごと、南北方向に約5km ごと、1 時間間隔で39 時間先までをカバーし、収録される気象要素は、海面更正気圧、地上気圧、風ベクトル（U、V）、気温、相対湿度、降水量、雲量(全雲量、上層雲量、中層雲量、下層雲量)、日射量と多彩です。気象庁は、これを1日に8回も配信しています。

　このほかにも、気象庁は多種多様なメッシュデータを配信していて、ホームページにはそれらの一覧が掲載されています。いくつかのプロダクトについては、ほんの数個ですがサンプルファイルも取得可能です。

https://www.data.jma.go.jp/developer/gpv_sample.html

　 気象庁のメッシュデータは、様々なビジネスで活用の可能性がありますが、残念なことに GRIB2 と呼ばれる特殊な形式で電子ファイル化されており、これからデータを取り出すには少々のスキルが必要とされます。
そこで、第１課では、上述の MSM-GPV を例に、気所庁が配信するGRIB2ファイルからデータを取り出す方法について学びます。

　なお、気象庁は、緯度経度に整然と並ぶ気象データを「GPVデータ」と呼んでいます。

1.1 GRIB2ファイル処理ツールwgrib2¶

1.1.1 コマンドラインプログラム¶

　Pythonには、GRIB2ファイルを操作するモジュール pygrib が存在しますが、残念ながらWindows向けには提供されていません。さらに、提供されているOSであったとしても、気象庁のプロダクトの一部は pygrib で扱うことができません。このため、現状としては、メッシュ気象データをPythonで処理するにしても、GRIB2ファイルの処理の部分に関しては、Python以外のプログラムを用いることが賢明です。
wgrib2 は、アメリカ大気海洋局(NOAA)の気候予測センター(CPC)が開発し公開するGRIB2ファイルの処理ツールで、その性能には定評があり、気象庁のプロダクトも正しく取り扱うことができます。そこで、「メッシュ気象データ分析チャレンジ！」では、GRIB2ファイルの処理に wgrib2 を使用します。

　wgrib2 は、コマンドラインプログラムという種類に分類されるプログラムです。Windows上でこのプログラムファイル wgrib2.exe のアイコンをダブルクリックしても、黒い窓が開いてなにやら高速に表示された後閉じて終わってしまい、利用できません。コマンドラインプログラムは、Windows でコマンドプロンプト（Macではターミナル）と呼ばれるプログラムを実行し、その中で動かすように作られています。そこで、少し寄り道をして、コマンドラインプログラムとコマンドプロンプトについてざっと学習します。

　コマンドプロンプトは、Windowsのスタートボタンから、[Windowsシステムツール]＞[コマンドプロンプト]を選択すると開くことができます（下図）。不愛想なウィンドウで、図中に示される太い下線の部分（これが本当のコマンドプロンプトです）に文字を打ち込むことができます。キーボードから処理を指示する文字列を打ち込み、最後に [Enter] キーを押すことで、指示した処理が実行されます。処理の実行中は文字を入力することはできません。処理が終了すると（本当の）コマンドプロンプトがまた現れます。

image.png

　wgrib2 は、GRIB2ファイルの中に含まれるデータの概要を書き出したり、データが使用する緯度経度座範囲を取り出したり、データをテキストとしてファイルに書き出したりと、様々な機能を有します。どのGRIB2ファイルを対象として、どのような処理をするかは、コマンドプロントの右側に一行の文字列で書き込んで指定します。この文字列のことをコマンドラインと呼びます。コンピューターに処理させる内容をコマンドラインで与えて実行させるので、"コマンドラインププログラム"なのです。
wgrib2で処理を行うの場合、コマンドラインは以下のような文字列になります。なお、｛操作の指示｝にあたる文字列のことをコマンドラインオプションと呼びます。

｛wgrib2の実行ファイル｝_半角空白_｛操作の指示｝_半角空白_｛対象とするGRIB2ファイル｝

　一番簡単な例として、wgrib2に、バージョンを表示さてみましょう。バージョンを表示させるコマンドラインオプションは、「-version」です。また、この処理にGRIB2ファイルは必要ないのでその部分はありません。したがって、コマンドラインは以下の通りです。

c:￥wgrib2￥wgrib2.exe -version

　結果は下図の通りです。バージョンと製作者名が表示されて終了します。

image.png

1.1.2 コマンドラインプログラムのPythonからの実行¶

　 本来はこのように使用するコマンドラインプログラムですが、Pythonのライブラリ subprocess を利用すると、Pythonプログラムの中でコマンドプロンプトを仮想的に作りだし、その中でコマンドラインプログラムを実行することができます。

　ライブラリ subprocess を利用するには、インポート文で利用を宣言する必要があります。以下のCellを実行してください。

In [1]:
# Script 1 #

import subprocess

　ライブラリ subprocess をインポートすると、関数 run が利用可能となるので、これを用いて以下のようにスクリプトを書きます。

rc = subprocess.run("command string",      
                    shell=True, text=True, capture_output=True)
for line in rc.stderr.splitlines():
    print(line)
for line in rc.stdout.splitlines():         
    print(line)

ここで、"command string" は、コマンドラインに打ち込む文字列を表します。

　したがって、Pythonで wgrib2のバージョンと製作者名を表示させるには、上の command string の部分にc:￥wgrib2￥wgrib2.exe -versionを書き込めばよいのです。

　ただし、Windows PC の場合は、文字列に少しだけ修正が必要です。日本語版Windowsにおいてパスを示す文字列は円マーク「￥」なのですが、"command string" に置く際には、Macと同じ半角スラッシュ[ / ]を使用する必要があります。すなわち、以下の文字列を書き込みます。
c:/wgrib2/wgrib2.exe -version
それでは実行してみましょう、下のセルに対し、これを踏まえ、以下のように操作してください。

  1. 下のセルをクリックして枠線を表示させます(青または緑の枠はセルが編集可能な状態にあることを示します)。
  2. "command string" を "c:/wgrib2/wgrib2.exe -version" に書き換えます。
  3. ツールバーの ［>|Run］ ボタンをクリックするか、Ctrlキーを押しながらEnterキーを押します。
In [2]:
# Script 2 #

rc = subprocess.run("command string", 
                    shell=True, text=True, capture_output=True)
for line in rc.stderr.splitlines():
    print(line)
for line in rc.stdout.splitlines():
    print(line)

'command' は、内部コマンドまたは外部コマンド、
操作可能なプログラムまたはバッチ ファイルとして認識されていません。

1.2 wgrib2によるGRIB2ファイルの操作¶

1.2.1 GRIB2ファイル内のデータのインベントリー¶

　wgrib2は、コマンドラインオプションを付けずに対象ファイル名だけを置くと、そのファイルに格納されているグリッドデータの一覧（インベントリー）を表示します。そこで、これを使って、気象庁メソ数値予報モデルGPV（MSM-GPV）ファイルの中にどのような気象データが記録されているかを見ることにしましょう。
2018年7月1日のMSM-GPVデータの中身を見ることにします。コマンドライン

｛wgrib2のファイル｝_半角空白_｛オプション｝_半角空白_｛対象のGRIB2ファイル｝

はどのようにしたらよいでしょうか。

  • {wgrib2のファイル} の部分： この部分の文字列は、置き場所も含めて「c:/wgrib2/wgrib2」です。
  • {オプション} の部分： この部分はありません。
  • {対象のGRIB2ファイル} の部分：
    GRIB2ファイルの置き場所は「jmadata/msm/2018/201807/」で、
    ファイル名は「Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin」なので、
    「jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin」です。

よって、コマンドラインは以下のようになります。
(長いのでウインドウ幅によっては折り返しが入っているかもしれませんが一行の文字列です)。

"c:/wgrib2/wgrib2 jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"

　以下を実行してください。
（Macの方は、「c:/wgrib2/wgrib2」を「~/work/grib2/wgrib2/wgrib2」に書き換えてください）

In [3]:
# Script 3 #

import subprocess

rc = subprocess.run("c:/wgrib2/wgrib2 jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin", 
                    shell=True, text=True, capture_output=True)
for line in rc.stderr.splitlines():
    print(line)
for line in rc.stdout.splitlines():
    print(line)

1.1:0:d=2018070106:PRMSL:mean sea level:anl:
1.2:0:d=2018070106:PRES:surface:anl:
1.3:0:d=2018070106:UGRD:10 m above ground:anl:
1.4:0:d=2018070106:VGRD:10 m above ground:anl:
1.5:0:d=2018070106:TMP:1.5 m above ground:anl:
1.6:0:d=2018070106:RH:1.5 m above ground:anl:
1.7:0:d=2018070106:LCDC:surface:anl:
1.8:0:d=2018070106:MCDC:surface:anl:
1.9:0:d=2018070106:HCDC:surface:anl:
1.10:0:d=2018070106:TCDC:surface:anl:
1.11:0:d=2018070106:PRMSL:mean sea level:1 hour fcst:
1.12:0:d=2018070106:PRES:surface:1 hour fcst:
1.13:0:d=2018070106:UGRD:10 m above ground:1 hour fcst:
1.14:0:d=2018070106:VGRD:10 m above ground:1 hour fcst:
1.15:0:d=2018070106:TMP:1.5 m above ground:1 hour fcst:
1.16:0:d=2018070106:RH:1.5 m above ground:1 hour fcst:
1.17:0:d=2018070106:LCDC:surface:1 hour fcst:
1.18:0:d=2018070106:MCDC:surface:1 hour fcst:
1.19:0:d=2018070106:HCDC:surface:1 hour fcst:
1.20:0:d=2018070106:TCDC:surface:1 hour fcst:
1.21:0:d=2018070106:APCP:surface:0-1 hour acc fcst:
1.22:0:d=2018070106:DSWRF:surface:0-1 hour ave fcst:
1.23:0:d=2018070106:PRMSL:mean sea level:2 hour fcst:
1.24:0:d=2018070106:PRES:surface:2 hour fcst:
1.25:0:d=2018070106:UGRD:10 m above ground:2 hour fcst:
1.26:0:d=2018070106:VGRD:10 m above ground:2 hour fcst:
1.27:0:d=2018070106:TMP:1.5 m above ground:2 hour fcst:
1.28:0:d=2018070106:RH:1.5 m above ground:2 hour fcst:
1.29:0:d=2018070106:LCDC:surface:2 hour fcst:
1.30:0:d=2018070106:MCDC:surface:2 hour fcst:
1.31:0:d=2018070106:HCDC:surface:2 hour fcst:
1.32:0:d=2018070106:TCDC:surface:2 hour fcst:
1.33:0:d=2018070106:APCP:surface:1-2 hour acc fcst:
1.34:0:d=2018070106:DSWRF:surface:1-2 hour ave fcst:
1.35:0:d=2018070106:PRMSL:mean sea level:3 hour fcst:
1.36:0:d=2018070106:PRES:surface:3 hour fcst:
1.37:0:d=2018070106:UGRD:10 m above ground:3 hour fcst:
1.38:0:d=2018070106:VGRD:10 m above ground:3 hour fcst:
1.39:0:d=2018070106:TMP:1.5 m above ground:3 hour fcst:
1.40:0:d=2018070106:RH:1.5 m above ground:3 hour fcst:
1.41:0:d=2018070106:LCDC:surface:3 hour fcst:
1.42:0:d=2018070106:MCDC:surface:3 hour fcst:
1.43:0:d=2018070106:HCDC:surface:3 hour fcst:
1.44:0:d=2018070106:TCDC:surface:3 hour fcst:
1.45:0:d=2018070106:APCP:surface:2-3 hour acc fcst:
1.46:0:d=2018070106:DSWRF:surface:2-3 hour ave fcst:
1.47:0:d=2018070106:PRMSL:mean sea level:4 hour fcst:
1.48:0:d=2018070106:PRES:surface:4 hour fcst:
1.49:0:d=2018070106:UGRD:10 m above ground:4 hour fcst:
1.50:0:d=2018070106:VGRD:10 m above ground:4 hour fcst:
1.51:0:d=2018070106:TMP:1.5 m above ground:4 hour fcst:
1.52:0:d=2018070106:RH:1.5 m above ground:4 hour fcst:
1.53:0:d=2018070106:LCDC:surface:4 hour fcst:
1.54:0:d=2018070106:MCDC:surface:4 hour fcst:
1.55:0:d=2018070106:HCDC:surface:4 hour fcst:
1.56:0:d=2018070106:TCDC:surface:4 hour fcst:
1.57:0:d=2018070106:APCP:surface:3-4 hour acc fcst:
1.58:0:d=2018070106:DSWRF:surface:3-4 hour ave fcst:
1.59:0:d=2018070106:PRMSL:mean sea level:5 hour fcst:
1.60:0:d=2018070106:PRES:surface:5 hour fcst:
1.61:0:d=2018070106:UGRD:10 m above ground:5 hour fcst:
1.62:0:d=2018070106:VGRD:10 m above ground:5 hour fcst:
1.63:0:d=2018070106:TMP:1.5 m above ground:5 hour fcst:
1.64:0:d=2018070106:RH:1.5 m above ground:5 hour fcst:
1.65:0:d=2018070106:LCDC:surface:5 hour fcst:
1.66:0:d=2018070106:MCDC:surface:5 hour fcst:
1.67:0:d=2018070106:HCDC:surface:5 hour fcst:
1.68:0:d=2018070106:TCDC:surface:5 hour fcst:
1.69:0:d=2018070106:APCP:surface:4-5 hour acc fcst:
1.70:0:d=2018070106:DSWRF:surface:4-5 hour ave fcst:
1.71:0:d=2018070106:PRMSL:mean sea level:6 hour fcst:
1.72:0:d=2018070106:PRES:surface:6 hour fcst:
1.73:0:d=2018070106:UGRD:10 m above ground:6 hour fcst:
1.74:0:d=2018070106:VGRD:10 m above ground:6 hour fcst:
1.75:0:d=2018070106:TMP:1.5 m above ground:6 hour fcst:
1.76:0:d=2018070106:RH:1.5 m above ground:6 hour fcst:
1.77:0:d=2018070106:LCDC:surface:6 hour fcst:
1.78:0:d=2018070106:MCDC:surface:6 hour fcst:
1.79:0:d=2018070106:HCDC:surface:6 hour fcst:
1.80:0:d=2018070106:TCDC:surface:6 hour fcst:
1.81:0:d=2018070106:APCP:surface:5-6 hour acc fcst:
1.82:0:d=2018070106:DSWRF:surface:5-6 hour ave fcst:
1.83:0:d=2018070106:PRMSL:mean sea level:7 hour fcst:
1.84:0:d=2018070106:PRES:surface:7 hour fcst:
1.85:0:d=2018070106:UGRD:10 m above ground:7 hour fcst:
1.86:0:d=2018070106:VGRD:10 m above ground:7 hour fcst:
1.87:0:d=2018070106:TMP:1.5 m above ground:7 hour fcst:
1.88:0:d=2018070106:RH:1.5 m above ground:7 hour fcst:
1.89:0:d=2018070106:LCDC:surface:7 hour fcst:
1.90:0:d=2018070106:MCDC:surface:7 hour fcst:
1.91:0:d=2018070106:HCDC:surface:7 hour fcst:
1.92:0:d=2018070106:TCDC:surface:7 hour fcst:
1.93:0:d=2018070106:APCP:surface:6-7 hour acc fcst:
1.94:0:d=2018070106:DSWRF:surface:6-7 hour ave fcst:
1.95:0:d=2018070106:PRMSL:mean sea level:8 hour fcst:
1.96:0:d=2018070106:PRES:surface:8 hour fcst:
1.97:0:d=2018070106:UGRD:10 m above ground:8 hour fcst:
1.98:0:d=2018070106:VGRD:10 m above ground:8 hour fcst:
1.99:0:d=2018070106:TMP:1.5 m above ground:8 hour fcst:
1.100:0:d=2018070106:RH:1.5 m above ground:8 hour fcst:
1.101:0:d=2018070106:LCDC:surface:8 hour fcst:
1.102:0:d=2018070106:MCDC:surface:8 hour fcst:
1.103:0:d=2018070106:HCDC:surface:8 hour fcst:
1.104:0:d=2018070106:TCDC:surface:8 hour fcst:
1.105:0:d=2018070106:APCP:surface:7-8 hour acc fcst:
1.106:0:d=2018070106:DSWRF:surface:7-8 hour ave fcst:
1.107:0:d=2018070106:PRMSL:mean sea level:9 hour fcst:
1.108:0:d=2018070106:PRES:surface:9 hour fcst:
1.109:0:d=2018070106:UGRD:10 m above ground:9 hour fcst:
1.110:0:d=2018070106:VGRD:10 m above ground:9 hour fcst:
1.111:0:d=2018070106:TMP:1.5 m above ground:9 hour fcst:
1.112:0:d=2018070106:RH:1.5 m above ground:9 hour fcst:
1.113:0:d=2018070106:LCDC:surface:9 hour fcst:
1.114:0:d=2018070106:MCDC:surface:9 hour fcst:
1.115:0:d=2018070106:HCDC:surface:9 hour fcst:
1.116:0:d=2018070106:TCDC:surface:9 hour fcst:
1.117:0:d=2018070106:APCP:surface:8-9 hour acc fcst:
1.118:0:d=2018070106:DSWRF:surface:8-9 hour ave fcst:
1.119:0:d=2018070106:PRMSL:mean sea level:10 hour fcst:
1.120:0:d=2018070106:PRES:surface:10 hour fcst:
1.121:0:d=2018070106:UGRD:10 m above ground:10 hour fcst:
1.122:0:d=2018070106:VGRD:10 m above ground:10 hour fcst:
1.123:0:d=2018070106:TMP:1.5 m above ground:10 hour fcst:
1.124:0:d=2018070106:RH:1.5 m above ground:10 hour fcst:
1.125:0:d=2018070106:LCDC:surface:10 hour fcst:
1.126:0:d=2018070106:MCDC:surface:10 hour fcst:
1.127:0:d=2018070106:HCDC:surface:10 hour fcst:
1.128:0:d=2018070106:TCDC:surface:10 hour fcst:
1.129:0:d=2018070106:APCP:surface:9-10 hour acc fcst:
1.130:0:d=2018070106:DSWRF:surface:9-10 hour ave fcst:
1.131:0:d=2018070106:PRMSL:mean sea level:11 hour fcst:
1.132:0:d=2018070106:PRES:surface:11 hour fcst:
1.133:0:d=2018070106:UGRD:10 m above ground:11 hour fcst:
1.134:0:d=2018070106:VGRD:10 m above ground:11 hour fcst:
1.135:0:d=2018070106:TMP:1.5 m above ground:11 hour fcst:
1.136:0:d=2018070106:RH:1.5 m above ground:11 hour fcst:
1.137:0:d=2018070106:LCDC:surface:11 hour fcst:
1.138:0:d=2018070106:MCDC:surface:11 hour fcst:
1.139:0:d=2018070106:HCDC:surface:11 hour fcst:
1.140:0:d=2018070106:TCDC:surface:11 hour fcst:
1.141:0:d=2018070106:APCP:surface:10-11 hour acc fcst:
1.142:0:d=2018070106:DSWRF:surface:10-11 hour ave fcst:
1.143:0:d=2018070106:PRMSL:mean sea level:12 hour fcst:
1.144:0:d=2018070106:PRES:surface:12 hour fcst:
1.145:0:d=2018070106:UGRD:10 m above ground:12 hour fcst:
1.146:0:d=2018070106:VGRD:10 m above ground:12 hour fcst:
1.147:0:d=2018070106:TMP:1.5 m above ground:12 hour fcst:
1.148:0:d=2018070106:RH:1.5 m above ground:12 hour fcst:
1.149:0:d=2018070106:LCDC:surface:12 hour fcst:
1.150:0:d=2018070106:MCDC:surface:12 hour fcst:
1.151:0:d=2018070106:HCDC:surface:12 hour fcst:
1.152:0:d=2018070106:TCDC:surface:12 hour fcst:
1.153:0:d=2018070106:APCP:surface:11-12 hour acc fcst:
1.154:0:d=2018070106:DSWRF:surface:11-12 hour ave fcst:
1.155:0:d=2018070106:PRMSL:mean sea level:13 hour fcst:
1.156:0:d=2018070106:PRES:surface:13 hour fcst:
1.157:0:d=2018070106:UGRD:10 m above ground:13 hour fcst:
1.158:0:d=2018070106:VGRD:10 m above ground:13 hour fcst:
1.159:0:d=2018070106:TMP:1.5 m above ground:13 hour fcst:
1.160:0:d=2018070106:RH:1.5 m above ground:13 hour fcst:
1.161:0:d=2018070106:LCDC:surface:13 hour fcst:
1.162:0:d=2018070106:MCDC:surface:13 hour fcst:
1.163:0:d=2018070106:HCDC:surface:13 hour fcst:
1.164:0:d=2018070106:TCDC:surface:13 hour fcst:
1.165:0:d=2018070106:APCP:surface:12-13 hour acc fcst:
1.166:0:d=2018070106:DSWRF:surface:12-13 hour ave fcst:
1.167:0:d=2018070106:PRMSL:mean sea level:14 hour fcst:
1.168:0:d=2018070106:PRES:surface:14 hour fcst:
1.169:0:d=2018070106:UGRD:10 m above ground:14 hour fcst:
1.170:0:d=2018070106:VGRD:10 m above ground:14 hour fcst:
1.171:0:d=2018070106:TMP:1.5 m above ground:14 hour fcst:
1.172:0:d=2018070106:RH:1.5 m above ground:14 hour fcst:
1.173:0:d=2018070106:LCDC:surface:14 hour fcst:
1.174:0:d=2018070106:MCDC:surface:14 hour fcst:
1.175:0:d=2018070106:HCDC:surface:14 hour fcst:
1.176:0:d=2018070106:TCDC:surface:14 hour fcst:
1.177:0:d=2018070106:APCP:surface:13-14 hour acc fcst:
1.178:0:d=2018070106:DSWRF:surface:13-14 hour ave fcst:
1.179:0:d=2018070106:PRMSL:mean sea level:15 hour fcst:
1.180:0:d=2018070106:PRES:surface:15 hour fcst:
1.181:0:d=2018070106:UGRD:10 m above ground:15 hour fcst:
1.182:0:d=2018070106:VGRD:10 m above ground:15 hour fcst:
1.183:0:d=2018070106:TMP:1.5 m above ground:15 hour fcst:
1.184:0:d=2018070106:RH:1.5 m above ground:15 hour fcst:
1.185:0:d=2018070106:LCDC:surface:15 hour fcst:
1.186:0:d=2018070106:MCDC:surface:15 hour fcst:
1.187:0:d=2018070106:HCDC:surface:15 hour fcst:
1.188:0:d=2018070106:TCDC:surface:15 hour fcst:
1.189:0:d=2018070106:APCP:surface:14-15 hour acc fcst:
1.190:0:d=2018070106:DSWRF:surface:14-15 hour ave fcst:

　このGRIB ファイルの中には非常に多くの気象データが格納されていて、総数は190 に上ることがわかります。

　wgrib2で取り出したインベントリーは、文字列がコロン「：」で区切られた形をして、いくつかは繰り返し登場します。そこで、Pythonに少し働いてもらい、どのような文字列が使われているかを調べてみようと思います。またこの際、関数 run に渡すコマンドライン文字列がとても長くてみずらいので、変数等の力も借りてスクリプトを少し見通し良くします。

　以下を実行してください。

In [4]:
# Script 4 #

wgrib2 = "c:/wgrib2/wgrib2"
# wgrib2 = "~/work/grib2/wgrib2/wgrib2"　# Macの場合　

grbdir = "jmadata/msm/2018/201807"
grbfile = "Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
grbpath = grbdir +"/"+ grbfile

rc = subprocess.run(f'{wgrib2} {grbpath}', 
                    shell=True, text=True, capture_output=True)
msgs = []
for line in rc.stdout.splitlines():
    msg = line.split(":")[1:-1]
    msgs.append(msg)
for i in range(len(msgs[0][:])):
    kw = [msg[i] for msg in msgs]    
    print(sorted(set(kw)))
    print("-"*10)

['0']
----------
['d=2018070106']
----------
['APCP', 'DSWRF', 'HCDC', 'LCDC', 'MCDC', 'PRES', 'PRMSL', 'RH', 'TCDC', 'TMP', 'UGRD', 'VGRD']
----------
['1.5 m above ground', '10 m above ground', 'mean sea level', 'surface']
----------
['0-1 hour acc fcst', '0-1 hour ave fcst', '1 hour fcst', '1-2 hour acc fcst', '1-2 hour ave fcst', '10 hour fcst', '10-11 hour acc fcst', '10-11 hour ave fcst', '11 hour fcst', '11-12 hour acc fcst', '11-12 hour ave fcst', '12 hour fcst', '12-13 hour acc fcst', '12-13 hour ave fcst', '13 hour fcst', '13-14 hour acc fcst', '13-14 hour ave fcst', '14 hour fcst', '14-15 hour acc fcst', '14-15 hour ave fcst', '15 hour fcst', '2 hour fcst', '2-3 hour acc fcst', '2-3 hour ave fcst', '3 hour fcst', '3-4 hour acc fcst', '3-4 hour ave fcst', '4 hour fcst', '4-5 hour acc fcst', '4-5 hour ave fcst', '5 hour fcst', '5-6 hour acc fcst', '5-6 hour ave fcst', '6 hour fcst', '6-7 hour acc fcst', '6-7 hour ave fcst', '7 hour fcst', '7-8 hour acc fcst', '7-8 hour ave fcst', '8 hour fcst', '8-9 hour acc fcst', '8-9 hour ave fcst', '9 hour fcst', '9-10 hour acc fcst', '9-10 hour ave fcst', 'anl']
----------

　インベントリーは、どうやら、連番、ゼロ、初期値の日時、気象要素記号、高度、予報時間でできているようです。

　GRIB2ファイルに対する処理は、「ファイルに格納されている どのデータに、何をする 」と一般的化して表現できます。wgrib2のオプション -match_fs "検索文字列" は、どのデータに を指定する役割を持ちます。より複雑な検索に対応できるよう、オプション -match "検索正規表現" も用意されています。このオプションは複数個使うことができ、「 且つ 」として扱われます。
「何をする」は極めて多様ですが、最も簡単なものとしては、より詳細にデータの内容表示させる -V があります。

　そこで、ここでは、インベントリの183番目(下から6番目)のデータの詳細情報を表示させて見ることにしましょう。下のスクリプトの8行目が「番頭1.183のデータの詳細を表示させる」を指示するコマンドラインオプションを定義する文です。

　それではこのCellを実行してください。

In [5]:
# Script 5 #

# wgrib2のファイル
wgrib2 = "c:/wgrib2/wgrib2"
#wgrib2 = "~/work/grib2/wgrib2/wgrib2"  # Macの場合

# オプション（操作の指示）
kwds = '-match_fs "1.183:" -V' #文字列中に「"」を入れたいときは、クォーテーションに「'」を使います

# 対象のファイル
grbdir = "jmadata/msm/2018/201807"  # GRIB2ファイルの置き場所
grbfile = "Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"  # GRIB2ファイルの名前
grbpath = grbdir +"/"+ grbfile  # ファイルへのフルパス

# 実行
rc = subprocess.run(f"{wgrib2} {kwds} {grbpath}",
                    shell=True, text=True, capture_output=True)
for line in rc.stdout.splitlines():
    print(line)

1.183:0:vt=2018070121:1.5 m above ground:15 hour fcst:TMP Temperature [K]:
    ndata=242905:undef=0:mean=296.075:min=278.445:max=303.39
    grid_template=0:winds(N/S):
	lat-lon grid:(481 x 505) units 1e-06 input WE:NS output WE:SN res 48
	lat 47.600000 to 22.400000 by 0.050000
	lon 120.000000 to 150.000000 by 0.062500 #points=242905

　詳細情報を読むと、このデータは、初期時刻2018年7月1日06UTC時の15時間後である、2018年7月1日21UTC時における地上1.5 mの気温を予測したもので、単位はK(ケルビン)、であり、北緯22.4～47.6度、東経120～150度の範囲を、緯度方向に0.05度、経度方向には0.0625度の間隔で張った格子点データであることがわかります。

1.2.2 GRIB2ファイルからのデータの取り出し¶

　オプション -csv 出力ファイル名 は、データをCSV形式のファイルとして出力させるためのオプションです。

　それでは、インベントリ番号1.183のデータをCSV ファイルに出力してみましょう。以下のCellを実行してください。フォルダchallenge に、ファイルTMP.csv が新規作成されているはずです。

In [6]:
# Script 6 #

# wgrib2のファイル
wgrib2 = "c:/wgrib2/wgrib2"  # Windowsの場合
# wgrib2 = "~/work/grib2/wgrib2/wgrib2"  # Macの場合

# オプション（操作の指示）
kwds = '-match_fs "1.183:" -csv TMP.csv'

# 対象のファイル
grbdir = "jmadata/msm/2018/201807"  # GRIB2ファイルの置き場所
grbfile = "Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"  # GRIB2ファイルの名前
grbpath = grbdir +"/"+ grbfile  # ファイルへのフルパス

rc = subprocess.run(f"{wgrib2} {kwds} {grbpath}",
                    shell=True, text=True, capture_output=True)
for line in rc.stdout.splitlines():
    print(line)

1.183:0:d=2018070106:TMP:1.5 m above ground:15 hour fcst:

　新規作成された TMP.csv を表計算ソフトやテキストエディタで表示し、中身を確認してください。下図のような内容になっているはずです。

image.png

　マイクロソフト・エクセルで表示させた方は、経度がE列、緯度がF列に示されます。シートの上の方にある低緯度の気温と、シート下の方のにある高緯度の気温を見比べて、低緯度の格子の気温の方が高緯度の格子の気温より高い傾向にあることを確認してください。

　コマンドラインオプション「 -match_fs 」で、インベントリに「 "TMP" 」が含まれるものを選択すれば、全ての予報期間の気温のデータをCSVファイルとして出力することができます。また、このオプションを用いずに、オプション -csv 出力ファイル名 を用いれば、全データがファイル出力されます。すなわち、GRIB2ファイルのフォーマット変換が行われることになります。

　wgrib2は、CSVのほかに、フラットバイナリやMySQL、NetCDFなどのフォーマットでデータを出力することができます。今度は、このファイルの全部のデータをNetCDFファイルで出力してみましょう。

　以下を実行してください。

In [7]:
# Script 7 #

# wgrib2のファイル
wgrib2 = "c:/wgrib2/wgrib2"  # Windowsの場合
# wgrib2 = "~/work/grib2/wgrib2/wgrib2"  # Macの場合

# 対象のファイル
grbdir = "jmadata/msm/2018/201807"  # GRIB2ファイルの置き場所
grbfile = "Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"  # GRIB2ファイルの名前
grbpath = grbdir +"/"+ grbfile  # ファイルへのフルパス

# オプション（操作の指示）
kwds = '-netcdf ' + grbfile + '.nc'

rc = subprocess.run(f"{wgrib2} {kwds} {grbpath}",
                    shell=True, text=True, capture_output=True)
#for line in rc.stdout.splitlines():
#    print(line)

　ファイルエクスプローラーで、フォルダ challenge の中を確認してください。長い名前で最後が「.nc 」で終わるファイルが新たにできているはずです。ファイルサイズを見ると、182Mバイトとあります。かなり大きいですね。

　wgrib2には、-match_fs や -V、-csv 以外にもたくさんのコマンドラインオプションがあるので、それらを表示させてみましょう。以下を実行してください。

In [8]:
# Script 8 #

# wgrib2のファイル
wgrib2 = "c:/wgrib2/wgrib2"  # Windowsの場合
# wgrib2 = "~/work/grib2/wgrib2/wgrib2"  # Macの場合

# オプション（操作の指示）
kwds = '-h'


rc = subprocess.run(f"{wgrib2} {kwds} ",
                    shell=True, text=True, capture_output=True)
for line in rc.stdout.splitlines():
    print(line)

wgrib2 v3.0.2 3/2021  Wesley Ebisuzaki, Reinoud Bokhorst, John Howard, Jaakko Hyvﾃ､tti, Dusan Jovic, Daniel Lee, Kristian Nilssen, Karl Pfeiffer, Pablo Romero, Manfred Schwarb, Gregor Schee, Arlindo da Silva, Niklas Sondell, Sam Trahan, George Trojan, Sergey Varlamov
   stock build
 -else            else         else, -if ... -else ... -endif
 -elseif          elif  X      elseif X (POSIX regular expression) conditional on match, -if ... -elseif ... -endif
 -elseif_fs       elif  X      elseif X (fixed string) conditional execution
 -elseif_n        elif  X      elseif (inv numbers in range), X=(start:end:step)
 -elseif_rec      elif  X      elseif (record numbers in range), X=(start:end:step)
 -elseif_reg      elif  X      elseif rpn registers defined, X = A, A:B, A:B:C, etc A = register number
 -endif           endif        terminates if block
 -if              if    X      if X (POSIX regular expression), conditional execution on match
 -if_delayed_error if           if delayed error
 -if_fs           if    X      if X (fixed string), conditional execution on match
 -if_n            if    X      if (inv numbers in range), X=(start:end:step)
 -if_rec          if    X      if (record numbers in range), X=(start:end:step)
 -if_reg          if    X      if rpn registers defined, X = A, A:B, A:B:C, etc A = register number
 -not_if          if    X      not_if X (regular expression), conditional execution on not match
 -not_if_fs       if    X      if X (fixed string) does not match, conditional execution up to next output/fi
 -0xSec           inv   X      Hex dump of section X (0..8)
 -aerosol_size    inv          optical properties of an aerosol
 -aerosol_wavelength inv          optical properties of an aerosol
 -bitmap          inv          bitmap mode
 -center          inv          center
 -checksum        inv   X      CRC checksum of section X (0..8), whole message (X = -1/message) or (X=data)
 -disc            inv          discipline (code table 0.0)
 -domain          inv          find rectangular domain for g2ctl/GrADS plots
 -end_ft          inv          verf time = reference_time + forecast_time + stat. proc time (YYYYMMDDHH) (same as -vt)
 -end_FT          inv          verf time = reference_time + forecast_time + stat. proc time (YYYYMMDDHHMMSS) (same as -VT)
 -ens             inv          ensemble information
 -ext_name        inv          extended name, var+qualifiers
 -ftime           inv          either ftime1 or ftime2 dep on version_ftime
 -ftime1          inv          forecast time
 -ftime2          inv          timestamp -- will replace -ftime in the future TESTING
 -ftn_api_fn0     inv          n npnts nx ny msg_no submsg i11,5(1x,i11)
 -full_name       inv          extended name, var+misc+lev (depreciated)
 -gdt             inv          contents of Grid Definition Template (g2c)
 -geolocation     inv          package (proj4,gctpc,internal,not_used) to get lat/lon of grid points
 -get_byte        inv   X Y Z  get bytes in Section X, Octet Y, number of bytes Z (decimal format)
 -get_hex         inv   X Y Z  get bytes in Section X, Octet Y, number of bytes Z (bytes in hexadecimal format)
 -get_ieee        inv   X Y Z  get ieee float in Section X, Octet Y, number of floats Z
 -get_int         inv   X Y Z  get 4-byte ints in Section X, Octet Y, number of ints Z
 -get_int2        inv   X Y Z  get 2-byte ints in Section X, Octet Y, number of ints Z
 -grib_max_bits   inv          maximum bits used in grib encoding
 -grid            inv          grid definition
 -grid_id         inv          show values from grid_id
 -hybrid          inv          shows vertical coordinate parameters from Sec4
 -ij              inv   X Y    value of field at grid(X,Y) X=1,..,nx Y=1,..,ny (WxText enabled)
 -ijlat           inv   X Y    lat,lon and grid value at grid(X,Y) X=1,..,nx Y=1,..,ny (WxText enabled)
 -ilat            inv   X      lat,lon and grid value at Xth grid point, X=1,..,npnts (WxText enabled)
 -JMA             inv          inventory for JMA locally defined PDT
 -lev             inv          level (code table 4.5)
 -ll2i            inv   X Y    x=lon y=lat, converts to (i), 1..ndata
 -ll2ij           inv   X Y    x=lon y=lat, converts lon-lat to (i,j) using gctpc
 -lon             inv   X Y    value at grid point nearest lon=X lat=Y (WxText enabled)
 -match_inv       inv          inventory used by -match, -not, -if and -not_if
 -Match_inv       inv          same as -match_inv except d=YYYYMMDDHH <-> D=YYYYMMDDHHmmss
 -max             inv          print maximum value
 -min             inv          print minimum value
 -misc            inv          variable name qualifiers like chemical, ensemble, probability, etc
 -MM              inv          reference time MM
 -model_version_date inv          prints model date code
 -n               inv          prints out inventory number
 -N_ens           inv          number of ensemble members
 -nl              inv          inserts new line into inventory
 -nlons           inv          number of longitudes for each latitude
 -npts            inv          number of grid points
 -nxny            inv          nx and ny of grid
 -packing         inv          shows the packing mode (use -v for more details)
 -pdt             inv          Product Definition Table (Code Table 4.0)
 -precision       inv          precision of packing
 -print           inv   X      inserts string (X) into inventory
 -prob            inv          probability information
 -process         inv          Process (code table 4.3)
 -processid       inv          process id (locally defined)
 -proj4_ij2ll     inv   X Y    X=x Y=y, converts to (i,j) to lon-lat using proj.4 (experimental) we:sn
 -proj4_ll2i      inv   X Y    x=lon y=lat, converts to (i) using proj.4 (experimental) 1..ndata
 -proj4_ll2ij     inv   X Y    x=lon y=lat, converts lon-lat (i,j) using proj.4 (experimental)
 -pyinv           inv          miscelaneous metadata for pywgrib2_XXX (experimental)
 -radius          inv          radius of Earth
 -range           inv          print out location of record in bytes, 0 = first byte
 -reset_delayed_error inv          clear reset_delayed_error flag
 -RT              inv          type of reference Time
 -s               inv          simple inventory
 -S               inv          simple inventory with minutes and seconds (subject to change)
 -s2              inv          simple inventory .. for testing ftime2
 -scale           inv          scale for packing
 -scaling         inv          scaling for packing (old format)
 -scan            inv          scan order of grid
 -Sec0            inv          contents of section0
 -Sec3            inv          contents of section 3 (Grid Definition Section)
 -Sec4            inv          Sec 4 values (Product definition section)
 -Sec5            inv          Sec 5 values (Data representation section)
 -Sec6            inv          show bit-map section
 -Sec_len         inv          length of various grib sections
 -spatial_proc    inv          show spacial processing, pdt=4.15
 -spectral_bands  inv          spectral bands for satellite, pdt=4.31 or 4.32
 -start_ft        inv          verf time = reference_time + forecast_time (YYYYMMDDHH) : no stat. proc time
 -start_FT        inv          verf time = reference_time + forecast_time (YYYYMMDDHHMMSS) - no stat. proc time
 -stats           inv          statistical summary of data values
 -subcenter       inv          subcenter
 -t               inv          reference time YYYYMMDDHH, -v2 for alt format
 -T               inv          reference time YYYYMMDDHHMMSS
 -table           inv          parameter table
 -timer           inv          reads OpenMP timer
 -unix_time       inv          print unix timestamp for rt & vt
 -V               inv          diagnostic output
 -var             inv          short variable name
 -varX            inv          raw variable name - discipline mastertab localtab center parmcat parmnum
 -vector_dir      inv          grid or earth relative winds
 -verf            inv          simple inventory using verification time
 -vt              inv          verf time = reference_time + forecast_time, -v2 for alt format
 -VT              inv          verf time = reference_time + forecast_time (YYYYMMDDHHMMSS)
 -warn_old_g2     inv          warn if old g2lib would have problem
 -wave_partition  inv          ocean surface wave partition (pdt=4.52)
 -YY              inv          reference time YYYY
 -inv_f77         inv>  X Y Z  match inventory written to Z with character*(Y) and X=(bin,ieee)
 -last            inv>  X      write last inv item to file X
 -last0           inv>  X      write last inv item to beginning of file X
 -nl_out          inv>  X      write new line in file X
 -print_out       inv>  X Y    prints string (X) in file (Y)
 -s_out           inv>  X      simple inventory written to X
 -big_endian      misc         sets ieee output to big endian (default is big endian)
 -box_ave         misc  X Y Z  box average X=odd integer (lon) Y=odd integer (lat) critical_weight
 -check_pdt_size  misc  X      check pdt size X=1 enable/default, X=0 disable
 -colon           misc  X      replace item deliminator (:) with X
 -config          misc         shows the configuration
 -count           misc         prints count, number times this -count was processed
 -end             misc         stop after first (sub)message (save time)
 -error_final     misc  X Y Z  error if at end X=count Y=ne,eq,le,lt,gt,ge Z=integer
 -export_lonlat   misc  X      save lon-lat data in binary file
 -fix_CFSv2_fcst  misc  X Y Z  fixes CFSv2 monthly fcst X=daily or 00/06/12/18 Y=pert no. Z=number ens fcsts v1.0
 -fix_ncep        misc         fix ncep PDT=8 headers produced by cnvgrib
 -gctpc           misc  X       X=0,1 use gctpc library (default=1)
 -grid_changes    misc         prints number of grid changes
 -grid_def        misc         read lon and lat data from grib file -- experimental
 -h               misc         help, shows common options
 -header          misc         f77 header or nx-ny header in text output (default)
 -help            misc  X      help [search string|all], -help all, shows all options
 -ijundefine      misc  X Y Z  sets grid point values to undefined X=(in-box|out-box) Y=ix0:ix1 Z=iy0:iy1 ix=(1..nx) iy=(1..ny)
 -import_bin      misc  X      read binary file (X) for data
 -import_grib     misc  X      read grib2 file (X) for data
 -import_grib_fs  misc  X Y    read grib2 file (Y) sequentially for record that matches X (fixed string)
 -import_ieee     misc  X      read ieee file (X) for data
 -import_lonlat   misc  X      read lon-lat data from binary file
 -import_netcdf   misc  X Y Z  TESTING X=file Y=var Z=hyper-cube specification
 -import_text     misc  X      read text file (X) for data
 -limit           misc  X      stops after X fields decoded
 -little_endian   misc         sets ieee output to little endian (default is big endian)
 -mem_del         misc  X      delete mem file X
 -mem_final       misc  X Y    write mem file X to file Y at cleanup step
 -mem_init        misc  X Y    read mem file X from file Y (on initialization)
 -ndate           misc  X Y    X=date Y=dt print date + dt
 -ndates          misc  X Y Z  X=date0 Y=(date1|dt1) Z=dt2 for (date=date0; date<(date1|date0+dt1); date+=dt2) print date
 -ndates_fmt      misc  X      X = C format for ndates option
 -new_grid_format misc  X      new_grid output format X=bin,ieee,grib
 -new_grid_interpolation misc  X      new_grid interpolation X=bilinear,bicubic,neighbor,budget
 -new_grid_ipopt  misc  X      new_grid ipopt values X=i1:i2..:iN N <= 20
 -new_grid_vectors misc  X      change fields to vector interpolate: X=none,default,UGRD:VGRD,(U:V list)
 -new_grid_winds  misc  X      new_grid wind orientation: X = grid, earth (no default)
 -no_header       misc         no f77 header or nx-ny header in text output
 -proj4           misc  X      X=0,1 use proj4 library for geolocation (testing)
 -read_sec        misc  X Y    read grib message section (0-8) X from binary file (Y)
 -rewind_final    misc  X      rewinds file X on cleanup step if already opened, CW2
 -rewind_proc     misc  X      rewinds file X on processing step if already opened, CW2
 -rpn             misc  X      reverse polish notation calculator
 -rpn_rcl         misc  X      data = register X .. same as -rpn rcl_X .. no geolocation calc needed
 -rpn_sto         misc  X      register X = data.. same as -rpn sto_X .. no geolocation calc needed
 -scaling_0001    misc         changes scaling testing (sample)
 -set             misc  X Y    set X = Y, X=local_table,etc (help: -set help help)
 -set_ave         misc  X      set ave/acc .. only use on pdt=4.0/4.8 (old code)
 -set_bin_prec    misc  X      X use X bits and ECMWF-style grib encoding
 -set_bitmap      misc  X      use bitmap when creating complex packed files X=1/0
 -set_byte        misc  X Y Z  set bytes in Section X, Octet Y, bytes Z (a|a:b:c)
 -set_date        misc  X      changes date code, X=(+|-)N(hr|dy|mo|yr), YYYYMMDDHHmmSS
 -set_ensm_derived_fcst misc  X Y    convert PDT 0,1,2 -> 2, 8,11,12 -> 12, X=code table 4.7 Y=num ens members
 -set_ens_num     misc  X Y Z  convert PDT 0,1 -> 1, 8,11 -> 11, X=code table 4.6 Y=pert num Z=num ens members -1=No Change
 -set_ftime       misc  X      either set_ftime1 or set_ftime2 dep on version_ftime
 -set_ftime1      misc  X      set ftime
 -set_ftime2      misc  X      set ftime2 .. will be replace -set_ftime/ave in the future -- TESTING ---
 -set_gds         misc  X      makes new gds (section 3), X=size in bytes
 -set_grib_max_bits misc  X      sets scaling so number of bits does not exceed N in (new) grib output
 -set_grib_type   misc  X      set grib type = jpeg, simple, ieee, complex(1|2|3), aec, same
 -set_hex         misc  X Y Z  set bytes in Section X, Octet Y, bytes Z (a|a:b:c|abc) in hexadecimal
 -set_ieee        misc  X Y Z  set ieee float in Section X, Octet Y, floats Z (a|a:b:c)
 -set_ijval       misc  X Y Z  sets grid point value X=ix Y=iy Z=val
 -set_int         misc  X Y Z  set 4-byte ints in Section X, Octet Y, signed integers Z (a|a:b:c)
 -set_int2        misc  X Y Z  set 2-byte ints in Section X, Octet Y, signed integers Z (a|a:b:c)
 -set_ival        misc  X Y    sets grid point value X=i1:i2:.. Y=va1:val2:.. grid[i1] = val1,etc i>0
 -set_lev         misc  X      changes level code .. not complete
 -set_metadata    misc  X      read meta-data for grib writing from file X
 -set_metadata_str misc  X      X = metadata string
 -set_pdt         misc  X      makes new pdt, X=(+)PDT_number or X=(+)PDT_number:size of PDT in octets, +=copy metadata
 -set_percentile  misc  X      convert PDT 0..6 -> 6, 8..15 -> 10, X=percentile (0..100)
 -set_prob        misc  5 args X/Y forecasts Z=Code Table 4.9 A=lower limit B=upper limit
 -set_radius      misc  X      set radius of Earth X= 0,2,4,5,6,8,9 (Code Table 3.2), 1:radius , 7:major:minor
 -set_scaling     misc  X Y    set decimal scaling=X/same binary scaling=Y/same new grib messages
 -set_sec_size    misc  X Y    resizes section , X=section number, Y=size in octets, DANGEROUS
 -set_ts_dates    misc  X Y Z  changes date code for time series X=YYYYMMDDHH(mmss) Y=dtime Z=#msgs/date
 -set_var         misc  X      changes variable name
 -start_timer     misc         starts OpenMP timer
 -status          misc  X      X X=file
 -submsg          misc  X      process submessage X (0=process all messages)
 -sys             misc  X      run system/shell command, X=shell command
 -text_col        misc  X      number of columns on text output
 -text_fmt        misc  X      format for text output (C)
 -udf             misc  X Y    run UDF, X=program+optional_args, Y=return file
 -udf_arg         misc  X Y    add grib-data to UDF argument file, X=file Y=name
 -undefine        misc  X Y Z  sets grid point values to undefined X=(in-box|out-box) Y=lon0:lon1 Z=lat0:lat1
 -undefine_val    misc  X      grid point set to undefined if X=val or X=low:high
 -v               misc         verbose (v=1)
 -v0              misc         not verbose (v=0)
 -v2              misc         really verbose (v=2)
 -version         misc         print version
 --version        misc         print version
 -AAIG            out          writes Ascii ArcInfo Grid file, lat-lon grid only (alpha)
 -AAIGlong        out          writes Ascii ArcInfo Grid file, lat-lon grid only long-name *.asc (alpha)
 -ave             out   X Y    average X=time step Y=output v2
 -ave0            out   X Y    average X=time step, Y=output grib file needs file is special order
 -ave_var         out   X Y    average/std dev/min/max X=time step, Y=output
 -bin             out   X      write binary data to X
 -cress_lola      out   X..Z,A lon-lat grid values X=lon0:nlon:dlon Y=lat0:nlat:dlat Z=file A=radius1:radius2:..:radiusN
 -csv             out   X      make comma separated file, X=file (WxText enabled)
 -csv_long        out   X      make comma separated file, X=file (WxText enabled)
 -cubeface2global out   X Y    write faces X as global cubed grid to Y: X=list of faces to exclude
 -ens_processing  out   X Y    ave/min/max/spread X=output Y=future use
 -fcst_ave        out   X Y    average X=time step Y=output v2
 -fcst_ave0       out   X Y    average X=time step, Y=output grib file needs file is special order
 -fi              out          depreceated, used in old IF structure
 -grib            out   X      writes GRIB record (one submessage) to X
 -GRIB            out   X      writes entire GRIB record (all submessages)
 -grib_ieee       out   X      writes data[] to X.grb, X.head, X.tail, and X.h
 -grib_out        out   X      writes decoded/modified data in grib-2 format to file X
 -grib_out_irr    out   X Y    writes irregular grid grib (GDT=130 not adopted) X=(all|defined) Y=(output file)
 -grib_out_irr2   out   5 args writes irregular grid grib GDT 101 X=npnts Y=grid_no Z=grid_ref A=UUID B=(output file)
 -gribtable_used  out   X      write out sample gribtable as derived from grib file, X=file
 -gridout         out   X      text file with grid: i j lat lon (1st record)
 -ieee            out   X      write (default:big-endian) IEEE data to X
 -ijbox           out   X..Z,A grid values in bounding box X=i1:i2[:di] Y=j1:j2[:dj] Z=file A=[bin|text|spread]
 -ijsmall_grib    out   X Y Z  make small domain grib file X=ix0:ix1 Y=iy0:iy1 Z=file
 -irr_grid        out   X Y Z  make irregular grid (GDT=130 not adopted), nearest neighbor, X=lon-lat list Y=radius (km) Z=output grib file
 -lola            out   X..Z,A lon-lat grid values X=lon0:nlon:dlon Y=lat0:nlat:dlat Z=file A=[bin|text|spread|grib]
 -merge_fcst      out   X Y    merge forecast ave/acc/min/max X=number to intervals to merge (0=every) Y=output grib file
 -mysql           out   5 args H=[host] U=[user] P=[password] D=[db] T=[table]
 -mysql_dump      out   7 args H=[host] U=[user] P=[password] D=[db] T=[table] W=[western_lons:0|1] PV=[remove unlikely:0|1]
 -mysql_speed     out   7 args H=[host] U=[user] P=[password] D=[db] T=[table] W=[western_lons:0|1] PV=[remove unlikely:0|1]
 -ncep_norm       out   X      normalize NCEP-type ave/acc X=output grib file
 -ncep_uv         out   X      combine U and V fields into one message like NCEP operations
 -netcdf          out   X      write netcdf data to X
 -new_grid        out   X..Z,A bilinear interpolate: X=projection Y=x0:nx:dx Z=y0:ny:dy A=grib_file alpha
 -new_grid_order  out   X Y    put in required order for -new_grid, X=out Y=out2 no matching vector
 -reduced_gaussian_grid out   X Y Z  reduced Gaussian grid, X=outputfile Y=-1 Z=(neighbor|linear)[-extrapolate]
 -small_grib      out   X Y Z  make small domain grib file X=lonW:lonE Y=latS:latN Z=file
 -spread          out   X      write text - spread sheet format into X (WxText enabled)
 -submsg_uv       out   X      combine vector fields into one message
 -text            out   X      write text data into X
 -time_processing out   X..Z,A average X=CodeTable 4.10 Y=CodeTable 4.11 Z=time step A=output
 -tosubmsg        out   X      convert GRIB message to submessage and write to file X
 -unmerge_fcst    out   X Y Z  unmerge_fcst X=output Y=fcst_time_0 Z: 0->result 1->+init 2->+all
 -wind_dir        out   X      calculate wind direction, X = output gribfile (direction in degrees, 0=wind from north, 90=wind from east)
 -wind_speed      out   X      calculate wind speed, X = output gribfile (U then V in datafile)
 -wind_uv         out   X      calculate UGRD/VGRD from speed/dir, X = output gribfile
 -write_sec       out   X Y    write grib msessage section X (0-8) to binary file Y
 -alarm           init  X      terminate after X seconds
 -append          init         append mode, write to existing output files
 -crlf            init         make the end of the inventory a crlf (windows) instead of newline (unix)
 -d               init  X      dump message X = n, n.m, n:offset, n.m:offset, only 1 -d allowed
 -egrep           init  X      egrep X | wgrib2 (X is POSIX regular expression)
 -egrep_v         init  X      egrep -v X | wgrib2 (X is POSIX regular expression)
 -eof_bin         init  X Y    send (binary) integer to file upon EOF: X=file Y=integer
 -eof_string      init  X Y    send string to file upon EOF: X=file Y=string
 -err_bin         init  X Y    send (binary) integer to file upon err exit: X=file Y=integer
 -err_string      init  X Y    send string to file upon err exit: X=file Y=string
 -fgrep           init  X      fgrep X | wgrib2
 -fgrep_v         init  X      fgrep -v X | wgrib2
 -fix_ncep_2      init         ncep bug fix 2, probability observation < -ve number
 -fix_ncep_3      init         sets flag to fix ncep bug 3 (constant fields)
 -fix_ncep_4      init         fixes NCEP grib2 files where DX and DY are undefined
 -fix_undef       init         set unused values to undef
 -for             init  X      process record numbers in range, X=(start:end:step), only one -for allowed
 -for_n           init  X      process inv numbers in range, X=(start:end:step), only one -for allowed
 -g2clib          init  X      X=0/1/2 0=WMO std 1=emulate g2clib 2=use g2clib
 -i               init         read Inventory from stdin
 -i_file          init  X      read Inventory from file
 -inv             init  X      write inventory to X
 -match           init  X      process data that matches X (POSIX regular expression)
 -match_fs        init  X      process data that matches X (fixed string)
 -match_inv_add   init  X Y Z  add new options to match_inventory
 -names           init  X      grib name convention, X=ecmwf, ncep
 -nc3             init         use netcdf3 (classic)
 -nc4             init         use netcdf4 (compressed, controlled endianness etc)
 -nc_grads        init         require netcdf file to be grads v1.9b4 compatible (fixed time step only)
 -nc_nlev         init  X      netcdf, X = max LEV dimension for {TIME,LEV,LAT,LON} data
 -nc_pack         init  X      pack/check limits of all NEW input variables, X=min:max[:byte|short|float]
 -ncpu            init  X      number of threads, default is environment variable OMP_NUM_THREADS/number of cpus
 -nc_table        init  X      X is conversion_to_netcdf_table file name
 -nc_time         init  X      netcdf, [[-]yyyymmddhhnnss]:[dt{s[ec]|m[in]|h[our]|d[ay]}], [-] is for time alignment only
 -no_append       init         not append mode, write to new output files (default)
 -no_nc_grads     init         netcdf file may be not grads v1.9b4 compatible, variable time step
 -no_nc_pack      init         no packing in netcdf for NEW variables
 -no_nc_table     init         disable previously defined conversion_to_netcdf_table
 -no_nc_time      init         netcdf, disable previously defined initial or relative date and time step
 -not             init  X      process data that does not match X (POSIX regular expression)
 -not_fs          init  X      process data that does not match X (fixed string)
 -one_line        init         puts all on one line (makes into inventory format)
 -order           init  X      decoded data in X (raw|we:sn|we:ns) order, we:sn is default
 -persistent      init  X      makes file X persistent if already opened (default on open), CW2
 -rewind_init     init  X      rewinds file X on initialization if already opened, CW2
 -set_ext_name    init  X      X=0/1 extended name on/off
 -set_ext_name_chars init  X Y    extended name characters X=field Y=space
 -set_regex       init  X      set regex mode X = 0:extended regex (default) 1:pattern 2:extended regex & quote metacharacters
 -set_version_ftime init  X      set version of ftime X=1, 2
 -tigge           init         use modified-TIGGE grib table
 -transient       init  X      make file X transient, CW2

1.3 WXBCオリジナルライブラリwxbcgribX¶

　WXBCオリジナルライブラリ wxbcgribX は、GRIB2ファイルからメッシュ気象データを取り出し、Pythonパッケージ Xarray が提供する処理機能を活用できるよう準備する機能を提供します。

  • ご注意
    Mac をお使いになる方は、フォルダ challenge 内のファイル wxbcgribx.py をエディタで開き、35行目の文をコメントアウトして以下のようにしてください。
wgrib2 = pathlib.Path("~/work/grib2/wgrib2/wgrib2") # Macの場合
# wgrib2 = pathlib.Path("C:/wgrib2/wgrib2.exe")       # Windowsの場合

　

1.3.1 Pythonパッケージ Xarray¶

　Pythonには、Xarray というライブラリパッケージが存在し、これをインポートすると、DataArray と呼ばれるデータ形式（コンピュータの用語ではオブジェクトクラス）を利用することができます。 DataArrayは、多次元のデータ配列、それが従う座標、データや座標に関する情報（属性と呼ばれる）を一纏めにしたオブジェクトです。

　気象データは、いつの、どの地点の(、どの高度の)ものであるかという情報と不可分であり、さらに、メッシュ気象データは、これが時空間上に規則正しく多数配置されたものなので、DataArrayは、Pythonでメッシュ気象データを取り扱うので最適なデータ型といえます。 さらに便利なことに、モジュール Xarray は、DataArrayの他に Dataset と呼ばれるオブジェクトクラスも提供しています。Datasetを用いると、座標が同じ複数のDataArrayをひとまとめにすることができます。これは、気象庁の数値予報モデルGPVのような、ファイルの中に複数のメッシュ気象データが同梱されているプロダクトを取り扱うのに場合に役立ちます。

1.3.2 Datasetオブジェクト¶

　Xarrayは、NetCDFファイルとの連携を強く意識して作られており、NetCDFファイルのデータをとても簡単にDatasetオブジェクトすることができます。
1.2 の最後でメソ数値予報モデルGPVのデータをNetCDFファイルとして保存したので、これを読み込んでみましょう。

　まず、ライブラリ xarray を、短縮名 xr を付けてインポートします。

In [9]:
# Script 9 #

import xarray as xr

　NetCDFファイルからデータを読み込んで ds という名前のDataSetオブジェクトにするには、xarrayの関数 open_dataset を用いて以下のようにします。

ds = xr.open_dataset("NetCDFファイル名")

　以下を実行してください。

In [10]:
# Script 10 #
ds = xr.open_dataset("Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin.nc")
ds
Out[10]:
<xarray.Dataset>
Dimensions:              (latitude: 505, longitude: 481, time: 16)
Coordinates:
  * latitude             (latitude) float64 22.4 22.45 22.5 ... 47.5 47.55 47.6
  * longitude            (longitude) float64 120.0 120.1 120.1 ... 149.9 150.0
  * time                 (time) datetime64[ns] 2018-07-01T06:00:00 ... 2018-0...
Data variables:
    PRMSL_meansealevel   (time, latitude, longitude) float32 ...
    PRES_surface         (time, latitude, longitude) float32 ...
    UGRD_10maboveground  (time, latitude, longitude) float32 ...
    VGRD_10maboveground  (time, latitude, longitude) float32 ...
    TMP_1D5maboveground  (time, latitude, longitude) float32 ...
    RH_1D5maboveground   (time, latitude, longitude) float32 ...
    LCDC_surface         (time, latitude, longitude) float32 ...
    MCDC_surface         (time, latitude, longitude) float32 ...
    HCDC_surface         (time, latitude, longitude) float32 ...
    TCDC_surface         (time, latitude, longitude) float32 ...
    APCP_surface         (time, latitude, longitude) float32 ...
    DSWRF_surface        (time, latitude, longitude) float32 ...
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 505
    • longitude: 481
    • time: 16
    • latitude
      (latitude)
      float64
      22.4 22.45 22.5 ... 47.5 47.55 47.6
      units :
      degrees_north
      long_name :
      latitude
      array([22.4 , 22.45, 22.5 , ..., 47.5 , 47.55, 47.6 ])
    • longitude
      (longitude)
      float64
      120.0 120.1 120.1 ... 149.9 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([120.    , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150.    ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00 ... 2018-07-...
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      3
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      forecasts or accumulated (including analyses), reference date is fixed
      time_step_setting :
      auto
      time_step :
      3600.0
      array(['2018-07-01T06:00:00.000000000', '2018-07-01T07:00:00.000000000',
       '2018-07-01T08:00:00.000000000', '2018-07-01T09:00:00.000000000',
       '2018-07-01T10:00:00.000000000', '2018-07-01T11:00:00.000000000',
       '2018-07-01T12:00:00.000000000', '2018-07-01T13:00:00.000000000',
       '2018-07-01T14:00:00.000000000', '2018-07-01T15:00:00.000000000',
       '2018-07-01T16:00:00.000000000', '2018-07-01T17:00:00.000000000',
       '2018-07-01T18:00:00.000000000', '2018-07-01T19:00:00.000000000',
       '2018-07-01T20:00:00.000000000', '2018-07-01T21:00:00.000000000'],
      dtype='datetime64[ns]')
    • PRMSL_meansealevel
      (time, latitude, longitude)
      float32
      ...
      short_name :
      PRMSL_meansealevel
      long_name :
      Pressure Reduced to MSL
      level :
      mean sea level
      units :
      Pa
      [3886480 values with dtype=float32]
    • PRES_surface
      (time, latitude, longitude)
      float32
      ...
      short_name :
      PRES_surface
      long_name :
      Pressure
      level :
      surface
      units :
      Pa
      [3886480 values with dtype=float32]
    • UGRD_10maboveground
      (time, latitude, longitude)
      float32
      ...
      short_name :
      UGRD_10maboveground
      long_name :
      U-Component of Wind
      level :
      10 m above ground
      units :
      m/s
      [3886480 values with dtype=float32]
    • VGRD_10maboveground
      (time, latitude, longitude)
      float32
      ...
      short_name :
      VGRD_10maboveground
      long_name :
      V-Component of Wind
      level :
      10 m above ground
      units :
      m/s
      [3886480 values with dtype=float32]
    • TMP_1D5maboveground
      (time, latitude, longitude)
      float32
      ...
      short_name :
      TMP_1D5maboveground
      long_name :
      Temperature
      level :
      1.5 m above ground
      units :
      K
      [3886480 values with dtype=float32]
    • RH_1D5maboveground
      (time, latitude, longitude)
      float32
      ...
      short_name :
      RH_1D5maboveground
      long_name :
      Relative Humidity
      level :
      1.5 m above ground
      units :
      percent
      [3886480 values with dtype=float32]
    • LCDC_surface
      (time, latitude, longitude)
      float32
      ...
      short_name :
      LCDC_surface
      long_name :
      Low Cloud Cover
      level :
      surface
      units :
      percent
      [3886480 values with dtype=float32]
    • MCDC_surface
      (time, latitude, longitude)
      float32
      ...
      short_name :
      MCDC_surface
      long_name :
      Medium Cloud Cover
      level :
      surface
      units :
      percent
      [3886480 values with dtype=float32]
    • HCDC_surface
      (time, latitude, longitude)
      float32
      ...
      short_name :
      HCDC_surface
      long_name :
      High Cloud Cover
      level :
      surface
      units :
      percent
      [3886480 values with dtype=float32]
    • TCDC_surface
      (time, latitude, longitude)
      float32
      ...
      short_name :
      TCDC_surface
      long_name :
      Total Cloud Cover
      level :
      surface
      units :
      percent
      [3886480 values with dtype=float32]
    • APCP_surface
      (time, latitude, longitude)
      float32
      ...
      short_name :
      APCP_surface
      long_name :
      Total Precipitation
      level :
      surface
      units :
      kg/m^2
      [3886480 values with dtype=float32]
    • DSWRF_surface
      (time, latitude, longitude)
      float32
      ...
      short_name :
      DSWRF_surface
      long_name :
      Downward Short-Wave Radiation Flux
      level :
      surface
      units :
      W/m^2
      [3886480 values with dtype=float32]
    • latitude
      PandasIndex
      PandasIndex(Float64Index([              22.4,              22.45,               22.5,
                           22.55,               22.6,              22.65,
                            22.7,              22.75,               22.8,
                           22.85,
              ...
               47.14999999999999, 47.199999999999996,  47.24999999999999,
                            47.3, 47.349999999999994,  47.39999999999999,
              47.449999999999996,  47.49999999999999,              47.55,
              47.599999999999994],
             dtype='float64', name='latitude', length=505))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([   120.0, 120.0625,  120.125, 120.1875,   120.25, 120.3125,
               120.375, 120.4375,    120.5, 120.5625,
              ...
              149.4375,    149.5, 149.5625,  149.625, 149.6875,   149.75,
              149.8125,  149.875, 149.9375,    150.0],
             dtype='float64', name='longitude', length=481))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00', '2018-07-01 07:00:00',
               '2018-07-01 08:00:00', '2018-07-01 09:00:00',
               '2018-07-01 10:00:00', '2018-07-01 11:00:00',
               '2018-07-01 12:00:00', '2018-07-01 13:00:00',
               '2018-07-01 14:00:00', '2018-07-01 15:00:00',
               '2018-07-01 16:00:00', '2018-07-01 17:00:00',
               '2018-07-01 18:00:00', '2018-07-01 19:00:00',
               '2018-07-01 20:00:00', '2018-07-01 21:00:00'],
              dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0

　データの次元、保持されているデータの種類、数値形式などが整然と表示されました。wgrib2で見たときはインベントリーは190個ありましたが、それらは時間について連結されています。
リストの右端に文書を模ったアイコンがあります。これをクリックするとそれぞれのデータの名称や単位など、さらに詳細な属性情報を見ることができます。重なったお皿のアイコンをクリックすると、データの概要が表示されます。

1.3.3 DataArrayオブジェクト¶

　DataSetに格納される個々の気象データは、DataArray オブジェクトとしてとりだすことができ、様々に解析することができます。上で作成した DataSet オブジェクト ds から、気温のデータを取り出して、 da という名前のDataArray オブジェクトを作ってみましょう。

　DataSetオブジェクト ds から特定の気象要素を取り出して DataArray オブジェクトda を作るには、気象要素のハンドルネームを使用して以下のようにします。

da = ds["気象要素のハンドルネーム"]

　気象要素のハンドルネームは、上の表の「Data Variables」のセクションの左端の列に示されているものです。GRIB2ファイルから特定の気象要素抽出するのに使用した略称とは別なので注意してください。気温のハンドルネームは「TMP_1D5maboveground」です。

　実際に取り込んでみましょう。以下を実行してください。

In [11]:
# Script 11 #

da = ds["TMP_1D5maboveground"]
da
Out[11]:
<xarray.DataArray 'TMP_1D5maboveground' (time: 16, latitude: 505, longitude: 481)>
[3886480 values with dtype=float32]
Coordinates:
  * latitude   (latitude) float64 22.4 22.45 22.5 22.55 ... 47.5 47.55 47.6
  * longitude  (longitude) float64 120.0 120.1 120.1 120.2 ... 149.9 149.9 150.0
  * time       (time) datetime64[ns] 2018-07-01T06:00:00 ... 2018-07-01T21:00:00
Attributes:
    short_name:  TMP_1D5maboveground
    long_name:   Temperature
    level:       1.5 m above ground
    units:       K
xarray.DataArray
'TMP_1D5maboveground'
  • time: 16
  • latitude: 505
  • longitude: 481
  • ...
    [3886480 values with dtype=float32]
    • latitude
      (latitude)
      float64
      22.4 22.45 22.5 ... 47.5 47.55 47.6
      units :
      degrees_north
      long_name :
      latitude
      array([22.4 , 22.45, 22.5 , ..., 47.5 , 47.55, 47.6 ])
    • longitude
      (longitude)
      float64
      120.0 120.1 120.1 ... 149.9 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([120.    , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150.    ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00 ... 2018-07-...
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      3
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      forecasts or accumulated (including analyses), reference date is fixed
      time_step_setting :
      auto
      time_step :
      3600.0
      array(['2018-07-01T06:00:00.000000000', '2018-07-01T07:00:00.000000000',
       '2018-07-01T08:00:00.000000000', '2018-07-01T09:00:00.000000000',
       '2018-07-01T10:00:00.000000000', '2018-07-01T11:00:00.000000000',
       '2018-07-01T12:00:00.000000000', '2018-07-01T13:00:00.000000000',
       '2018-07-01T14:00:00.000000000', '2018-07-01T15:00:00.000000000',
       '2018-07-01T16:00:00.000000000', '2018-07-01T17:00:00.000000000',
       '2018-07-01T18:00:00.000000000', '2018-07-01T19:00:00.000000000',
       '2018-07-01T20:00:00.000000000', '2018-07-01T21:00:00.000000000'],
      dtype='datetime64[ns]')
    • latitude
      PandasIndex
      PandasIndex(Float64Index([              22.4,              22.45,               22.5,
                           22.55,               22.6,              22.65,
                            22.7,              22.75,               22.8,
                           22.85,
              ...
               47.14999999999999, 47.199999999999996,  47.24999999999999,
                            47.3, 47.349999999999994,  47.39999999999999,
              47.449999999999996,  47.49999999999999,              47.55,
              47.599999999999994],
             dtype='float64', name='latitude', length=505))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([   120.0, 120.0625,  120.125, 120.1875,   120.25, 120.3125,
               120.375, 120.4375,    120.5, 120.5625,
              ...
              149.4375,    149.5, 149.5625,  149.625, 149.6875,   149.75,
              149.8125,  149.875, 149.9375,    150.0],
             dtype='float64', name='longitude', length=481))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00', '2018-07-01 07:00:00',
               '2018-07-01 08:00:00', '2018-07-01 09:00:00',
               '2018-07-01 10:00:00', '2018-07-01 11:00:00',
               '2018-07-01 12:00:00', '2018-07-01 13:00:00',
               '2018-07-01 14:00:00', '2018-07-01 15:00:00',
               '2018-07-01 16:00:00', '2018-07-01 17:00:00',
               '2018-07-01 18:00:00', '2018-07-01 19:00:00',
               '2018-07-01 20:00:00', '2018-07-01 21:00:00'],
              dtype='datetime64[ns]', name='time', freq=None))
  • short_name :
    TMP_1D5maboveground
    long_name :
    Temperature
    level :
    1.5 m above ground
    units :
    K

　先ほどと同じような表が示されましたが、表の左上に、今後は「 xarray.DataArray 」と表示されています。また、「Data variables」 のセクションがないところが違っています。

　DataArrayオブジェクトは、データの数値だけでなくそれが従う座標と時刻、さらにデータの単位等もすべて一塊になっているので、例えば分布図を可視化するなどの時に威力を発揮します。DataArrayオブジェクトのメソッド plot を使って、15時間先の気温予報図を描いてみましょう。
da は、時刻、緯度、経度の次元からなる3次元データですが、分布図は緯度経度だけの２次元の図なので，da[15,:,:]として、時刻だけを固定し、緯度度経度は「すべて」を指定します。これにより、描画のメソッドには緯度経度の２次元のデータが渡され、分布図が描けます。

In [12]:
# Script 12 #

da[15,:,:].plot()
Out[12]:
<matplotlib.collections.QuadMesh at 0x1f928831030>
No description has been provided for this image

　取り出した格子の緯度経度や、横軸縦軸もしっかりと作られていますね。

　特定の位置の気温の時系列グラフも簡単に作ることができます。時系列図は時間方向のみの１次元の図なので，da[:,200,180]として緯度経度を固定し、時刻だけは「すべて」を指定します。これにより、描画のメソッドには時刻の１次元のデータが渡され、時系列図が描けます。
以下を実行してください。

In [13]:
# Script 13 #

da[:,200,180].plot.line(x="time",figsize=(4,3))
Out[13]:
[<matplotlib.lines.Line2D at 0x1f929561780>]
No description has been provided for this image

1.3.4 データ読み出し関数 getgpv¶

　ここまでの学習で、GRIB2ファイルに保存されているメッシュ気象データは、下図の手順を踏むことで、Pythonで処理できることが分かりました。

image.png

　図中、破線で囲まれた部分は、GRIB2ファイルで提供された気象データを扱う際には必ず踏まなければならない手順なので、関数にしてしまった方が便利です。そこで、気象ビジネス推進コンソーシアム人材育成ワーキンググルプ内勉強会「気象データ×IT勉強会」は、これまで学習したことに基づいて、この部分を一気に処理する関数 getgpv を作り、簡単に使えるよう、ライブラリ wxbcgribX に納めました。このライブラリに収められている関数を使うと、GRIB2ファイルからのデータの取り出しがとても簡単に行えます。なお、このライブラリの実体は、フォルダchallenge に置かれているテキストファイル wxbcgribx.py です。

　それでは、関数 getgpv の使い方を説明します。まず最初に、ライブラリを短縮名 wx でインポートします。

In [14]:
# Script 14 #

import  wxbcgribx as wx

ライブラリ wxbcgribx をインポートすると、関数 getgpv が利用できるようになります。この関数は、以下の構文で用います。

ds = wx.getgpv( grbpath, element, ncdir="./nc", to_netcdf=True, from_netcdf=True, lalomima=None, verbose=False)

　引数は以下のとおりです。

  • grbpath:（文字列または文字列のリスト）読み込むべきGRIB2ファイルへのパス
  • element:（文字列または文字列のリスト）読み出す気象要素　GRIB2ファイルのインベントリで使用されている略称で指定
  • to_netcdf: （True/False）作成したDataDetオブジェクトをNetCDFファイルとして保管する場合にTrueを指定、デフォルトはTrue（保管する）
  • ncdir: (文字）上記NetCDFファイルを保管する場所、デフォルトは「./nc」
  • from_netcdf: （True/False）以前に保管したNetCDFファイルがあればGRIB2ファイルでなくそれから読み込むときにTrueを指定、デフォルトはTrue（NetCDFを読みGRIB2からは読まない）
  • lalomima: （4要素リスト）緯度経度範囲を与えると、Datasetオブジェクトにする緯度経度範囲を指定する。指定しないと全範囲を取得する
  • verbose: （True/False）処理の進捗を表示するときにTrue　デフォルトはFalse（表示しない）
  • 戻り値: Xarray Dataset オブジェクト

注意：この関数をWindowsで使用する場合、日本語のフォルダ名を含む場所にファイルを保存できないことがあります。

　この関数を使って、気象庁MSM-GPVのGRIB2ファイルから気温と相対湿度のデータを改めて取り出してみましょう。

　関数 getgpv は、第1引数(gribpath)で指定したGRIB2ファイルから、第2引数(elements)で指定した略号に対応する気象データを取り出して、Datasetオブジェクトにするので、Script 15 のようにすると、MSM-GPVのGRIB2ファイルからデータが読み取られ、気温と湿度が格納されたDatasetオブジェクト ds が作られます。
下のCellを実行してください。

In [15]:
# Script 15 #

grbdir = "jmadata/msm/2018/201807"
grbfile = "Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
grbpath = grbdir +"/"+ grbfile
elements = ["TMP","RH"]

ds = wx.getgpv(grbpath, elements, to_netcdf=True, verbose=True)

ds

1.5:0:d=2018070106:TMP:1.5 m above ground:anl:
1.6:0:d=2018070106:RH:1.5 m above ground:anl:
1.15:0:d=2018070106:TMP:1.5 m above ground:1 hour fcst:
1.16:0:d=2018070106:RH:1.5 m above ground:1 hour fcst:
1.27:0:d=2018070106:TMP:1.5 m above ground:2 hour fcst:
1.28:0:d=2018070106:RH:1.5 m above ground:2 hour fcst:
1.39:0:d=2018070106:TMP:1.5 m above ground:3 hour fcst:
1.40:0:d=2018070106:RH:1.5 m above ground:3 hour fcst:
1.51:0:d=2018070106:TMP:1.5 m above ground:4 hour fcst:
1.52:0:d=2018070106:RH:1.5 m above ground:4 hour fcst:
1.63:0:d=2018070106:TMP:1.5 m above ground:5 hour fcst:
1.64:0:d=2018070106:RH:1.5 m above ground:5 hour fcst:
1.75:0:d=2018070106:TMP:1.5 m above ground:6 hour fcst:
1.76:0:d=2018070106:RH:1.5 m above ground:6 hour fcst:
1.87:0:d=2018070106:TMP:1.5 m above ground:7 hour fcst:
1.88:0:d=2018070106:RH:1.5 m above ground:7 hour fcst:
1.99:0:d=2018070106:TMP:1.5 m above ground:8 hour fcst:
1.100:0:d=2018070106:RH:1.5 m above ground:8 hour fcst:
1.111:0:d=2018070106:TMP:1.5 m above ground:9 hour fcst:
1.112:0:d=2018070106:RH:1.5 m above ground:9 hour fcst:
1.123:0:d=2018070106:TMP:1.5 m above ground:10 hour fcst:
1.124:0:d=2018070106:RH:1.5 m above ground:10 hour fcst:
1.135:0:d=2018070106:TMP:1.5 m above ground:11 hour fcst:
1.136:0:d=2018070106:RH:1.5 m above ground:11 hour fcst:
1.147:0:d=2018070106:TMP:1.5 m above ground:12 hour fcst:
1.148:0:d=2018070106:RH:1.5 m above ground:12 hour fcst:
1.159:0:d=2018070106:TMP:1.5 m above ground:13 hour fcst:
1.160:0:d=2018070106:RH:1.5 m above ground:13 hour fcst:
1.171:0:d=2018070106:TMP:1.5 m above ground:14 hour fcst:
1.172:0:d=2018070106:RH:1.5 m above ground:14 hour fcst:
1.183:0:d=2018070106:TMP:1.5 m above ground:15 hour fcst:
1.184:0:d=2018070106:RH:1.5 m above ground:15 hour fcst:
Out[15]:
<xarray.Dataset>
Dimensions:              (latitude: 505, longitude: 481, time: 16)
Coordinates:
  * latitude             (latitude) float64 22.4 22.45 22.5 ... 47.5 47.55 47.6
  * longitude            (longitude) float64 120.0 120.1 120.1 ... 149.9 150.0
  * time                 (time) datetime64[ns] 2018-07-01T06:00:00 ... 2018-0...
Data variables:
    TMP_1D5maboveground  (time, latitude, longitude) float32 302.1 ... 280.1
    RH_1D5maboveground   (time, latitude, longitude) float32 81.88 ... 98.99
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 505
    • longitude: 481
    • time: 16
    • latitude
      (latitude)
      float64
      22.4 22.45 22.5 ... 47.5 47.55 47.6
      units :
      degrees_north
      long_name :
      latitude
      array([22.4 , 22.45, 22.5 , ..., 47.5 , 47.55, 47.6 ])
    • longitude
      (longitude)
      float64
      120.0 120.1 120.1 ... 149.9 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([120.    , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150.    ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00 ... 2018-07-...
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      3
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      forecasts or accumulated (including analyses), reference date is fixed
      time_step_setting :
      auto
      time_step :
      3600.0
      array(['2018-07-01T06:00:00.000000000', '2018-07-01T07:00:00.000000000',
       '2018-07-01T08:00:00.000000000', '2018-07-01T09:00:00.000000000',
       '2018-07-01T10:00:00.000000000', '2018-07-01T11:00:00.000000000',
       '2018-07-01T12:00:00.000000000', '2018-07-01T13:00:00.000000000',
       '2018-07-01T14:00:00.000000000', '2018-07-01T15:00:00.000000000',
       '2018-07-01T16:00:00.000000000', '2018-07-01T17:00:00.000000000',
       '2018-07-01T18:00:00.000000000', '2018-07-01T19:00:00.000000000',
       '2018-07-01T20:00:00.000000000', '2018-07-01T21:00:00.000000000'],
      dtype='datetime64[ns]')
    • TMP_1D5maboveground
      (time, latitude, longitude)
      float32
      302.1 302.1 302.1 ... 280.2 280.1
      short_name :
      TMP_1D5maboveground
      long_name :
      Temperature
      level :
      1.5 m above ground
      units :
      K
      array([[[302.1476 , 302.1398 , 302.1398 , ..., 301.65543, 301.632  ,
         301.61636],
        [302.12418, 302.12418, 302.11636, ..., 301.6476 , 301.632  ,
         301.61636],
        [302.10855, 302.10074, 302.10074, ..., 301.6476 , 301.632  ,
         301.60855],
        ...,
        [300.1476 , 299.9523 , 299.9601 , ..., 279.55386, 279.4523 ,
         279.34293],
        [299.8898 , 299.7023 , 299.67105, ..., 279.55386, 279.43668,
         279.3195 ],
        [299.99136, 299.71793, 299.65543, ..., 279.55386, 279.43668,
         279.3195 ]],

       [[302.191  , 302.191  , 302.15195, ..., 301.6285 , 301.6207 ,
         301.60507],
        [302.15195, 302.13632, 302.11288, ..., 301.6285 , 301.6207 ,
         301.60507],
        [302.10507, 302.10507, 302.09726, ..., 301.6285 , 301.61288,
         301.59726],
...
        [289.45142, 289.04517, 288.86548, ..., 280.22485, 280.16235,
         280.09204],
        [289.7561 , 289.30298, 289.08423, ..., 280.24048, 280.17798,
         280.10767],
        [290.1858 , 289.78735, 289.38892, ..., 280.2561 , 280.1936 ,
         280.11548]],

       [[302.15546, 302.15546, 302.16327, ..., 301.53046, 301.48358,
         301.44452],
        [302.1789 , 302.17108, 302.17108, ..., 301.52264, 301.48358,
         301.4367 ],
        [302.19452, 302.1867 , 302.1867 , ..., 301.50702, 301.46796,
         301.4289 ],
        ...,
        [289.34296, 289.00702, 288.89764, ..., 280.21796, 280.16327,
         280.13202],
        [289.60858, 289.21014, 289.03046, ..., 280.22577, 280.1789 ,
         280.13202],
        [290.00702, 289.63983, 289.3039 , ..., 280.2492 , 280.19452,
         280.13983]]], dtype=float32)
    • RH_1D5maboveground
      (time, latitude, longitude)
      float32
      81.88 81.78 81.69 ... 98.94 98.99
      short_name :
      RH_1D5maboveground
      long_name :
      Relative Humidity
      level :
      1.5 m above ground
      units :
      percent
      array([[[81.87808 , 81.78433 , 81.69058 , ..., 76.59683 , 76.72183 ,
         76.87808 ],
        [81.94058 , 81.84683 , 81.81558 , ..., 76.62808 , 76.78433 ,
         76.94058 ],
        [81.94058 , 81.90933 , 81.90933 , ..., 76.65933 , 76.84683 ,
         76.97183 ],
        ...,
        [40.815586, 41.753086, 42.034336, ..., 99.25308 , 99.06558 ,
         98.81558 ],
        [41.878086, 42.190586, 42.534336, ..., 99.22183 , 99.00308 ,
         98.78433 ],
        [41.971836, 42.565586, 42.565586, ..., 99.15933 , 98.97183 ,
         98.75308 ]],

       [[81.07036 , 81.38286 , 82.00786 , ..., 77.28911 , 77.32036 ,
         77.35161 ],
        [81.72661 , 82.10161 , 82.44536 , ..., 77.25786 , 77.28911 ,
         77.35161 ],
        [82.44536 , 82.47661 , 82.41411 , ..., 77.22661 , 77.28911 ,
         77.32036 ],
...
        [87.815216, 90.08084 , 90.846466, ..., 98.98709 , 99.01834 ,
         99.04959 ],
        [86.690216, 89.002716, 90.065216, ..., 98.95584 , 98.98709 ,
         99.01834 ],
        [85.08084 , 87.04959 , 89.096466, ..., 98.940216, 98.971466,
         99.002716]],

       [[82.16208 , 82.146454, 82.13083 , ..., 75.94333 , 76.16208 ,
         76.41208 ],
        [82.240204, 82.240204, 82.240204, ..., 75.97458 , 76.19333 ,
         76.427704],
        [82.28708 , 82.28708 , 82.302704, ..., 76.00583 , 76.240204,
         76.458954],
        ...,
        [87.677704, 89.44333 , 89.896454, ..., 98.927704, 98.990204,
         99.052704],
        [86.615204, 88.50583 , 89.34958 , ..., 98.91208 , 98.958954,
         99.021454],
        [85.00583 , 86.677704, 88.333954, ..., 98.896454, 98.94333 ,
         98.990204]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([              22.4,              22.45,               22.5,
                           22.55,               22.6,              22.65,
                            22.7,              22.75,               22.8,
                           22.85,
              ...
               47.14999999999999, 47.199999999999996,  47.24999999999999,
                            47.3, 47.349999999999994,  47.39999999999999,
              47.449999999999996,  47.49999999999999,              47.55,
              47.599999999999994],
             dtype='float64', name='latitude', length=505))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([   120.0, 120.0625,  120.125, 120.1875,   120.25, 120.3125,
               120.375, 120.4375,    120.5, 120.5625,
              ...
              149.4375,    149.5, 149.5625,  149.625, 149.6875,   149.75,
              149.8125,  149.875, 149.9375,    150.0],
             dtype='float64', name='longitude', length=481))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00', '2018-07-01 07:00:00',
               '2018-07-01 08:00:00', '2018-07-01 09:00:00',
               '2018-07-01 10:00:00', '2018-07-01 11:00:00',
               '2018-07-01 12:00:00', '2018-07-01 13:00:00',
               '2018-07-01 14:00:00', '2018-07-01 15:00:00',
               '2018-07-01 16:00:00', '2018-07-01 17:00:00',
               '2018-07-01 18:00:00', '2018-07-01 19:00:00',
               '2018-07-01 20:00:00', '2018-07-01 21:00:00'],
              dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0

　気温と湿度だけを内容とするDatasetオブジェクトが作成されました。

　関数 getgpv は、キーワード引数 verbose に True を与えて実行すると、関数が実際にアクセスした場所を表示します。上のスクリプトでは、True に設定しているので、アクセスしたインベントリも表示されました。

　関数 getgpv は、キーワード引数 to_netcdf にTrueを設定すると、GRIB2ファイルから読み出したデータをNetCDFファイルで保存します。そして、キーワード引数 from_netcdf にTrueを設定すると、これから読もうとするデータがすでにNetcdfファイルとして保存されていれば、そちらからデータを取得します。Netcdfファイルからの読み込み速度は、GRIB2ファイルからよりもずっと高速です。
Script 15 では、 to_netcdf を設定しているので、NetCDFファイルが作られています。そこで、今度は、キーワード引数 from_netcdf にTrueを設定して同じデータを読みなおしてみましょう。以下を実行してください。

In [16]:
# Script 16 #

grbdir = "jmadata/msm/2018/201807"
grbfile = "Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
grbpath = grbdir +"/"+ grbfile

elements = ["TMP","RH"]
ds = wx.getgpv(grbpath, elements, 
               ncdir="./nc", from_netcdf=True, verbose=True)
ds

reading from nc\Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.TMP_RH.nc
Out[16]:
<xarray.Dataset>
Dimensions:              (latitude: 505, longitude: 481, time: 16)
Coordinates:
  * latitude             (latitude) float64 22.4 22.45 22.5 ... 47.5 47.55 47.6
  * longitude            (longitude) float64 120.0 120.1 120.1 ... 149.9 150.0
  * time                 (time) datetime64[ns] 2018-07-01T06:00:00 ... 2018-0...
Data variables:
    TMP_1D5maboveground  (time, latitude, longitude) float32 302.1 ... 280.1
    RH_1D5maboveground   (time, latitude, longitude) float32 81.88 ... 98.99
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 505
    • longitude: 481
    • time: 16
    • latitude
      (latitude)
      float64
      22.4 22.45 22.5 ... 47.5 47.55 47.6
      units :
      degrees_north
      long_name :
      latitude
      array([22.4 , 22.45, 22.5 , ..., 47.5 , 47.55, 47.6 ])
    • longitude
      (longitude)
      float64
      120.0 120.1 120.1 ... 149.9 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([120.    , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150.    ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00 ... 2018-07-...
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      3
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      forecasts or accumulated (including analyses), reference date is fixed
      time_step_setting :
      auto
      time_step :
      3600.0
      array(['2018-07-01T06:00:00.000000000', '2018-07-01T07:00:00.000000000',
       '2018-07-01T08:00:00.000000000', '2018-07-01T09:00:00.000000000',
       '2018-07-01T10:00:00.000000000', '2018-07-01T11:00:00.000000000',
       '2018-07-01T12:00:00.000000000', '2018-07-01T13:00:00.000000000',
       '2018-07-01T14:00:00.000000000', '2018-07-01T15:00:00.000000000',
       '2018-07-01T16:00:00.000000000', '2018-07-01T17:00:00.000000000',
       '2018-07-01T18:00:00.000000000', '2018-07-01T19:00:00.000000000',
       '2018-07-01T20:00:00.000000000', '2018-07-01T21:00:00.000000000'],
      dtype='datetime64[ns]')
    • TMP_1D5maboveground
      (time, latitude, longitude)
      float32
      302.1 302.1 302.1 ... 280.2 280.1
      short_name :
      TMP_1D5maboveground
      long_name :
      Temperature
      level :
      1.5 m above ground
      units :
      K
      array([[[302.1476 , 302.1398 , 302.1398 , ..., 301.65543, 301.632  ,
         301.61636],
        [302.12418, 302.12418, 302.11636, ..., 301.6476 , 301.632  ,
         301.61636],
        [302.10855, 302.10074, 302.10074, ..., 301.6476 , 301.632  ,
         301.60855],
        ...,
        [300.1476 , 299.9523 , 299.9601 , ..., 279.55386, 279.4523 ,
         279.34293],
        [299.8898 , 299.7023 , 299.67105, ..., 279.55386, 279.43668,
         279.3195 ],
        [299.99136, 299.71793, 299.65543, ..., 279.55386, 279.43668,
         279.3195 ]],

       [[302.191  , 302.191  , 302.15195, ..., 301.6285 , 301.6207 ,
         301.60507],
        [302.15195, 302.13632, 302.11288, ..., 301.6285 , 301.6207 ,
         301.60507],
        [302.10507, 302.10507, 302.09726, ..., 301.6285 , 301.61288,
         301.59726],
...
        [289.45142, 289.04517, 288.86548, ..., 280.22485, 280.16235,
         280.09204],
        [289.7561 , 289.30298, 289.08423, ..., 280.24048, 280.17798,
         280.10767],
        [290.1858 , 289.78735, 289.38892, ..., 280.2561 , 280.1936 ,
         280.11548]],

       [[302.15546, 302.15546, 302.16327, ..., 301.53046, 301.48358,
         301.44452],
        [302.1789 , 302.17108, 302.17108, ..., 301.52264, 301.48358,
         301.4367 ],
        [302.19452, 302.1867 , 302.1867 , ..., 301.50702, 301.46796,
         301.4289 ],
        ...,
        [289.34296, 289.00702, 288.89764, ..., 280.21796, 280.16327,
         280.13202],
        [289.60858, 289.21014, 289.03046, ..., 280.22577, 280.1789 ,
         280.13202],
        [290.00702, 289.63983, 289.3039 , ..., 280.2492 , 280.19452,
         280.13983]]], dtype=float32)
    • RH_1D5maboveground
      (time, latitude, longitude)
      float32
      81.88 81.78 81.69 ... 98.94 98.99
      short_name :
      RH_1D5maboveground
      long_name :
      Relative Humidity
      level :
      1.5 m above ground
      units :
      percent
      array([[[81.87808 , 81.78433 , 81.69058 , ..., 76.59683 , 76.72183 ,
         76.87808 ],
        [81.94058 , 81.84683 , 81.81558 , ..., 76.62808 , 76.78433 ,
         76.94058 ],
        [81.94058 , 81.90933 , 81.90933 , ..., 76.65933 , 76.84683 ,
         76.97183 ],
        ...,
        [40.815586, 41.753086, 42.034336, ..., 99.25308 , 99.06558 ,
         98.81558 ],
        [41.878086, 42.190586, 42.534336, ..., 99.22183 , 99.00308 ,
         98.78433 ],
        [41.971836, 42.565586, 42.565586, ..., 99.15933 , 98.97183 ,
         98.75308 ]],

       [[81.07036 , 81.38286 , 82.00786 , ..., 77.28911 , 77.32036 ,
         77.35161 ],
        [81.72661 , 82.10161 , 82.44536 , ..., 77.25786 , 77.28911 ,
         77.35161 ],
        [82.44536 , 82.47661 , 82.41411 , ..., 77.22661 , 77.28911 ,
         77.32036 ],
...
        [87.815216, 90.08084 , 90.846466, ..., 98.98709 , 99.01834 ,
         99.04959 ],
        [86.690216, 89.002716, 90.065216, ..., 98.95584 , 98.98709 ,
         99.01834 ],
        [85.08084 , 87.04959 , 89.096466, ..., 98.940216, 98.971466,
         99.002716]],

       [[82.16208 , 82.146454, 82.13083 , ..., 75.94333 , 76.16208 ,
         76.41208 ],
        [82.240204, 82.240204, 82.240204, ..., 75.97458 , 76.19333 ,
         76.427704],
        [82.28708 , 82.28708 , 82.302704, ..., 76.00583 , 76.240204,
         76.458954],
        ...,
        [87.677704, 89.44333 , 89.896454, ..., 98.927704, 98.990204,
         99.052704],
        [86.615204, 88.50583 , 89.34958 , ..., 98.91208 , 98.958954,
         99.021454],
        [85.00583 , 86.677704, 88.333954, ..., 98.896454, 98.94333 ,
         98.990204]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([              22.4,              22.45,               22.5,
                           22.55,               22.6,              22.65,
                            22.7,              22.75,               22.8,
                           22.85,
              ...
               47.14999999999999, 47.199999999999996,  47.24999999999999,
                            47.3, 47.349999999999994,  47.39999999999999,
              47.449999999999996,  47.49999999999999,              47.55,
              47.599999999999994],
             dtype='float64', name='latitude', length=505))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([   120.0, 120.0625,  120.125, 120.1875,   120.25, 120.3125,
               120.375, 120.4375,    120.5, 120.5625,
              ...
              149.4375,    149.5, 149.5625,  149.625, 149.6875,   149.75,
              149.8125,  149.875, 149.9375,    150.0],
             dtype='float64', name='longitude', length=481))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00', '2018-07-01 07:00:00',
               '2018-07-01 08:00:00', '2018-07-01 09:00:00',
               '2018-07-01 10:00:00', '2018-07-01 11:00:00',
               '2018-07-01 12:00:00', '2018-07-01 13:00:00',
               '2018-07-01 14:00:00', '2018-07-01 15:00:00',
               '2018-07-01 16:00:00', '2018-07-01 17:00:00',
               '2018-07-01 18:00:00', '2018-07-01 19:00:00',
               '2018-07-01 20:00:00', '2018-07-01 21:00:00'],
              dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0

　あっというまにDeatasetオブジェクトが作られました。表示から、先ほどと読み取り先が違うことも分かります。

　関数 getgpv は、引数 grbpath にファイルのリストを与えると、これらをすべて処理し接合したDatasetオブジェクトを一つ作成します。これは、複数のGRIB2ファイルで構成されているGPVファイルを処理する場合に便利です。

　MSM-GPVのもそのようなプロダクトの一つです。MSM-GPVは予報期間が異なる3つのファイルで構成されているので、下のCellを実行して、3つのファイルを一度に読んで、時系列ブラフを描いてください。Script 13の実行結果と比較すると、予報期間が、39時間先まで伸びていつことがわかります。

In [17]:
# Script 17 #

grblist= ["jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin",
          "jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH16-33_grib2.bin",
          "jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin"]

elements = ["TMP","RH"]
ds = wx.getgpv(grblist, elements, 
               ncdir="./nc")

da = ds["TMP_1D5maboveground"]
da[:,200,180].plot.line(x="time",figsize=(15,5))
Out[17]:
[<matplotlib.lines.Line2D at 0x1f929644a90>]
No description has been provided for this image

　関数 getgpv は、キーワード引数 lalomima に、緯度経度範囲を与えると、Datasetオブジェクトにする緯度経度範囲を指定することができます。この機能は、推計気象分布や解析雨量など、格子間隔が狭いGPVプロダクトを多数積み重ねる時に便利です。
範囲の指定は、緯度と経度それぞれの上限と下限からなる４要素のリスト（タプルも可）を用いて、下のように行います。

lalomima = [ 緯度の最小値, 緯度の最大値, 経度の最小値, 経度の最大値 ]

　それでは、緯度経度の範囲を指定して、MSM-GPVデータの北海道の部分だけのDatasetオブジェクトを作ります。にします。以下のCellを実行して下さい。

In [18]:
# Script 18 #

area = [41.350, 45.533, 139.327, 145.998] #北海道の範囲

grblist= ["jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin",
          "jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH16-33_grib2.bin",
          "jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin"]

elements = ["TMP","RH"]
ds = wx.getgpv(grblist, elements, 
               ncdir="./nc", lalomima=area)

da = ds["TMP_1D5maboveground"]
da[0,:,:].plot()

ds
Out[18]:
<xarray.Dataset>
Dimensions:              (latitude: 85, longitude: 108, time: 40)
Coordinates:
  * latitude             (latitude) float64 41.35 41.4 41.45 ... 45.5 45.55
  * longitude            (longitude) float64 139.3 139.4 139.4 ... 145.9 146.0
  * time                 (time) datetime64[ns] 2018-07-01T06:00:00 ... 2018-0...
Data variables:
    TMP_1D5maboveground  (time, latitude, longitude) float32 291.1 ... 282.8
    RH_1D5maboveground   (time, latitude, longitude) float32 99.16 ... 99.08
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 85
    • longitude: 108
    • time: 40
    • latitude
      (latitude)
      float64
      41.35 41.4 41.45 ... 45.5 45.55
      units :
      degrees_north
      long_name :
      latitude
      array([41.35, 41.4 , 41.45, 41.5 , 41.55, 41.6 , 41.65, 41.7 , 41.75, 41.8 ,
       41.85, 41.9 , 41.95, 42.  , 42.05, 42.1 , 42.15, 42.2 , 42.25, 42.3 ,
       42.35, 42.4 , 42.45, 42.5 , 42.55, 42.6 , 42.65, 42.7 , 42.75, 42.8 ,
       42.85, 42.9 , 42.95, 43.  , 43.05, 43.1 , 43.15, 43.2 , 43.25, 43.3 ,
       43.35, 43.4 , 43.45, 43.5 , 43.55, 43.6 , 43.65, 43.7 , 43.75, 43.8 ,
       43.85, 43.9 , 43.95, 44.  , 44.05, 44.1 , 44.15, 44.2 , 44.25, 44.3 ,
       44.35, 44.4 , 44.45, 44.5 , 44.55, 44.6 , 44.65, 44.7 , 44.75, 44.8 ,
       44.85, 44.9 , 44.95, 45.  , 45.05, 45.1 , 45.15, 45.2 , 45.25, 45.3 ,
       45.35, 45.4 , 45.45, 45.5 , 45.55])
    • longitude
      (longitude)
      float64
      139.3 139.4 139.4 ... 145.9 146.0
      units :
      degrees_east
      long_name :
      longitude
      array([139.3125, 139.375 , 139.4375, 139.5   , 139.5625, 139.625 , 139.6875,
       139.75  , 139.8125, 139.875 , 139.9375, 140.    , 140.0625, 140.125 ,
       140.1875, 140.25  , 140.3125, 140.375 , 140.4375, 140.5   , 140.5625,
       140.625 , 140.6875, 140.75  , 140.8125, 140.875 , 140.9375, 141.    ,
       141.0625, 141.125 , 141.1875, 141.25  , 141.3125, 141.375 , 141.4375,
       141.5   , 141.5625, 141.625 , 141.6875, 141.75  , 141.8125, 141.875 ,
       141.9375, 142.    , 142.0625, 142.125 , 142.1875, 142.25  , 142.3125,
       142.375 , 142.4375, 142.5   , 142.5625, 142.625 , 142.6875, 142.75  ,
       142.8125, 142.875 , 142.9375, 143.    , 143.0625, 143.125 , 143.1875,
       143.25  , 143.3125, 143.375 , 143.4375, 143.5   , 143.5625, 143.625 ,
       143.6875, 143.75  , 143.8125, 143.875 , 143.9375, 144.    , 144.0625,
       144.125 , 144.1875, 144.25  , 144.3125, 144.375 , 144.4375, 144.5   ,
       144.5625, 144.625 , 144.6875, 144.75  , 144.8125, 144.875 , 144.9375,
       145.    , 145.0625, 145.125 , 145.1875, 145.25  , 145.3125, 145.375 ,
       145.4375, 145.5   , 145.5625, 145.625 , 145.6875, 145.75  , 145.8125,
       145.875 , 145.9375, 146.    ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00 ... 2018-07-...
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      3
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      forecasts or accumulated (including analyses), reference date is fixed
      time_step_setting :
      auto
      time_step :
      3600.0
      array(['2018-07-01T06:00:00.000000000', '2018-07-01T07:00:00.000000000',
       '2018-07-01T08:00:00.000000000', '2018-07-01T09:00:00.000000000',
       '2018-07-01T10:00:00.000000000', '2018-07-01T11:00:00.000000000',
       '2018-07-01T12:00:00.000000000', '2018-07-01T13:00:00.000000000',
       '2018-07-01T14:00:00.000000000', '2018-07-01T15:00:00.000000000',
       '2018-07-01T16:00:00.000000000', '2018-07-01T17:00:00.000000000',
       '2018-07-01T18:00:00.000000000', '2018-07-01T19:00:00.000000000',
       '2018-07-01T20:00:00.000000000', '2018-07-01T21:00:00.000000000',
       '2018-07-01T22:00:00.000000000', '2018-07-01T23:00:00.000000000',
       '2018-07-02T00:00:00.000000000', '2018-07-02T01:00:00.000000000',
       '2018-07-02T02:00:00.000000000', '2018-07-02T03:00:00.000000000',
       '2018-07-02T04:00:00.000000000', '2018-07-02T05:00:00.000000000',
       '2018-07-02T06:00:00.000000000', '2018-07-02T07:00:00.000000000',
       '2018-07-02T08:00:00.000000000', '2018-07-02T09:00:00.000000000',
       '2018-07-02T10:00:00.000000000', '2018-07-02T11:00:00.000000000',
       '2018-07-02T12:00:00.000000000', '2018-07-02T13:00:00.000000000',
       '2018-07-02T14:00:00.000000000', '2018-07-02T15:00:00.000000000',
       '2018-07-02T16:00:00.000000000', '2018-07-02T17:00:00.000000000',
       '2018-07-02T18:00:00.000000000', '2018-07-02T19:00:00.000000000',
       '2018-07-02T20:00:00.000000000', '2018-07-02T21:00:00.000000000'],
      dtype='datetime64[ns]')
    • TMP_1D5maboveground
      (time, latitude, longitude)
      float32
      291.1 291.2 291.3 ... 282.8 282.8
      short_name :
      TMP_1D5maboveground
      long_name :
      Temperature
      level :
      1.5 m above ground
      units :
      K
      array([[[291.1398 , 291.21793, 291.31168, ..., 286.8976 , 286.92105,
         286.9601 ],
        [291.007  , 291.0851 , 291.21793, ..., 286.99136, 287.0695 ,
         287.12418],
        [290.93668, 290.96793, 291.10074, ..., 287.0695 , 287.12418,
         287.16324],
        ...,
        [288.10855, 288.11636, 288.10855, ..., 283.6476 , 283.65543,
         283.7023 ],
        [288.05386, 288.0695 , 288.0695 , ..., 283.66324, 283.67105,
         283.7101 ],
        [287.97574, 288.007  , 288.0226 , ..., 283.67886, 283.68668,
         283.7101 ]],

       [[291.22226, 291.21445, 291.21445, ..., 286.91757, 287.19882,
         287.42538],
        [291.11288, 291.13632, 291.1832 , ..., 286.98007, 287.21445,
         287.48007],
        [291.01132, 291.02695, 291.10507, ..., 287.21445, 287.22226,
         287.441  ],
...
        [287.30627, 287.36096, 287.41565, ..., 282.8141 , 282.7594 ,
         282.7516 ],
        [287.1969 , 287.2516 , 287.30627, ..., 282.7672 , 282.72034,
         282.6891 ],
        [287.0797 , 287.1266 , 287.17346, ..., 282.71252, 282.67346,
         282.65002]],

       [[292.15015, 292.17358, 292.15796, ..., 287.02515, 287.04077,
         287.03296],
        [292.1189 , 292.15015, 292.15015, ..., 287.13452, 287.11108,
         287.1267 ],
        [292.0642 , 292.07202, 292.1189 , ..., 287.1814 , 287.17358,
         287.20483],
        ...,
        [287.70483, 287.77515, 287.82202, ..., 282.73608, 282.72046,
         282.73608],
        [287.61108, 287.67358, 287.72046, ..., 282.7439 , 282.73608,
         282.7517 ],
        [287.47046, 287.53296, 287.57202, ..., 282.75952, 282.75952,
         282.7517 ]]], dtype=float32)
    • RH_1D5maboveground
      (time, latitude, longitude)
      float32
      99.16 99.13 99.13 ... 99.1 99.08
      short_name :
      RH_1D5maboveground
      long_name :
      Relative Humidity
      level :
      1.5 m above ground
      units :
      percent
      array([[[99.15933 , 99.12808 , 99.12808 , ..., 99.94058 , 99.90933 ,
         99.87808 ],
        [99.19058 , 99.19058 , 99.15933 , ..., 99.90933 , 99.84683 ,
         99.78433 ],
        [99.25308 , 99.28433 , 99.22183 , ..., 99.90933 , 99.87808 ,
         99.81558 ],
        ...,
        [99.15933 , 99.15933 , 99.15933 , ..., 99.31558 , 99.31558 ,
         99.31558 ],
        [99.15933 , 99.15933 , 99.15933 , ..., 99.31558 , 99.31558 ,
         99.31558 ],
        [99.15933 , 99.15933 , 99.15933 , ..., 99.31558 , 99.31558 ,
         99.31558 ]],

       [[99.07036 , 99.07036 , 99.07036 , ..., 99.88286 , 99.72661 ,
         99.53911 ],
        [99.07036 , 99.07036 , 99.07036 , ..., 99.88286 , 99.72661 ,
         99.57036 ],
        [99.07036 , 99.07036 , 99.07036 , ..., 99.82036 , 99.82036 ,
         99.66411 ],
...
        [98.93364 , 98.918015, 98.90239 , ..., 99.293015, 99.136765,
         99.05864 ],
        [98.949265, 98.949265, 98.949265, ..., 99.105515, 99.08989 ,
         99.136765],
        [98.949265, 98.96489 , 98.980515, ..., 99.15239 , 99.199265,
         99.27739 ]],

       [[92.39242 , 92.28304 , 92.31429 , ..., 99.68929 , 99.53304 ,
         99.37679 ],
        [92.42367 , 92.31429 , 92.26742 , ..., 99.43929 , 99.25179 ,
         99.17367 ],
        [92.61117 , 92.51742 , 92.34554 , ..., 99.17367 , 99.09554 ,
         99.09554 ],
        ...,
        [98.81429 , 98.89242 , 98.95492 , ..., 99.09554 , 99.12679 ,
         99.15804 ],
        [98.90804 , 98.97054 , 99.01742 , ..., 99.11117 , 99.12679 ,
         99.15804 ],
        [98.95492 , 99.00179 , 99.03304 , ..., 99.09554 , 99.09554 ,
         99.07992 ]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([41.349999999999994,               41.4, 41.449999999999996,
               41.49999999999999,              41.55, 41.599999999999994,
                           41.65, 41.699999999999996,  41.74999999999999,
                            41.8, 41.849999999999994,               41.9,
              41.949999999999996,  41.99999999999999,              42.05,
              42.099999999999994,              42.15, 42.199999999999996,
               42.24999999999999,               42.3, 42.349999999999994,
               42.39999999999999, 42.449999999999996,  42.49999999999999,
                           42.55, 42.599999999999994,  42.64999999999999,
              42.699999999999996,  42.74999999999999,               42.8,
              42.849999999999994,  42.89999999999999, 42.949999999999996,
               42.99999999999999,              43.05, 43.099999999999994,
               43.14999999999999, 43.199999999999996,  43.24999999999999,
                            43.3, 43.349999999999994,  43.39999999999999,
              43.449999999999996,  43.49999999999999,              43.55,
              43.599999999999994,  43.64999999999999, 43.699999999999996,
               43.74999999999999,               43.8, 43.849999999999994,
               43.89999999999999, 43.949999999999996,  43.99999999999999,
                           44.05, 44.099999999999994,  44.14999999999999,
              44.199999999999996,  44.24999999999999,               44.3,
              44.349999999999994,  44.39999999999999, 44.449999999999996,
               44.49999999999999,              44.55, 44.599999999999994,
               44.64999999999999, 44.699999999999996,  44.74999999999999,
                            44.8, 44.849999999999994,  44.89999999999999,
              44.949999999999996,  44.99999999999999,              45.05,
              45.099999999999994,  45.14999999999999, 45.199999999999996,
               45.24999999999999,               45.3, 45.349999999999994,
               45.39999999999999, 45.449999999999996,  45.49999999999999,
                           45.55],
             dtype='float64', name='latitude'))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([139.3125,  139.375, 139.4375,    139.5, 139.5625,  139.625,
              139.6875,   139.75, 139.8125,  139.875,
              ...
              145.4375,    145.5, 145.5625,  145.625, 145.6875,   145.75,
              145.8125,  145.875, 145.9375,    146.0],
             dtype='float64', name='longitude', length=108))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00', '2018-07-01 07:00:00',
               '2018-07-01 08:00:00', '2018-07-01 09:00:00',
               '2018-07-01 10:00:00', '2018-07-01 11:00:00',
               '2018-07-01 12:00:00', '2018-07-01 13:00:00',
               '2018-07-01 14:00:00', '2018-07-01 15:00:00',
               '2018-07-01 16:00:00', '2018-07-01 17:00:00',
               '2018-07-01 18:00:00', '2018-07-01 19:00:00',
               '2018-07-01 20:00:00', '2018-07-01 21:00:00',
               '2018-07-01 22:00:00', '2018-07-01 23:00:00',
               '2018-07-02 00:00:00', '2018-07-02 01:00:00',
               '2018-07-02 02:00:00', '2018-07-02 03:00:00',
               '2018-07-02 04:00:00', '2018-07-02 05:00:00',
               '2018-07-02 06:00:00', '2018-07-02 07:00:00',
               '2018-07-02 08:00:00', '2018-07-02 09:00:00',
               '2018-07-02 10:00:00', '2018-07-02 11:00:00',
               '2018-07-02 12:00:00', '2018-07-02 13:00:00',
               '2018-07-02 14:00:00', '2018-07-02 15:00:00',
               '2018-07-02 16:00:00', '2018-07-02 17:00:00',
               '2018-07-02 18:00:00', '2018-07-02 19:00:00',
               '2018-07-02 20:00:00', '2018-07-02 21:00:00'],
              dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0
No description has been provided for this image

　なお、GRIB2ファイルのインベントリで使われている気象要素の略号をリストで得たいときは、ライブラリwxbcgribxの関数 getgrvarname を使用します。

In [19]:
# Script 19 #

grbpath= "jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
wx.getgrvarname(grbpath)
Out[19]:
['APCP',
 'DSWRF',
 'HCDC',
 'LCDC',
 'MCDC',
 'PRES',
 'PRMSL',
 'RH',
 'TCDC',
 'TMP',
 'UGRD',
 'VGRD']

wxbcgribXのライセンス¶

ライブラリ **wxbcgribX** は、気象ビジネス推進コンソーシアム人材育成ワーキンググルプ内勉強会「気象データ×IT勉強会」が著作権を持ち、以下の条件で使用を認めます。なお、この条件は、「MITライセンス」と呼ばれているものと同一です。

Copyright (c) 2022 気象データ×IT勉強会

以下に定める条件に従い、本ソフトウェアおよび関連文書のファイル（以下「ソフトウェア」）の複製を取得するすべての人に対し、ソフトウェアを無制限に扱うことを無償で許可します。これには、ソフトウェアの複製を使用、複写、変更、結合、掲載、頒布、サブライセンス、および/または販売する権利、およびソフトウェアを提供する相手に同じことを許可する権利も無制限に含まれます。

上記の著作権表示および本許諾表示を、ソフトウェアのすべての複製または重要な部分に記載するものとします。

ソフトウェアは「現状のまま」で、明示であるか暗黙であるかを問わず、何らの保証もなく提供されます。ここでいう保証とは、商品性、特定の目的への適合性、および権利非侵害についての保証も含みますが、それに限定されるものではありません。 作者または著作権者は、契約行為、不法行為、またはそれ以外であろうと、ソフトウェアに起因または関連し、あるいはソフトウェアの使用またはその他の扱いによって生じる一切の請求、損害、その他の義務について何らの責任も負わないものとします。

おわりに¶

　この課では、気象庁が様々なメッシュデータをGRIB2というフォーマットのファイルで配信していることを学習しました。そして、アメリカ大気海洋局の気候予測センターが公開するプログラムwgrib2 をPythonスクリプトで実行し、GRIB2ファイルの中身を見たりデータを取り出したりできることを学習しました。

　さらに、PythonにはXarrayと呼ばれるライブラリがあり、これを使うとメッシュ気象データが簡単に取り扱え（そうであ）ることを学びました。また、ライブラリ wxbcgribx を使うと、GRIB2ファイルの中のメッシュ気象データを簡単にDatasetオブジェクトやDataArrayオブジェクトにできることを学習しました。

　第２課では、Xarray.DataArrayオブジェクトメソッドを駆使した。メッシュ気象データの処理テクニックを詳しく学習します。

　ところで、はじめに　で、気象庁はメッシュ気象データのことを GPV データと呼んでいると記しましたが、両者は厳密には同一ではありません。GPVは、気象庁によれば、 GPV とは、Grid Point Value の略であり、縦の格子と横の格子の交わった点における値の集合で事物の分布を近似すること意図したデータです。一方メッシュデータとは、連続的に分布する事物を、縦の格子と横の格子で区切られる矩形の範囲で計量して一つの値とし、その集合により事物の分布を表現しようとする考え方のデータです。確かに、気象庁推計気象分布（気温）などは、矩形の中心点の値を矩形の代表値と考えても大きな違いはなく、GPVとメッシュデータとの違いをとやかく言う必要がない場合は少なく在りませんが、例えば、人口の分布など、メッシュデータとしては表現できるがGPVでは表現できない量もあることは、知っておいてください。

ANNEX¶

　せっかくですので、wxbcgribX の関数 getgpv を用いて、サンプルとして用意した気象庁のGPVプロダクトを概観してみましょう。使用するスクリプトの詳細については、第2課で学習します。

In [20]:
# Script A1 #

import wxbcgribx as wx
from pathlib import Path

data_root = Path("./jmadata")
yyyymmdd = "20180701"

A 全球数値予報モデルGPV　(GSM-GPV)¶

In [21]:
# Script A2 #

product = "gsm"
grbdir = data_root/product/yyyymmdd[0:4]/yyyymmdd[0:6]
frange = ["FD0000-0312","FD0315-0800","FD0803-1100"]
hh = "12"
grblist= [Path(f"{grbdir}/Z__C_RJTD_{yyyymmdd}{hh}0000_GSM_GPV_Rjp_Lsurf_{oo}_grib2.bin")
         for oo in frange]
print(wx.getgrvarname(grblist[0]))
grblist

['APCP', 'DSWRF', 'HCDC', 'LCDC', 'MCDC', 'PRES', 'PRMSL', 'RH', 'TCDC', 'TMP', 'UGRD', 'VGRD']
Out[21]:
[WindowsPath('jmadata/gsm/2018/201807/Z__C_RJTD_20180701120000_GSM_GPV_Rjp_Lsurf_FD0000-0312_grib2.bin'),
 WindowsPath('jmadata/gsm/2018/201807/Z__C_RJTD_20180701120000_GSM_GPV_Rjp_Lsurf_FD0315-0800_grib2.bin'),
 WindowsPath('jmadata/gsm/2018/201807/Z__C_RJTD_20180701120000_GSM_GPV_Rjp_Lsurf_FD0803-1100_grib2.bin')]
In [22]:
# Script A3 #

elements = ["TMP"]
ds = wx.getgpv(grblist,elements,to_netcdf=False)
ds
Out[22]:
<xarray.Dataset>
Dimensions:            (latitude: 151, longitude: 121, time: 145)
Coordinates:
  * latitude           (latitude) float64 20.0 20.2 20.4 20.6 ... 49.6 49.8 50.0
  * longitude          (longitude) float64 120.0 120.2 120.5 ... 149.8 150.0
  * time               (time) datetime64[ns] 2018-07-01T12:00:00 ... 2018-07-...
Data variables:
    TMP_2maboveground  (time, latitude, longitude) float32 301.0 301.0 ... 280.8
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 151
    • longitude: 121
    • time: 145
    • latitude
      (latitude)
      float64
      20.0 20.2 20.4 ... 49.6 49.8 50.0
      units :
      degrees_north
      long_name :
      latitude
      array([20. , 20.2, 20.4, 20.6, 20.8, 21. , 21.2, 21.4, 21.6, 21.8, 22. , 22.2,
       22.4, 22.6, 22.8, 23. , 23.2, 23.4, 23.6, 23.8, 24. , 24.2, 24.4, 24.6,
       24.8, 25. , 25.2, 25.4, 25.6, 25.8, 26. , 26.2, 26.4, 26.6, 26.8, 27. ,
       27.2, 27.4, 27.6, 27.8, 28. , 28.2, 28.4, 28.6, 28.8, 29. , 29.2, 29.4,
       29.6, 29.8, 30. , 30.2, 30.4, 30.6, 30.8, 31. , 31.2, 31.4, 31.6, 31.8,
       32. , 32.2, 32.4, 32.6, 32.8, 33. , 33.2, 33.4, 33.6, 33.8, 34. , 34.2,
       34.4, 34.6, 34.8, 35. , 35.2, 35.4, 35.6, 35.8, 36. , 36.2, 36.4, 36.6,
       36.8, 37. , 37.2, 37.4, 37.6, 37.8, 38. , 38.2, 38.4, 38.6, 38.8, 39. ,
       39.2, 39.4, 39.6, 39.8, 40. , 40.2, 40.4, 40.6, 40.8, 41. , 41.2, 41.4,
       41.6, 41.8, 42. , 42.2, 42.4, 42.6, 42.8, 43. , 43.2, 43.4, 43.6, 43.8,
       44. , 44.2, 44.4, 44.6, 44.8, 45. , 45.2, 45.4, 45.6, 45.8, 46. , 46.2,
       46.4, 46.6, 46.8, 47. , 47.2, 47.4, 47.6, 47.8, 48. , 48.2, 48.4, 48.6,
       48.8, 49. , 49.2, 49.4, 49.6, 49.8, 50. ])
    • longitude
      (longitude)
      float64
      120.0 120.2 120.5 ... 149.8 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([120.  , 120.25, 120.5 , 120.75, 121.  , 121.25, 121.5 , 121.75, 122.  ,
       122.25, 122.5 , 122.75, 123.  , 123.25, 123.5 , 123.75, 124.  , 124.25,
       124.5 , 124.75, 125.  , 125.25, 125.5 , 125.75, 126.  , 126.25, 126.5 ,
       126.75, 127.  , 127.25, 127.5 , 127.75, 128.  , 128.25, 128.5 , 128.75,
       129.  , 129.25, 129.5 , 129.75, 130.  , 130.25, 130.5 , 130.75, 131.  ,
       131.25, 131.5 , 131.75, 132.  , 132.25, 132.5 , 132.75, 133.  , 133.25,
       133.5 , 133.75, 134.  , 134.25, 134.5 , 134.75, 135.  , 135.25, 135.5 ,
       135.75, 136.  , 136.25, 136.5 , 136.75, 137.  , 137.25, 137.5 , 137.75,
       138.  , 138.25, 138.5 , 138.75, 139.  , 139.25, 139.5 , 139.75, 140.  ,
       140.25, 140.5 , 140.75, 141.  , 141.25, 141.5 , 141.75, 142.  , 142.25,
       142.5 , 142.75, 143.  , 143.25, 143.5 , 143.75, 144.  , 144.25, 144.5 ,
       144.75, 145.  , 145.25, 145.5 , 145.75, 146.  , 146.25, 146.5 , 146.75,
       147.  , 147.25, 147.5 , 147.75, 148.  , 148.25, 148.5 , 148.75, 149.  ,
       149.25, 149.5 , 149.75, 150.  ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T12:00:00 ... 2018-07-...
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530446400.0
      reference_time_type :
      3
      reference_date :
      2018.07.01 12:00:00 UTC
      reference_time_description :
      forecasts or accumulated (including analyses), reference date is fixed
      time_step_setting :
      auto
      time_step :
      3600.0
      array(['2018-07-01T12:00:00.000000000', '2018-07-01T13:00:00.000000000',
       '2018-07-01T14:00:00.000000000', '2018-07-01T15:00:00.000000000',
       '2018-07-01T16:00:00.000000000', '2018-07-01T17:00:00.000000000',
       '2018-07-01T18:00:00.000000000', '2018-07-01T19:00:00.000000000',
       '2018-07-01T20:00:00.000000000', '2018-07-01T21:00:00.000000000',
       '2018-07-01T22:00:00.000000000', '2018-07-01T23:00:00.000000000',
       '2018-07-02T00:00:00.000000000', '2018-07-02T01:00:00.000000000',
       '2018-07-02T02:00:00.000000000', '2018-07-02T03:00:00.000000000',
       '2018-07-02T04:00:00.000000000', '2018-07-02T05:00:00.000000000',
       '2018-07-02T06:00:00.000000000', '2018-07-02T07:00:00.000000000',
       '2018-07-02T08:00:00.000000000', '2018-07-02T09:00:00.000000000',
       '2018-07-02T10:00:00.000000000', '2018-07-02T11:00:00.000000000',
       '2018-07-02T12:00:00.000000000', '2018-07-02T13:00:00.000000000',
       '2018-07-02T14:00:00.000000000', '2018-07-02T15:00:00.000000000',
       '2018-07-02T16:00:00.000000000', '2018-07-02T17:00:00.000000000',
       '2018-07-02T18:00:00.000000000', '2018-07-02T19:00:00.000000000',
       '2018-07-02T20:00:00.000000000', '2018-07-02T21:00:00.000000000',
       '2018-07-02T22:00:00.000000000', '2018-07-02T23:00:00.000000000',
       '2018-07-03T00:00:00.000000000', '2018-07-03T01:00:00.000000000',
       '2018-07-03T02:00:00.000000000', '2018-07-03T03:00:00.000000000',
       '2018-07-03T04:00:00.000000000', '2018-07-03T05:00:00.000000000',
       '2018-07-03T06:00:00.000000000', '2018-07-03T07:00:00.000000000',
       '2018-07-03T08:00:00.000000000', '2018-07-03T09:00:00.000000000',
       '2018-07-03T10:00:00.000000000', '2018-07-03T11:00:00.000000000',
       '2018-07-03T12:00:00.000000000', '2018-07-03T13:00:00.000000000',
       '2018-07-03T14:00:00.000000000', '2018-07-03T15:00:00.000000000',
       '2018-07-03T16:00:00.000000000', '2018-07-03T17:00:00.000000000',
       '2018-07-03T18:00:00.000000000', '2018-07-03T19:00:00.000000000',
       '2018-07-03T20:00:00.000000000', '2018-07-03T21:00:00.000000000',
       '2018-07-03T22:00:00.000000000', '2018-07-03T23:00:00.000000000',
       '2018-07-04T00:00:00.000000000', '2018-07-04T01:00:00.000000000',
       '2018-07-04T02:00:00.000000000', '2018-07-04T03:00:00.000000000',
       '2018-07-04T04:00:00.000000000', '2018-07-04T05:00:00.000000000',
       '2018-07-04T06:00:00.000000000', '2018-07-04T07:00:00.000000000',
       '2018-07-04T08:00:00.000000000', '2018-07-04T09:00:00.000000000',
       '2018-07-04T10:00:00.000000000', '2018-07-04T11:00:00.000000000',
       '2018-07-04T12:00:00.000000000', '2018-07-04T13:00:00.000000000',
       '2018-07-04T14:00:00.000000000', '2018-07-04T15:00:00.000000000',
       '2018-07-04T16:00:00.000000000', '2018-07-04T17:00:00.000000000',
       '2018-07-04T18:00:00.000000000', '2018-07-04T19:00:00.000000000',
       '2018-07-04T20:00:00.000000000', '2018-07-04T21:00:00.000000000',
       '2018-07-04T22:00:00.000000000', '2018-07-04T23:00:00.000000000',
       '2018-07-05T00:00:00.000000000', '2018-07-05T03:00:00.000000000',
       '2018-07-05T06:00:00.000000000', '2018-07-05T09:00:00.000000000',
       '2018-07-05T12:00:00.000000000', '2018-07-05T15:00:00.000000000',
       '2018-07-05T18:00:00.000000000', '2018-07-05T21:00:00.000000000',
       '2018-07-06T00:00:00.000000000', '2018-07-06T03:00:00.000000000',
       '2018-07-06T06:00:00.000000000', '2018-07-06T09:00:00.000000000',
       '2018-07-06T12:00:00.000000000', '2018-07-06T15:00:00.000000000',
       '2018-07-06T18:00:00.000000000', '2018-07-06T21:00:00.000000000',
       '2018-07-07T00:00:00.000000000', '2018-07-07T03:00:00.000000000',
       '2018-07-07T06:00:00.000000000', '2018-07-07T09:00:00.000000000',
       '2018-07-07T12:00:00.000000000', '2018-07-07T15:00:00.000000000',
       '2018-07-07T18:00:00.000000000', '2018-07-07T21:00:00.000000000',
       '2018-07-08T00:00:00.000000000', '2018-07-08T03:00:00.000000000',
       '2018-07-08T06:00:00.000000000', '2018-07-08T09:00:00.000000000',
       '2018-07-08T12:00:00.000000000', '2018-07-08T15:00:00.000000000',
       '2018-07-08T18:00:00.000000000', '2018-07-08T21:00:00.000000000',
       '2018-07-09T00:00:00.000000000', '2018-07-09T03:00:00.000000000',
       '2018-07-09T06:00:00.000000000', '2018-07-09T09:00:00.000000000',
       '2018-07-09T12:00:00.000000000', '2018-07-09T15:00:00.000000000',
       '2018-07-09T18:00:00.000000000', '2018-07-09T21:00:00.000000000',
       '2018-07-10T00:00:00.000000000', '2018-07-10T03:00:00.000000000',
       '2018-07-10T06:00:00.000000000', '2018-07-10T09:00:00.000000000',
       '2018-07-10T12:00:00.000000000', '2018-07-10T15:00:00.000000000',
       '2018-07-10T18:00:00.000000000', '2018-07-10T21:00:00.000000000',
       '2018-07-11T00:00:00.000000000', '2018-07-11T03:00:00.000000000',
       '2018-07-11T06:00:00.000000000', '2018-07-11T09:00:00.000000000',
       '2018-07-11T12:00:00.000000000', '2018-07-11T15:00:00.000000000',
       '2018-07-11T18:00:00.000000000', '2018-07-11T21:00:00.000000000',
       '2018-07-12T00:00:00.000000000', '2018-07-12T03:00:00.000000000',
       '2018-07-12T06:00:00.000000000', '2018-07-12T09:00:00.000000000',
       '2018-07-12T12:00:00.000000000'], dtype='datetime64[ns]')
    • TMP_2maboveground
      (time, latitude, longitude)
      float32
      301.0 301.0 301.0 ... 280.8 280.8
      short_name :
      TMP_2maboveground
      long_name :
      Temperature
      level :
      2 m above ground
      units :
      K
      array([[[301.0494 , 301.0416 , 301.01035, ..., 301.50253, 301.47128,
         301.44785],
        [301.02597, 301.01035, 300.9791 , ..., 301.46347, 301.43222,
         301.40878],
        [301.07285, 301.0416 , 300.9791 , ..., 301.45566, 301.4244 ,
         301.39316],
        ...,
        [296.1041 , 296.05722, 295.93222, ..., 280.13535, 279.95566,
         279.82285],
        [295.50253, 295.40878, 295.15097, ..., 280.06503, 279.87753,
         279.74472],
        [295.13535, 294.82285, 294.36972, ..., 280.03378, 279.86972,
         279.7369 ]],

       [[301.1366 , 301.09753, 301.0194 , ..., 301.37878, 301.36316,
         301.35535],
        [301.12097, 301.08972, 301.02722, ..., 301.35535, 301.3319 ,
         301.3241 ],
        [301.15222, 301.12097, 301.04285, ..., 301.36316, 301.33972,
         301.3241 ],
...
        [298.75604, 298.24042, 297.81854, ..., 280.70917, 280.75604,
         280.74042],
        [297.90448, 297.38104, 296.8498 , ..., 280.77167, 280.82635,
         280.84198],
        [297.4748 , 296.90448, 296.27948, ..., 280.83417, 280.9123 ,
         280.96698]],

       [[301.19238, 301.22363, 301.28613, ..., 300.81738, 300.7627 ,
         300.73145],
        [301.2002 , 301.23145, 301.29395, ..., 300.80957, 300.67676,
         300.62207],
        [301.23926, 301.25488, 301.30957, ..., 300.72363, 300.59082,
         300.52832],
        ...,
        [295.48145, 294.9502 , 294.5752 , ..., 280.59082, 280.583  ,
         280.55957],
        [295.35645, 294.62988, 293.84863, ..., 280.65332, 280.66895,
         280.65332],
        [295.1377 , 294.35645, 293.36426, ..., 280.74707, 280.78613,
         280.7705 ]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([              20.0,               20.2,               20.4,
              20.599999999999998, 20.799999999999997,               21.0,
                            21.2,               21.4, 21.599999999999998,
              21.799999999999997,
              ...
                            48.2,               48.4,               48.6,
                            48.8,               49.0,               49.2,
                            49.4,               49.6,               49.8,
                            50.0],
             dtype='float64', name='latitude', length=151))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([ 120.0, 120.25,  120.5, 120.75,  121.0, 121.25,  121.5, 121.75,
               122.0, 122.25,
              ...
              147.75,  148.0, 148.25,  148.5, 148.75,  149.0, 149.25,  149.5,
              149.75,  150.0],
             dtype='float64', name='longitude', length=121))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 12:00:00', '2018-07-01 13:00:00',
               '2018-07-01 14:00:00', '2018-07-01 15:00:00',
               '2018-07-01 16:00:00', '2018-07-01 17:00:00',
               '2018-07-01 18:00:00', '2018-07-01 19:00:00',
               '2018-07-01 20:00:00', '2018-07-01 21:00:00',
               ...
               '2018-07-11 09:00:00', '2018-07-11 12:00:00',
               '2018-07-11 15:00:00', '2018-07-11 18:00:00',
               '2018-07-11 21:00:00', '2018-07-12 00:00:00',
               '2018-07-12 03:00:00', '2018-07-12 06:00:00',
               '2018-07-12 09:00:00', '2018-07-12 12:00:00'],
              dtype='datetime64[ns]', name='time', length=145, freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0
In [23]:
# Script A4 #

da = ds["TMP_2maboveground"]
da[-1,:,:].plot(cmap='RdYlGn_r',figsize=(5,5))
Out[23]:
<matplotlib.collections.QuadMesh at 0x1f929f73dc0>
No description has been provided for this image
In [24]:
# Script A5 #

da.sel(latitude=35.6916,longitude=139.7516, method='nearest').plot.line(x="time",figsize=(12,3))
Out[24]:
[<matplotlib.lines.Line2D at 0x1f92a00cd30>]
No description has been provided for this image

B 局地数値予報モデルGPV　(LFM-GPV)¶

In [25]:
# Script A6 #

product = "lfm"
grbdir = data_root/product/yyyymmdd[0:4]/yyyymmdd[0:6]
dtr = wx.dt_range("202001010000", 19, minutes=30)
frange = [oo.strftime('FH%H%M') for oo in dtr]
hh = "06"
grblist= [Path(f"{grbdir}/Z__C_RJTD_{yyyymmdd}{hh}0000_LFM_GPV_Rjp_Lsurf_{oo}_grib2.bin")
         for oo in frange]
print(wx.getvarname(grblist[0]))
grblist

['HCDC', 'LCDC', 'MCDC', 'PRES', 'PRMSL', 'RH', 'TCDC', 'TMP', 'UGRD', 'VGRD']
Out[25]:
[WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0000_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0030_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0100_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0130_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0200_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0230_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0300_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0330_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0400_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0430_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0500_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0530_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0600_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0630_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0700_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0730_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0800_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0830_grib2.bin'),
 WindowsPath('jmadata/lfm/2018/201807/Z__C_RJTD_20180701060000_LFM_GPV_Rjp_Lsurf_FH0900_grib2.bin')]
In [26]:
# Script A7 #

elements = ["TMP"]
ds = wx.getgpv(grblist[0],elements,to_netcdf=False) #最初の１つしか読んでいません
ds
Out[26]:
<xarray.Dataset>
Dimensions:              (latitude: 1261, longitude: 1201, time: 1)
Coordinates:
  * latitude             (latitude) float64 22.4 22.42 22.44 ... 47.58 47.6
  * longitude            (longitude) float64 120.0 120.0 120.1 ... 150.0 150.0
  * time                 (time) datetime64[ns] 2018-07-01T06:00:00
Data variables:
    TMP_1D5maboveground  (time, latitude, longitude) float32 302.1 302.1 ... nan
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 1261
    • longitude: 1201
    • time: 1
    • latitude
      (latitude)
      float64
      22.4 22.42 22.44 ... 47.58 47.6
      units :
      degrees_north
      long_name :
      latitude
      array([22.4 , 22.42, 22.44, ..., 47.56, 47.58, 47.6 ])
    • longitude
      (longitude)
      float64
      120.0 120.0 120.1 ... 150.0 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([120.   , 120.025, 120.05 , ..., 149.95 , 149.975, 150.   ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      1
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      analyses, reference date is fixed
      time_step_setting :
      auto
      time_step :
      0.0
      array(['2018-07-01T06:00:00.000000000'], dtype='datetime64[ns]')
    • TMP_1D5maboveground
      (time, latitude, longitude)
      float32
      302.1 302.1 302.1 ... nan nan nan
      short_name :
      TMP_1D5maboveground
      long_name :
      Temperature
      level :
      1.5 m above ground
      units :
      K
      array([[[302.13507, 302.12726, 302.11945, ...,       nan,       nan,
               nan],
        [302.12726, 302.11945, 302.11163, ...,       nan,       nan,
               nan],
        [302.12726, 302.11945, 302.11163, ...,       nan,       nan,
               nan],
        ...,
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([              22.4, 22.419999999999998, 22.439999999999998,
              22.459999999999997, 22.479999999999997, 22.499999999999996,
                           22.52,              22.54,              22.56,
                           22.58,
              ...
              47.419999999999995,              47.44, 47.459999999999994,
                           47.48,  47.49999999999999, 47.519999999999996,
               47.53999999999999, 47.559999999999995,  47.57999999999999,
              47.599999999999994],
             dtype='float64', name='latitude', length=1261))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([             120.0,            120.025, 120.05000000000001,
              120.07500000000002, 120.10000000000002, 120.12500000000003,
              120.15000000000003, 120.17500000000004, 120.20000000000005,
              120.22500000000005,
              ...
              149.77500000000677, 149.80000000000678, 149.82500000000678,
               149.8500000000068,  149.8750000000068,  149.9000000000068,
               149.9250000000068,  149.9500000000068, 149.97500000000682,
              150.00000000000682],
             dtype='float64', name='longitude', length=1201))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00'], dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0
In [27]:
# Script A8 #

da = ds["TMP_1D5maboveground"]
da[0,:,:].plot(cmap='RdYlGn_r',figsize=(5,5))
Out[27]:
<matplotlib.collections.QuadMesh at 0x1f92a0a29b0>
No description has been provided for this image
In [28]:
# Script A9 #

# 全域をオブジェクト化するとサイズが巨大になるので領域を限定して全ファイルを読みます
area = (35.6916,35.6916,139.7516,139.7516)  #アメダス「東京」の緯度経度
ds = wx.getgpv(grblist,elements,to_netcdf=False,lalomima=area)

da = ds["TMP_1D5maboveground"]
da.sel(latitude=35.6916,longitude=139.7516, method='nearest').plot.line(x="time",figsize=(12,3))
Out[28]:
[<matplotlib.lines.Line2D at 0x1f92a79dff0>]
No description has been provided for this image

C メソ数値予報モデルGPV　(MSM-GPV)¶

In [29]:
# Script A10 #

product = "msm"
grbdir = data_root/product/yyyymmdd[0:4]/yyyymmdd[0:6]
frange = ["FH00-15","FH16-33","FH34-39"]
hh = "06"
grblist= [Path(f"{grbdir}/Z__C_RJTD_{yyyymmdd}{hh}0000_MSM_GPV_Rjp_Lsurf_{oo}_grib2.bin")
         for oo in frange]
print(wx.getgrvarname(grblist[0]))
grblist

['APCP', 'DSWRF', 'HCDC', 'LCDC', 'MCDC', 'PRES', 'PRMSL', 'RH', 'TCDC', 'TMP', 'UGRD', 'VGRD']
Out[29]:
[WindowsPath('jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin'),
 WindowsPath('jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH16-33_grib2.bin'),
 WindowsPath('jmadata/msm/2018/201807/Z__C_RJTD_20180701060000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin')]
In [30]:
# Script A11 #

elements = ["TMP","RH","APCP"]
ds = wx.getgpv(grblist,elements,to_netcdf=False)
ds
Out[30]:
<xarray.Dataset>
Dimensions:              (latitude: 505, longitude: 481, time: 40)
Coordinates:
  * latitude             (latitude) float64 22.4 22.45 22.5 ... 47.5 47.55 47.6
  * longitude            (longitude) float64 120.0 120.1 120.1 ... 149.9 150.0
  * time                 (time) datetime64[ns] 2018-07-01T06:00:00 ... 2018-0...
Data variables:
    TMP_1D5maboveground  (time, latitude, longitude) float32 302.1 ... 279.5
    RH_1D5maboveground   (time, latitude, longitude) float32 81.88 ... 98.3
    APCP_surface         (time, latitude, longitude) float32 nan nan ... 0.0 0.0
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 505
    • longitude: 481
    • time: 40
    • latitude
      (latitude)
      float64
      22.4 22.45 22.5 ... 47.5 47.55 47.6
      units :
      degrees_north
      long_name :
      latitude
      array([22.4 , 22.45, 22.5 , ..., 47.5 , 47.55, 47.6 ])
    • longitude
      (longitude)
      float64
      120.0 120.1 120.1 ... 149.9 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([120.    , 120.0625, 120.125 , ..., 149.875 , 149.9375, 150.    ])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00 ... 2018-07-...
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      3
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      forecasts or accumulated (including analyses), reference date is fixed
      time_step_setting :
      auto
      time_step :
      3600.0
      array(['2018-07-01T06:00:00.000000000', '2018-07-01T07:00:00.000000000',
       '2018-07-01T08:00:00.000000000', '2018-07-01T09:00:00.000000000',
       '2018-07-01T10:00:00.000000000', '2018-07-01T11:00:00.000000000',
       '2018-07-01T12:00:00.000000000', '2018-07-01T13:00:00.000000000',
       '2018-07-01T14:00:00.000000000', '2018-07-01T15:00:00.000000000',
       '2018-07-01T16:00:00.000000000', '2018-07-01T17:00:00.000000000',
       '2018-07-01T18:00:00.000000000', '2018-07-01T19:00:00.000000000',
       '2018-07-01T20:00:00.000000000', '2018-07-01T21:00:00.000000000',
       '2018-07-01T22:00:00.000000000', '2018-07-01T23:00:00.000000000',
       '2018-07-02T00:00:00.000000000', '2018-07-02T01:00:00.000000000',
       '2018-07-02T02:00:00.000000000', '2018-07-02T03:00:00.000000000',
       '2018-07-02T04:00:00.000000000', '2018-07-02T05:00:00.000000000',
       '2018-07-02T06:00:00.000000000', '2018-07-02T07:00:00.000000000',
       '2018-07-02T08:00:00.000000000', '2018-07-02T09:00:00.000000000',
       '2018-07-02T10:00:00.000000000', '2018-07-02T11:00:00.000000000',
       '2018-07-02T12:00:00.000000000', '2018-07-02T13:00:00.000000000',
       '2018-07-02T14:00:00.000000000', '2018-07-02T15:00:00.000000000',
       '2018-07-02T16:00:00.000000000', '2018-07-02T17:00:00.000000000',
       '2018-07-02T18:00:00.000000000', '2018-07-02T19:00:00.000000000',
       '2018-07-02T20:00:00.000000000', '2018-07-02T21:00:00.000000000'],
      dtype='datetime64[ns]')
    • TMP_1D5maboveground
      (time, latitude, longitude)
      float32
      302.1 302.1 302.1 ... 279.5 279.5
      short_name :
      TMP_1D5maboveground
      long_name :
      Temperature
      level :
      1.5 m above ground
      units :
      K
      array([[[302.1476 , 302.1398 , 302.1398 , ..., 301.65543, 301.632  ,
         301.61636],
        [302.12418, 302.12418, 302.11636, ..., 301.6476 , 301.632  ,
         301.61636],
        [302.10855, 302.10074, 302.10074, ..., 301.6476 , 301.632  ,
         301.60855],
        ...,
        [300.1476 , 299.9523 , 299.9601 , ..., 279.55386, 279.4523 ,
         279.34293],
        [299.8898 , 299.7023 , 299.67105, ..., 279.55386, 279.43668,
         279.3195 ],
        [299.99136, 299.71793, 299.65543, ..., 279.55386, 279.43668,
         279.3195 ]],

       [[302.191  , 302.191  , 302.15195, ..., 301.6285 , 301.6207 ,
         301.60507],
        [302.15195, 302.13632, 302.11288, ..., 301.6285 , 301.6207 ,
         301.60507],
        [302.10507, 302.10507, 302.09726, ..., 301.6285 , 301.61288,
         301.59726],
...
        [289.3219 , 289.0797 , 289.0172 , ..., 279.5641 , 279.5172 ,
         279.47034],
        [289.29065, 288.93127, 288.7672 , ..., 279.58752, 279.54065,
         279.49377],
        [289.41565, 289.03284, 288.7672 , ..., 279.60315, 279.55627,
         279.5172 ]],

       [[302.15796, 302.15796, 302.15015, ..., 301.33765, 301.3064 ,
         301.27515],
        [302.15796, 302.15796, 302.15796, ..., 301.34546, 301.3142 ,
         301.28296],
        [302.15796, 302.15796, 302.15015, ..., 301.34546, 301.32202,
         301.28296],
        ...,
        [289.1267 , 288.92358, 288.90796, ..., 279.51733, 279.47046,
         279.4314 ],
        [289.08765, 288.77515, 288.65796, ..., 279.53296, 279.48608,
         279.44702],
        [289.22827, 288.88452, 288.67358, ..., 279.5564 , 279.50952,
         279.46265]]], dtype=float32)
    • RH_1D5maboveground
      (time, latitude, longitude)
      float32
      81.88 81.78 81.69 ... 98.31 98.3
      short_name :
      RH_1D5maboveground
      long_name :
      Relative Humidity
      level :
      1.5 m above ground
      units :
      percent
      array([[[81.87808 , 81.78433 , 81.69058 , ..., 76.59683 , 76.72183 ,
         76.87808 ],
        [81.94058 , 81.84683 , 81.81558 , ..., 76.62808 , 76.78433 ,
         76.94058 ],
        [81.94058 , 81.90933 , 81.90933 , ..., 76.65933 , 76.84683 ,
         76.97183 ],
        ...,
        [40.815586, 41.753086, 42.034336, ..., 99.25308 , 99.06558 ,
         98.81558 ],
        [41.878086, 42.190586, 42.534336, ..., 99.22183 , 99.00308 ,
         98.78433 ],
        [41.971836, 42.565586, 42.565586, ..., 99.15933 , 98.97183 ,
         98.75308 ]],

       [[81.07036 , 81.38286 , 82.00786 , ..., 77.28911 , 77.32036 ,
         77.35161 ],
        [81.72661 , 82.10161 , 82.44536 , ..., 77.25786 , 77.28911 ,
         77.35161 ],
        [82.44536 , 82.47661 , 82.41411 , ..., 77.22661 , 77.28911 ,
         77.32036 ],
...
        [83.77739 , 85.02739 , 85.543015, ..., 98.699265, 98.71489 ,
         98.730515],
        [83.30864 , 84.90239 , 85.824265, ..., 98.62114 , 98.636765,
         98.636765],
        [82.230515, 83.918015, 85.199265, ..., 98.52739 , 98.543015,
         98.55864 ]],

       [[79.65804 , 79.54867 , 79.50179 , ..., 76.78304 , 76.86117 ,
         76.95492 ],
        [80.32992 , 80.06429 , 79.90804 , ..., 76.70492 , 76.78304 ,
         76.87679 ],
        [81.01742 , 80.70492 , 80.45492 , ..., 76.64242 , 76.72054 ,
         76.81429 ],
        ...,
        [86.40804 , 87.43929 , 87.64242 , ..., 98.57992 , 98.57992 ,
         98.54867 ],
        [86.07992 , 87.40804 , 88.01742 , ..., 98.48617 , 98.43929 ,
         98.39242 ],
        [84.97054 , 86.42367 , 87.42367 , ..., 98.36117 , 98.31429 ,
         98.29867 ]]], dtype=float32)
    • APCP_surface
      (time, latitude, longitude)
      float32
      nan nan nan nan ... 0.0 0.0 0.0 0.0
      short_name :
      APCP_surface
      long_name :
      Total Precipitation
      level :
      surface
      units :
      kg/m^2
      array([[[     nan,      nan,      nan, ...,      nan,      nan,
              nan],
        [     nan,      nan,      nan, ...,      nan,      nan,
              nan],
        [     nan,      nan,      nan, ...,      nan,      nan,
              nan],
        ...,
        [     nan,      nan,      nan, ...,      nan,      nan,
              nan],
        [     nan,      nan,      nan, ...,      nan,      nan,
              nan],
        [     nan,      nan,      nan, ...,      nan,      nan,
              nan]],

       [[0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ],
        [0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ],
        [0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ],
...
        [0.      , 0.      , 0.      , ..., 0.015625, 0.015625,
         0.015625],
        [0.      , 0.      , 0.      , ..., 0.015625, 0.015625,
         0.015625],
        [0.      , 0.      , 0.      , ..., 0.015625, 0.015625,
         0.015625]],

       [[0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ],
        [0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ],
        [0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ],
        ...,
        [0.      , 0.      , 0.      , ..., 0.015625, 0.      ,
         0.015625],
        [0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ],
        [0.      , 0.      , 0.      , ..., 0.      , 0.      ,
         0.      ]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([              22.4,              22.45,               22.5,
                           22.55,               22.6,              22.65,
                            22.7,              22.75,               22.8,
                           22.85,
              ...
               47.14999999999999, 47.199999999999996,  47.24999999999999,
                            47.3, 47.349999999999994,  47.39999999999999,
              47.449999999999996,  47.49999999999999,              47.55,
              47.599999999999994],
             dtype='float64', name='latitude', length=505))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([   120.0, 120.0625,  120.125, 120.1875,   120.25, 120.3125,
               120.375, 120.4375,    120.5, 120.5625,
              ...
              149.4375,    149.5, 149.5625,  149.625, 149.6875,   149.75,
              149.8125,  149.875, 149.9375,    150.0],
             dtype='float64', name='longitude', length=481))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00', '2018-07-01 07:00:00',
               '2018-07-01 08:00:00', '2018-07-01 09:00:00',
               '2018-07-01 10:00:00', '2018-07-01 11:00:00',
               '2018-07-01 12:00:00', '2018-07-01 13:00:00',
               '2018-07-01 14:00:00', '2018-07-01 15:00:00',
               '2018-07-01 16:00:00', '2018-07-01 17:00:00',
               '2018-07-01 18:00:00', '2018-07-01 19:00:00',
               '2018-07-01 20:00:00', '2018-07-01 21:00:00',
               '2018-07-01 22:00:00', '2018-07-01 23:00:00',
               '2018-07-02 00:00:00', '2018-07-02 01:00:00',
               '2018-07-02 02:00:00', '2018-07-02 03:00:00',
               '2018-07-02 04:00:00', '2018-07-02 05:00:00',
               '2018-07-02 06:00:00', '2018-07-02 07:00:00',
               '2018-07-02 08:00:00', '2018-07-02 09:00:00',
               '2018-07-02 10:00:00', '2018-07-02 11:00:00',
               '2018-07-02 12:00:00', '2018-07-02 13:00:00',
               '2018-07-02 14:00:00', '2018-07-02 15:00:00',
               '2018-07-02 16:00:00', '2018-07-02 17:00:00',
               '2018-07-02 18:00:00', '2018-07-02 19:00:00',
               '2018-07-02 20:00:00', '2018-07-02 21:00:00'],
              dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0
In [31]:
# Script A12 #

da = ds["RH_1D5maboveground"]
da[-1,:,:].plot(cmap='RdYlGn_r',figsize=(5,5))
Out[31]:
<matplotlib.collections.QuadMesh at 0x1f92aa47490>
No description has been provided for this image
In [32]:
# Script A13 #

da.sel(latitude=35.6916,longitude=139.7516, method='nearest').plot.line(x="time",figsize=(12,3))
Out[32]:
[<matplotlib.lines.Line2D at 0x1f9295b97b0>]
No description has been provided for this image

D 推計気象分布¶

In [33]:
# Script A14 #

product = "obs_gpv"
grbdir = data_root/product/yyyymmdd[0:4]/yyyymmdd[0:6]
dtr = wx.dt_range("201807010600", 40, hours=1)
yyyymmddhhmm = [oo.strftime('%Y%m%d%H%M') for oo in dtr]
grblist= [Path(f"{grbdir}/Z__C_RJTD_{oo}00_OBS_GPV_Rjp_Ggis1km_Ptt_A{oo}_grib2.bin")
         for oo in yyyymmddhhmm]
print(wx.getgrvarname(grblist[0]))
grblist

['TMP']
Out[33]:
[WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701060000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807010600_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701070000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807010700_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701080000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807010800_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701090000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807010900_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701100000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011000_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701110000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011100_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701120000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011200_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701130000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011300_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701140000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011400_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701150000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011500_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701160000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011600_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701170000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011700_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701180000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011800_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701190000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807011900_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701200000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807012000_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701210000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807012100_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701220000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807012200_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180701230000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807012300_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702000000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020000_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702010000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020100_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702020000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020200_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702030000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020300_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702040000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020400_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702050000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020500_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702060000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020600_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702070000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020700_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702080000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020800_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702090000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807020900_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702100000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021000_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702110000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021100_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702120000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021200_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702130000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021300_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702140000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021400_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702150000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021500_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702160000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021600_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702170000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021700_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702180000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021800_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702190000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807021900_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702200000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807022000_grib2.bin'),
 WindowsPath('jmadata/obs_gpv/2018/201807/Z__C_RJTD_20180702210000_OBS_GPV_Rjp_Ggis1km_Ptt_A201807022100_grib2.bin')]
In [34]:
# Script A15 #

elements = ["TMP"]
ds = wx.getgpv(grblist[0],elements,to_netcdf=False) #最初の１つしか読んでいません
ds
Out[34]:
<xarray.Dataset>
Dimensions:      (latitude: 3360, longitude: 2560, time: 1)
Coordinates:
  * latitude     (latitude) float64 20.0 20.01 20.02 20.03 ... 47.98 47.99 48.0
  * longitude    (longitude) float64 118.0 118.0 118.0 ... 150.0 150.0 150.0
  * time         (time) datetime64[ns] 2018-07-01T06:00:00
Data variables:
    TMP_surface  (time, latitude, longitude) float32 nan nan nan ... nan nan nan
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 3360
    • longitude: 2560
    • time: 1
    • latitude
      (latitude)
      float64
      20.0 20.01 20.02 ... 47.99 48.0
      units :
      degrees_north
      long_name :
      latitude
      array([20.004176, 20.012509, 20.020843, ..., 47.979166, 47.9875  , 47.995833])
    • longitude
      (longitude)
      float64
      118.0 118.0 118.0 ... 150.0 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([118.00625, 118.01875, 118.03125, ..., 149.96875, 149.98125, 149.99375])
    • time
      (time)
      datetime64[ns]
      2018-07-01T06:00:00
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530424800.0
      reference_time_type :
      1
      reference_date :
      2018.07.01 06:00:00 UTC
      reference_time_description :
      analyses, reference date is fixed
      time_step_setting :
      auto
      time_step :
      0.0
      array(['2018-07-01T06:00:00.000000000'], dtype='datetime64[ns]')
    • TMP_surface
      (time, latitude, longitude)
      float32
      nan nan nan nan ... nan nan nan nan
      short_name :
      TMP_surface
      long_name :
      Temperature
      level :
      surface
      units :
      K
      array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([20.004175999999998,  20.01250933045549, 20.020842660910983,
              20.029175991366476, 20.037509321821968,  20.04584265227746,
              20.054175982732954, 20.062509313188446,  20.07084264364394,
               20.07917597409943,
              ...
               47.92083302590056, 47.929166356356056,  47.93749968681155,
               47.94583301726704, 47.954166347722534,  47.96249967817803,
               47.97083300863352,  47.97916633908901, 47.987499669544505,
                       47.995833],
             dtype='float64', name='latitude', length=3360))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([         118.00625,          118.01875,          118.03125,
                       118.04375,          118.05625, 118.06875000000001,
              118.08125000000001, 118.09375000000001, 118.10625000000002,
              118.11875000000002,
              ...
              149.88125000000724, 149.89375000000723, 149.90625000000725,
              149.91875000000726, 149.93125000000725, 149.94375000000724,
              149.95625000000726, 149.96875000000728, 149.98125000000726,
              149.99375000000725],
             dtype='float64', name='longitude', length=2560))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-01 06:00:00'], dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0
In [35]:
# Script A16 #

da = ds["TMP_surface"]
da[0,:,:].plot(cmap='RdYlGn_r',figsize=(5,5))
Out[35]:
<matplotlib.collections.QuadMesh at 0x1f92895a7d0>
No description has been provided for this image
In [36]:
# Script A17 #

# 全域をオブジェクト化するとサイズが巨大になるので領域を限定して全ファイルを読みます
area = (35.6916,35.6916,139.7516,139.7516)  #アメダス「東京」の緯度経度
ds = wx.getgpv(grblist,elements,to_netcdf=False,lalomima=area)

da = ds["TMP_surface"]
da.sel(latitude=35.6916,longitude=139.7516, method='nearest').plot.line(x="time",figsize=(12,3))
Out[36]:
[<matplotlib.lines.Line2D at 0x1f92ad6f190>]
No description has been provided for this image

D 解析雨量¶

In [37]:
# Script A18 #

product = "ra"
grbdir = data_root/product/yyyymmdd[0:4]/yyyymmdd[0:6]
dtr = wx.dt_range("201807060000", 10, minutes=30)
yyyymmddhhmm = [oo.strftime('%Y%m%d%H%M') for oo in dtr]
grblist= [Path(f"{grbdir}/Z__C_RJTD_{oo}00_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin")
         for oo in yyyymmddhhmm]
print(wx.getgrvarname(grblist[0]))
grblist

['var discipline=0 center=34 local_table=1 parmcat=1 parm=200']
Out[37]:
[WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706000000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706003000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706010000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706013000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706020000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706023000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706030000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706033000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706040000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin'),
 WindowsPath('jmadata/ra/2018/201807/Z__C_RJTD_20180706043000_SRF_GPV_Ggis1km_Prr60lv_ANAL_grib2.bin')]
In [38]:
# Script A19 #

elements = ["var discipline=0 center=34 local_table=1 parmcat=1 parm=200"]
ds = wx.getgpv(grblist[0],elements,to_netcdf=False) #最初の１つしか読んでいません
ds
Out[38]:
<xarray.Dataset>
Dimensions:             (latitude: 3360, longitude: 2560, time: 1)
Coordinates:
  * latitude            (latitude) float64 20.0 20.01 20.02 ... 47.98 47.99 48.0
  * longitude           (longitude) float64 118.0 118.0 118.0 ... 150.0 150.0
  * time                (time) datetime64[ns] 2018-07-06
Data variables:
    var0_1_200_surface  (time, latitude, longitude) float32 nan nan ... nan nan
Attributes:
    Conventions:          COARDS
    History:              created by wgrib2
    GRIB2_grid_template:  0
xarray.Dataset
    • latitude: 3360
    • longitude: 2560
    • time: 1
    • latitude
      (latitude)
      float64
      20.0 20.01 20.02 ... 47.99 48.0
      units :
      degrees_north
      long_name :
      latitude
      array([20.004167, 20.0125  , 20.020834, ..., 47.979166, 47.9875  , 47.995833])
    • longitude
      (longitude)
      float64
      118.0 118.0 118.0 ... 150.0 150.0
      units :
      degrees_east
      long_name :
      longitude
      array([118.00625, 118.01875, 118.03125, ..., 149.96875, 149.98125, 149.99375])
    • time
      (time)
      datetime64[ns]
      2018-07-06
      long_name :
      verification time generated by wgrib2 function verftime()
      reference_time :
      1530835200.0
      reference_time_type :
      1
      reference_date :
      2018.07.06 00:00:00 UTC
      reference_time_description :
      analyses, reference date is fixed
      time_step_setting :
      auto
      time_step :
      0.0
      array(['2018-07-06T00:00:00.000000000'], dtype='datetime64[ns]')
    • var0_1_200_surface
      (time, latitude, longitude)
      float32
      nan nan nan nan ... nan nan nan nan
      short_name :
      var0_1_200_surface
      long_name :
      desc
      level :
      surface
      units :
      unit
      array([[[nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        ...,
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan],
        [nan, nan, nan, ..., nan, nan, nan]]], dtype=float32)
    • latitude
      PandasIndex
      PandasIndex(Float64Index([20.004167000000002, 20.012500333134863, 20.020833666269723,
              20.029166999404588, 20.037500332539448,  20.04583366567431,
              20.054166998809173, 20.062500331944033, 20.070833665078894,
              20.079166998213754,
              ...
               47.92083300178624, 47.929166334921106,  47.93749966805596,
               47.94583300119083,  47.95416633432569,  47.96249966746055,
               47.97083300059541,  47.97916633373028,  47.98749966686513,
                       47.995833],
             dtype='float64', name='latitude', length=3360))
    • longitude
      PandasIndex
      PandasIndex(Float64Index([         118.00625,          118.01875,          118.03125,
                       118.04375,          118.05625, 118.06875000000001,
              118.08125000000001, 118.09375000000001, 118.10625000000002,
              118.11875000000002,
              ...
              149.88125000000724, 149.89375000000723, 149.90625000000725,
              149.91875000000726, 149.93125000000725, 149.94375000000724,
              149.95625000000726, 149.96875000000728, 149.98125000000726,
              149.99375000000725],
             dtype='float64', name='longitude', length=2560))
    • time
      PandasIndex
      PandasIndex(DatetimeIndex(['2018-07-06'], dtype='datetime64[ns]', name='time', freq=None))
  • Conventions :
    COARDS
    History :
    created by wgrib2
    GRIB2_grid_template :
    0
In [39]:
# Script A20 #

da = ds["var0_1_200_surface"]
da[0,:,:].plot(cmap='RdYlGn_r',figsize=(5,5))
Out[39]:
<matplotlib.collections.QuadMesh at 0x1f931947310>
No description has been provided for this image
In [40]:
# Script A21 #

# 全域をオブジェクト化するとサイズが巨大になるので領域を限定して全ファイルを読みます
area = (35.6916,35.6916,139.7516,139.7516)  #アメダス「東京」の緯度経度
ds = wx.getgpv(grblist,elements,to_netcdf=False,lalomima=area)

da = ds["var0_1_200_surface"]
da.sel(latitude=35.6916,longitude=139.7516, method='nearest').plot.line(x="time",figsize=(12,3))
Out[40]:
[<matplotlib.lines.Line2D at 0x1f931a2c700>]
No description has been provided for this image

著作権について¶

Copyright 2023 気象ビジネス推進コンソーシアム
(C) 2023 WXBC.

＜利用条件＞
本書は、本書に記載した要件・技術・方式に関する内容が変更されないこと、および出典を明示いただくことを前提に、無償でその全部または一部を複製、翻案、翻訳、転記、引用、公衆送信等して利用できます。なお、全体または一部を複製、翻案、翻訳された場合は、本書にある著作権表示および利用条件を明示してください。

＜免責事項＞
本書の著作権者は、本書の記載内容に関して、その正確性、商品性、利用目的への適合性等に関して保証するものではなく、特許権、著作権、その他の権利を侵害していないことを保証するものでもありません。本書の利用により生じた損害について、本書の著作権者は、法律上のいかなる責任も負いません。