2023年度 WXBC 人材育成WG テクノロジー研修
気象庁GPVデータ分析チャレンジ!基礎編¶
はじめに¶
WXBCテクノロジー研修「気象庁GPVデータ分析チャレンジ!入門」では、気象庁が対外的に作成するGPVプロダクトを紹介し、そのいくつかについて、配信媒体であるGRIBファイルからデータを取り出し、可視化する方法を学びました。気象庁のプロダクトへの理解を深めることに重点を置いたので、GPVを取り扱う上で必要な知識や処理テクニックに関して省いた事柄は少なくありません。この研修では、これらのうち、基礎的なものについて掘り下げて説明します。
1 GRIBファイル処理ツール wgrib2¶
気象庁のGPVプロダクトは、GRIB(GRIB/GRIB2)と呼ばれる特殊なフォーマットでファイルに収められています。気象庁は、プロダクトそれぞれについて、フォーマットの詳細を公開しており、また、これをデコードするためのC言語のサンプルコードを求めに応じ無保証で提供してくれる場合がありますが、基本的には、利用者が資料に基づいてプログラムを開発しデコードするものという立場をとっています。しかし、個々の利用者が同じものに個別にコストを掛けるのは全く非効率的なので、この研修では、アメリカ大気海洋局(NOAA)の気候予測センター(CPC)が提供する、wgrib2 というツール( https://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/index.html )でデコードをする方法を説明します。
wgrib2 は、以下の機能を持ちます。
- GRIBファイルの作成と読み出し
- データの一部取り出し
- 特定領域の取り出し
- 各種ファイル形式への変換(ieee, text, binary, CSV, netcdf, mysql)
- 新規データの追記
2 章で説明するライブラリwxbcbribXを使用すれば、私たちは wgrib2 を目にすることなく GRIB ファイルを取り扱えますが、その裏で働く大変重要なツールなので、この研修においては、最低限の使用法について学習しておきます。
GRIBファイルを操作するソフトウェアについては、このほかに、Pythonライブラリ pygrib が存在しますが、残念ながらWindowsでは利用できません(Windows用Fortranコンパイラを用意しソースコードを自らコンパイルしない限り)。そして、利用可能なOS上であっても、気象庁のプロダクトの一部(解析雨量、推計気象分布等)が取り扱えない(異常終了します)ため、選択肢には入りません。
ほかに、上述のNOAA/CPCでは、wgrib2 のPythonインターフェスの開発が進められているようです( https://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/pywgrib2.html )。まだ安定性が十分でないようですが、大いに期待されるところです。
1.1 コマンドラインでの実行¶
wgrib2 は、コマンドラインプログラムという種類に分類されるプログラムです。そこで、まずは、コマンドラインプログラム本来の使い方をざっと学習します。
コマンドラインプログラムは、プログラムのファイル名で始まる文字列を、コマンドプロンプト (Mac では ターミナル )と呼ばれる黒いウィンドウに書き込んで、最後に [Enter] キーを押すことで実行します。 プログラムのファイル名で始まる文字列は、コマンドライン と呼ばれています。
wgrib2で処理を行う場合、コマンドラインは概ね以下のような構造になります。なお、{操作の指示}の部分の文字列は、コマンドラインオプション と呼ばれることもあります。
"{wgrib2のファイル名}半角空白{操作の指示}半角空白{対象のGRIB2ファイル名}"
一番簡単な例として、wgrib2に、バージョンを表示させる場合を考えましょう。バージョンを表示させるコマンドラインオプションは、「-version」です。また、この処理にGRIB2ファイルは必要ないのでその部分はありません。よって、コマンドラインは以下となります。
c:\wgrib2\wgrib2.exe -version
これを実行すると、結果は下図のようになります。バージョンと製作者名が表示されて終了します。可能な人は、コマンドプロンプトを起動して上のテキストをコピペし、実行してみてください。 ここで、「wgrib2のファイル」として、「wgrib2.exe」ではなく「c:\wgrib2\wgrib2.exe」としたのは、このファイルがコンピューターのどの場所に置かれているかまでを間違いなく示すためです。このような「 置き場所+ファイル名 」のことを、パス(path) と言います。
wgrib2は、コマンドラインオプションを付けずに、対象のGRIBファイル名 だけを置くと、そのファイルに格納されているグリッドデータの一覧(インベントリー)を表示します。これを使って、サンプルとして用意したメソ数値予報モデルGPV(MSM-GPV)の GRIB ファイルを一つ開き、どのような気象データが中に記録されているを見ることにしましょう。
2023年6月1日0時UTCを初期値とする地上気象のGRIBファイルの中身を見ることにします。
{wgrib2のファイル名}半角空白{操作の指示}半角空白{対象のGRIB2ファイル名}
{wgrib2のファイル名} の部分:
この部分の文字列は、置き場所も含めて「c:\wgrib2\wgrib2.exe」です。{操作の指示} の部分:
この部分はありません。{対象のGRIBファイル名} の部分:
サンプルファイルの名は、 Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin で、
置かれている場所は、 C:\Users\user_name\challenge4\jmadata\msm\2023\202306 です。よって、全体としては以下となります。
C:\Users\user_name\challenge4\jmadata\msm\2023\202306\Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin
従って、コマンドラインは、これらを全部繋げた以下となります。これをコマンドラインに打ち込みます(長いですが続けて打ち込みます)。
c:\wgrib2\wgrib2.exe C:\Users\user_name\challenge4\jmadata\msm\2023\202306\Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin
その結果は下図の通りです。190ものインベントリーが表示されます(はじめの方は流れ去ってしまいます)。
可能な方は、これについても試してみてください。試すときは、コマンドライン中にある「user_name」の部分をご自身のアカウント名に変更して実行してください。
1.2 Python からの実行¶
ご覧(体験)いただいたとおり、コマンドラインプログラムを素で使うのはなかなかシンドイ作業です。そこで、Python の力を借りて、もっと楽に wgrib2 を制御しましょう。
それには、ライブラリ subprocess が提供する関数 run を利用します。
まず、ライブラリ subprocess をインポートします。以下を実行してください。
import subprocess
ライブラリ subprocess をインポートすると、関数 run が利用可能となるので、これを用いて以下のようにスクリプトを書きます。
rc = subprocess.run("command line 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 line string" は、コマンドラインに打ち込む文字列を表します。基本的には、コマンドラインをそのままクォーテーションマークで囲めばOKですが、Windows で使われている円マーク(バックスラッシュ)「¥」は半角スラッシュ[ / ]に置き換えてください。
スクリプトにはfor文が2つありますが、これらは、実行の過程でプログラムが発生させたメッセージやエラーをpython上で表示させるためのものです。
先に、コマンドプロンプトから実行したのと同じ処理をさせる python スクリプトは以下のようになります。
実行して結果を確認してください。
rc = subprocess.run("c:/wgrib2/wgrib2 jmadata/msm/2023/202306/Z__C_RJTD_20230601000000_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=2023060100:PRMSL:mean sea level:anl: 1.2:0:d=2023060100:PRES:surface:anl: 1.3:0:d=2023060100:UGRD:10 m above ground:anl: 1.4:0:d=2023060100:VGRD:10 m above ground:anl: 1.5:0:d=2023060100:TMP:1.5 m above ground:anl: 1.6:0:d=2023060100:RH:1.5 m above ground:anl: 1.7:0:d=2023060100:LCDC:surface:anl: 1.8:0:d=2023060100:MCDC:surface:anl: 1.9:0:d=2023060100:HCDC:surface:anl: 1.10:0:d=2023060100:TCDC:surface:anl: 1.11:0:d=2023060100:PRMSL:mean sea level:1 hour fcst: 1.12:0:d=2023060100:PRES:surface:1 hour fcst: 1.13:0:d=2023060100:UGRD:10 m above ground:1 hour fcst: 1.14:0:d=2023060100:VGRD:10 m above ground:1 hour fcst: 1.15:0:d=2023060100:TMP:1.5 m above ground:1 hour fcst: 1.16:0:d=2023060100:RH:1.5 m above ground:1 hour fcst: 1.17:0:d=2023060100:LCDC:surface:1 hour fcst: 1.18:0:d=2023060100:MCDC:surface:1 hour fcst: 1.19:0:d=2023060100:HCDC:surface:1 hour fcst: 1.20:0:d=2023060100:TCDC:surface:1 hour fcst: 1.21:0:d=2023060100:APCP:surface:0-1 hour acc fcst: 1.22:0:d=2023060100:DSWRF:surface:0-1 hour ave fcst: 1.23:0:d=2023060100:PRMSL:mean sea level:2 hour fcst: 1.24:0:d=2023060100:PRES:surface:2 hour fcst: 1.25:0:d=2023060100:UGRD:10 m above ground:2 hour fcst: 1.26:0:d=2023060100:VGRD:10 m above ground:2 hour fcst: 1.27:0:d=2023060100:TMP:1.5 m above ground:2 hour fcst: 1.28:0:d=2023060100:RH:1.5 m above ground:2 hour fcst: 1.29:0:d=2023060100:LCDC:surface:2 hour fcst: 1.30:0:d=2023060100:MCDC:surface:2 hour fcst: 1.31:0:d=2023060100:HCDC:surface:2 hour fcst: 1.32:0:d=2023060100:TCDC:surface:2 hour fcst: 1.33:0:d=2023060100:APCP:surface:1-2 hour acc fcst: 1.34:0:d=2023060100:DSWRF:surface:1-2 hour ave fcst: 1.35:0:d=2023060100:PRMSL:mean sea level:3 hour fcst: 1.36:0:d=2023060100:PRES:surface:3 hour fcst: 1.37:0:d=2023060100:UGRD:10 m above ground:3 hour fcst: 1.38:0:d=2023060100:VGRD:10 m above ground:3 hour fcst: 1.39:0:d=2023060100:TMP:1.5 m above ground:3 hour fcst: 1.40:0:d=2023060100:RH:1.5 m above ground:3 hour fcst: 1.41:0:d=2023060100:LCDC:surface:3 hour fcst: 1.42:0:d=2023060100:MCDC:surface:3 hour fcst: 1.43:0:d=2023060100:HCDC:surface:3 hour fcst: 1.44:0:d=2023060100:TCDC:surface:3 hour fcst: 1.45:0:d=2023060100:APCP:surface:2-3 hour acc fcst: 1.46:0:d=2023060100:DSWRF:surface:2-3 hour ave fcst: 1.47:0:d=2023060100:PRMSL:mean sea level:4 hour fcst: 1.48:0:d=2023060100:PRES:surface:4 hour fcst: 1.49:0:d=2023060100:UGRD:10 m above ground:4 hour fcst: 1.50:0:d=2023060100:VGRD:10 m above ground:4 hour fcst: 1.51:0:d=2023060100:TMP:1.5 m above ground:4 hour fcst: 1.52:0:d=2023060100:RH:1.5 m above ground:4 hour fcst: 1.53:0:d=2023060100:LCDC:surface:4 hour fcst: 1.54:0:d=2023060100:MCDC:surface:4 hour fcst: 1.55:0:d=2023060100:HCDC:surface:4 hour fcst: 1.56:0:d=2023060100:TCDC:surface:4 hour fcst: 1.57:0:d=2023060100:APCP:surface:3-4 hour acc fcst: 1.58:0:d=2023060100:DSWRF:surface:3-4 hour ave fcst: 1.59:0:d=2023060100:PRMSL:mean sea level:5 hour fcst: 1.60:0:d=2023060100:PRES:surface:5 hour fcst: 1.61:0:d=2023060100:UGRD:10 m above ground:5 hour fcst: 1.62:0:d=2023060100:VGRD:10 m above ground:5 hour fcst: 1.63:0:d=2023060100:TMP:1.5 m above ground:5 hour fcst: 1.64:0:d=2023060100:RH:1.5 m above ground:5 hour fcst: 1.65:0:d=2023060100:LCDC:surface:5 hour fcst: 1.66:0:d=2023060100:MCDC:surface:5 hour fcst: 1.67:0:d=2023060100:HCDC:surface:5 hour fcst: 1.68:0:d=2023060100:TCDC:surface:5 hour fcst: 1.69:0:d=2023060100:APCP:surface:4-5 hour acc fcst: 1.70:0:d=2023060100:DSWRF:surface:4-5 hour ave fcst: 1.71:0:d=2023060100:PRMSL:mean sea level:6 hour fcst: 1.72:0:d=2023060100:PRES:surface:6 hour fcst: 1.73:0:d=2023060100:UGRD:10 m above ground:6 hour fcst: 1.74:0:d=2023060100:VGRD:10 m above ground:6 hour fcst: 1.75:0:d=2023060100:TMP:1.5 m above ground:6 hour fcst: 1.76:0:d=2023060100:RH:1.5 m above ground:6 hour fcst: 1.77:0:d=2023060100:LCDC:surface:6 hour fcst: 1.78:0:d=2023060100:MCDC:surface:6 hour fcst: 1.79:0:d=2023060100:HCDC:surface:6 hour fcst: 1.80:0:d=2023060100:TCDC:surface:6 hour fcst: 1.81:0:d=2023060100:APCP:surface:5-6 hour acc fcst: 1.82:0:d=2023060100:DSWRF:surface:5-6 hour ave fcst: 1.83:0:d=2023060100:PRMSL:mean sea level:7 hour fcst: 1.84:0:d=2023060100:PRES:surface:7 hour fcst: 1.85:0:d=2023060100:UGRD:10 m above ground:7 hour fcst: 1.86:0:d=2023060100:VGRD:10 m above ground:7 hour fcst: 1.87:0:d=2023060100:TMP:1.5 m above ground:7 hour fcst: 1.88:0:d=2023060100:RH:1.5 m above ground:7 hour fcst: 1.89:0:d=2023060100:LCDC:surface:7 hour fcst: 1.90:0:d=2023060100:MCDC:surface:7 hour fcst: 1.91:0:d=2023060100:HCDC:surface:7 hour fcst: 1.92:0:d=2023060100:TCDC:surface:7 hour fcst: 1.93:0:d=2023060100:APCP:surface:6-7 hour acc fcst: 1.94:0:d=2023060100:DSWRF:surface:6-7 hour ave fcst: 1.95:0:d=2023060100:PRMSL:mean sea level:8 hour fcst: 1.96:0:d=2023060100:PRES:surface:8 hour fcst: 1.97:0:d=2023060100:UGRD:10 m above ground:8 hour fcst: 1.98:0:d=2023060100:VGRD:10 m above ground:8 hour fcst: 1.99:0:d=2023060100:TMP:1.5 m above ground:8 hour fcst: 1.100:0:d=2023060100:RH:1.5 m above ground:8 hour fcst: 1.101:0:d=2023060100:LCDC:surface:8 hour fcst: 1.102:0:d=2023060100:MCDC:surface:8 hour fcst: 1.103:0:d=2023060100:HCDC:surface:8 hour fcst: 1.104:0:d=2023060100:TCDC:surface:8 hour fcst: 1.105:0:d=2023060100:APCP:surface:7-8 hour acc fcst: 1.106:0:d=2023060100:DSWRF:surface:7-8 hour ave fcst: 1.107:0:d=2023060100:PRMSL:mean sea level:9 hour fcst: 1.108:0:d=2023060100:PRES:surface:9 hour fcst: 1.109:0:d=2023060100:UGRD:10 m above ground:9 hour fcst: 1.110:0:d=2023060100:VGRD:10 m above ground:9 hour fcst: 1.111:0:d=2023060100:TMP:1.5 m above ground:9 hour fcst: 1.112:0:d=2023060100:RH:1.5 m above ground:9 hour fcst: 1.113:0:d=2023060100:LCDC:surface:9 hour fcst: 1.114:0:d=2023060100:MCDC:surface:9 hour fcst: 1.115:0:d=2023060100:HCDC:surface:9 hour fcst: 1.116:0:d=2023060100:TCDC:surface:9 hour fcst: 1.117:0:d=2023060100:APCP:surface:8-9 hour acc fcst: 1.118:0:d=2023060100:DSWRF:surface:8-9 hour ave fcst: 1.119:0:d=2023060100:PRMSL:mean sea level:10 hour fcst: 1.120:0:d=2023060100:PRES:surface:10 hour fcst: 1.121:0:d=2023060100:UGRD:10 m above ground:10 hour fcst: 1.122:0:d=2023060100:VGRD:10 m above ground:10 hour fcst: 1.123:0:d=2023060100:TMP:1.5 m above ground:10 hour fcst: 1.124:0:d=2023060100:RH:1.5 m above ground:10 hour fcst: 1.125:0:d=2023060100:LCDC:surface:10 hour fcst: 1.126:0:d=2023060100:MCDC:surface:10 hour fcst: 1.127:0:d=2023060100:HCDC:surface:10 hour fcst: 1.128:0:d=2023060100:TCDC:surface:10 hour fcst: 1.129:0:d=2023060100:APCP:surface:9-10 hour acc fcst: 1.130:0:d=2023060100:DSWRF:surface:9-10 hour ave fcst: 1.131:0:d=2023060100:PRMSL:mean sea level:11 hour fcst: 1.132:0:d=2023060100:PRES:surface:11 hour fcst: 1.133:0:d=2023060100:UGRD:10 m above ground:11 hour fcst: 1.134:0:d=2023060100:VGRD:10 m above ground:11 hour fcst: 1.135:0:d=2023060100:TMP:1.5 m above ground:11 hour fcst: 1.136:0:d=2023060100:RH:1.5 m above ground:11 hour fcst: 1.137:0:d=2023060100:LCDC:surface:11 hour fcst: 1.138:0:d=2023060100:MCDC:surface:11 hour fcst: 1.139:0:d=2023060100:HCDC:surface:11 hour fcst: 1.140:0:d=2023060100:TCDC:surface:11 hour fcst: 1.141:0:d=2023060100:APCP:surface:10-11 hour acc fcst: 1.142:0:d=2023060100:DSWRF:surface:10-11 hour ave fcst: 1.143:0:d=2023060100:PRMSL:mean sea level:12 hour fcst: 1.144:0:d=2023060100:PRES:surface:12 hour fcst: 1.145:0:d=2023060100:UGRD:10 m above ground:12 hour fcst: 1.146:0:d=2023060100:VGRD:10 m above ground:12 hour fcst: 1.147:0:d=2023060100:TMP:1.5 m above ground:12 hour fcst: 1.148:0:d=2023060100:RH:1.5 m above ground:12 hour fcst: 1.149:0:d=2023060100:LCDC:surface:12 hour fcst: 1.150:0:d=2023060100:MCDC:surface:12 hour fcst: 1.151:0:d=2023060100:HCDC:surface:12 hour fcst: 1.152:0:d=2023060100:TCDC:surface:12 hour fcst: 1.153:0:d=2023060100:APCP:surface:11-12 hour acc fcst: 1.154:0:d=2023060100:DSWRF:surface:11-12 hour ave fcst: 1.155:0:d=2023060100:PRMSL:mean sea level:13 hour fcst: 1.156:0:d=2023060100:PRES:surface:13 hour fcst: 1.157:0:d=2023060100:UGRD:10 m above ground:13 hour fcst: 1.158:0:d=2023060100:VGRD:10 m above ground:13 hour fcst: 1.159:0:d=2023060100:TMP:1.5 m above ground:13 hour fcst: 1.160:0:d=2023060100:RH:1.5 m above ground:13 hour fcst: 1.161:0:d=2023060100:LCDC:surface:13 hour fcst: 1.162:0:d=2023060100:MCDC:surface:13 hour fcst: 1.163:0:d=2023060100:HCDC:surface:13 hour fcst: 1.164:0:d=2023060100:TCDC:surface:13 hour fcst: 1.165:0:d=2023060100:APCP:surface:12-13 hour acc fcst: 1.166:0:d=2023060100:DSWRF:surface:12-13 hour ave fcst: 1.167:0:d=2023060100:PRMSL:mean sea level:14 hour fcst: 1.168:0:d=2023060100:PRES:surface:14 hour fcst: 1.169:0:d=2023060100:UGRD:10 m above ground:14 hour fcst: 1.170:0:d=2023060100:VGRD:10 m above ground:14 hour fcst: 1.171:0:d=2023060100:TMP:1.5 m above ground:14 hour fcst: 1.172:0:d=2023060100:RH:1.5 m above ground:14 hour fcst: 1.173:0:d=2023060100:LCDC:surface:14 hour fcst: 1.174:0:d=2023060100:MCDC:surface:14 hour fcst: 1.175:0:d=2023060100:HCDC:surface:14 hour fcst: 1.176:0:d=2023060100:TCDC:surface:14 hour fcst: 1.177:0:d=2023060100:APCP:surface:13-14 hour acc fcst: 1.178:0:d=2023060100:DSWRF:surface:13-14 hour ave fcst: 1.179:0:d=2023060100:PRMSL:mean sea level:15 hour fcst: 1.180:0:d=2023060100:PRES:surface:15 hour fcst: 1.181:0:d=2023060100:UGRD:10 m above ground:15 hour fcst: 1.182:0:d=2023060100:VGRD:10 m above ground:15 hour fcst: 1.183:0:d=2023060100:TMP:1.5 m above ground:15 hour fcst: 1.184:0:d=2023060100:RH:1.5 m above ground:15 hour fcst: 1.185:0:d=2023060100:LCDC:surface:15 hour fcst: 1.186:0:d=2023060100:MCDC:surface:15 hour fcst: 1.187:0:d=2023060100:HCDC:surface:15 hour fcst: 1.188:0:d=2023060100:TCDC:surface:15 hour fcst: 1.189:0:d=2023060100:APCP:surface:14-15 hour acc fcst: 1.190:0:d=2023060100:DSWRF:surface:14-15 hour ave fcst:
コマンドプロンプトに打ち込んだ時と同じ結果が得られました。ライブラリ subprocess のおかげで wgrib2 の制御が楽になったので、これを利用して wgrib2 の機能についていろいろとみてゆきましょう。
1.3 GPVのインベントリー¶
まず、GRIBファイルから取り出したインベントリー情報について見直してみましょう。インベントリーとは、そのファイルに格納されているGPVデータの概要のことです。出力結果を見ると、GRIB ファイルの中には非常に多くのGPVデータが格納されていて、総数は190 に上ることがわかります。
それぞれのインベントリーは、文字列がコロン「:」で区切られた形をして、いくつかは繰り返し登場ようなので、Pythonに少し働いてもらい、それぞれのセクションに現れる文字列を集めて重複を除き、どのような文字列が使われているかを調べてみようと思います。またこの際、関数 run に渡すコマンドライン文字列がとても長くてみずらいので、変数等の力も借りてスクリプトを少し見通し良くします。
以下を実行してください。
wgrib2 = "c:/wgrib2/wgrib2"
# wgrib2 = "~/work/grib2/wgrib2/wgrib2" # Macの場合
grbdir = "jmadata/msm/2023/202306"
grbfile = "Z__C_RJTD_20230601000000_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=2023060100'] ---------- ['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'] ----------
このGRIBファイルのインベントリーは、「連番:ゼロ:d=初期値の日時:気象要素記号:高度:予報時間」であることがわかりました。
1.4 GPVの選択¶
wgrib2 の機能において、GRIBファイルの中の特定のデータに関する処理の場合は、どのデータに と 何をする の両方を {操作の指示} に込めます。
{wgrib2のファイル名}半角空白{操作の指示}半角空白{対象のGRIB2ファイル名}
wgrib2 では、コマンドラインオプション「 -match_fs "検索文字列" 」を使用すると どのデータに を指定することができます。より複雑な検索に対応できるよう、オプション -match "検索正規表現" も用意されています。これらのオプションは複数個使うことができ、「 且つ 」として扱われます。 これを使って、インベントリー183番目(下から6番目)のGPVデータを抜き出してみましょう。抜き出して「何をする」か、ですが、これについては、最も簡単なものとして、-V を使ってみましょう。GPVの概要を表示します。
この場合のコマンドラインオプション({操作の指示})は、「 -match_fs "1.183:" -V 」になります。
従って、MSM-GPVのGRIBファイルから、第183番目のGPVを抜き出し、その概要を表示するスクリプトは以下のようになります。実行してください。
# wgrib2のファイル
wgrib2 = 'c:/wgrib2/wgrib2.exe'
#wgrib2 = "~/work/grib2/wgrib2/wgrib2" # Macの場合
# オプション(操作の指示)
kwds = '-match_fs "1.183:" -V' #文字列中に「"」を入れたいときは、クォーテーションに「'」を使います
# 対象のファイル
grbdir = 'jmadata/msm/2023/202306'
grbfile = 'Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin'
grbpath = grbdir +'/'+ grbfile # ファイルへのパス
# wgrib2の実行
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=2023060115:1.5 m above ground:15 hour fcst:TMP Temperature [K]: ndata=242905:undef=0:mean=292.568:min=275.627:max=302.439 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
出力を読むと、このデータは、初期時刻(2023年6月1日00UTC時)の15時間後における地上1.5 mの気温を予測したもので、単位はK(ケルビン)、であり、北緯22.4~47.6度、東経120~150度の範囲を、緯度方向に0.05度、経度方向には0.0625度の間隔で張ったGPVであることがわかります。
1.5 GPVの取り出し¶
オプション -csv 出力ファイル名 は、データをCSV形式のファイルとして出力させる指示を wgrib2 に与えます。
そこで、選択の指示とCSV出力の指示を並べて、インベントリ番号1.183のデータをCSV ファイルに出力してみましょう。以下のようになります。実行すると、フォルダchallenge4 に、ファイル TMP.csv が新規作成されます。
# wgrib2のファイル
wgrib2 = 'c:/wgrib2/wgrib2.exe' # Windowsの場合
# wgrib2 = '~/work/grib2/wgrib2/wgrib2' # Macの場合
# オプション(操作の指示)
kwds = '-match_fs "1.183:" -csv TMP.csv'
# 対象のファイル
grbdir = 'jmadata/msm/2023/202306'
grbfile = 'Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin'
grbpath = grbdir +'/'+ grbfile
# wgrib2の実行
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=2023060100:TMP:1.5 m above ground:15 hour fcst:
新規作成された TMP.csv を表計算ソフトやテキストエディタで表示し、中身を確認してください。下図のような内容になっているはずです。
マイクロソフト・エクセルで表示させた方は、経度がE列、緯度がF列に示されます。シートの上の方にある低緯度の気温と、シート下の方のにある高緯度の気温を見比べて、低緯度の格子の気温の方が高緯度の格子の気温より高い傾向にあることを確認してください。
コマンドラインオプション「 -match_fs 」で、インベントリに「 "TMP" 」が含まれるものを選択すれば、全ての予報期間の気温のデータをCSVファイルとして出力することができます。また、このオプションを用いずに、オプション -csv 出力ファイル名 を用いれば、全データがファイル出力されます。すなわち、GRIB2ファイルのフォーマット変換が行われることになります。
wgrib2は、CSVのほかに、フラットバイナリやMySQL、NetCDFなどのフォーマットでデータを出力することができます。今度は、このファイルの全部のデータをNetCDFファイルで出力してみましょう。
以下を実行してください。
# wgrib2のファイル
wgrib2 = 'c:/wgrib2/wgrib2.exe' # Windowsの場合
# wgrib2 = "~/work/grib2/wgrib2/wgrib2" # Macの場合
# 対象のファイル
grbdir = 'jmadata/msm/2023/202306'
grbfile = 'Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin'
grbpath = grbdir +'/'+ grbfile
# オプション(操作の指示)
kwds = f'-netcdf nc/{grbfile}.nc'
# wgrib2の実行
rc = subprocess.run(f"{wgrib2} {kwds} {grbpath}",
shell=True, text=True, capture_output=True)
#for line in rc.stdout.splitlines():
# print(line)
ファイルエクスプローラーで、フォルダ challenge4/nc の中を確認してください。長い名前で最後が「 .nc 」で終わるファイルが新たにできているはずです。ファイルサイズを見ると、182Mバイトとあります。元のGRIBファイルも67Mバイトとなかなかの大きさですが、NetCDFファイルにするとさらに大きなサイズになります。
wgrib2には、今説明した -match_fs や -V、-csv 以外にもたくさんのコマンドラインオプションがあるので、最後に、それらを表示させてみましょう。以下を実行してください。
# wgrib2のファイル
wgrib2 = "c:/wgrib2/wgrib2" # Windowsの場合
# wgrib2 = "~/work/grib2/wgrib2/wgrib2" # Macの場合
# オプション(操作の指示)
kwds = '-h'
# wgrib2の実行
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
2 WXBCオリジナルライブラリ wxbcgribX¶
wxbcgribX は、GPVデータを取り扱う上で必要な定型的な処理をひとまとめにした関数を提供するライブラリで、WXBC人材育成WGの有志が開発したものです。wxbcgribX は小規模なライブラリなので、その実体である ファイル wxbcgribx.py をプログラムと同じフォルダ内に置いておけば、インポートして利用できます。この講習の事前準備において、いくつかのPythonライブラリをインストールしましたが、そのような作業は不要です。
ただし、wxbcgribX のいくつかの関数は、内部的にソフトウエア wgrib2 を使用するため、これをインストールしたうえで、ファイル wxbcgribx.py の35行目に wgrib2 がインストールされているパスを書き込む必要があります。
wxbcgribX のライセンスは下記のとおりです。
Copyright (c) 2022, 2023, 2024 気象ビジネス推進コンソーシアム
以下に定める条件に従い、本ソフトウェアおよび関連文書のファイル(以下「ソフトウェア」)の複製を取得するすべての人に対し、ソフトウェアを無制限に扱うことを無償で許可します。これには、ソフトウェアの複製を使用、複写、変更、結合、掲載、頒布、サブライセンス、および/または販売する権利、およびソフトウェアを提供する相手に同じことを許可する権利も無制限に含まれます。
上記の著作権表示および本許諾表示を、ソフトウェアのすべての複製または重要な部分に記載するものとします。
ソフトウェアは「現状のまま」で、明示であるか暗黙であるかを問わず、何らの保証もなく提供されます。ここでいう保証とは、商品性、特定の目的への適合性、および権利非侵害についての保証も含みますが、それに限定されるものではありません。 作者または著作権者は、契約行為、不法行為、またはそれ以外であろうと、ソフトウェアに起因または関連し、あるいはソフトウェアの使用またはその他の扱いによって生じる一切の請求、損害、その他の義務について何らの責任も負わないものとします。
以下に、ライブラリが提供する主な関数とその概要を示します。
- getvarname(src)
指定したGRIBファイルのインベントリで使用されている気象要素の記号を得る関数- 引数:
- src: 対象とするGRIBファイルへのパス
- 戻り値:気象要素の識別子のリスト
- 引数:
- getensname(src)
指定したGRIBファイルに格納されるアンサンブルメンバーの識別子を得る関数- 引数:
- src: 対象とするGRIBファイルへのパス
- 戻り値:メンバー識別子のリスト
- 引数:
- getgpv(grblist, elements,
to_netcdf=False, from_netcdf=False, ncdir="./nc",
lalomima=None, timezone=None, verbose=False)
指定したGRIBファイルから指定した気象要素のGPVを取り出す関数- 引数:
- grblist: 対象とするGRIBファイルへのパスのリスト
- elements: GRIBファイルから取り出す気象要素
インベントリで使用されている文字列をリストで与える(['TMP','RH']など)。 - ncdir: 作成したDataDetオブジェクトをNetCDFファイルとして保管する場所
- to_netcdf: 作成したDataDetオブジェクトをNetCDFファイルとして保管するかどうかの指定
引数を省略するかFalseを与えると保管しない。Trueを指定すると保管される。 - from_netcdf: 以前に保管したNetCDFファイルが存在するときそれを読み込むかどうかの指定
引数を省略するかFalseを与えると読み込まずGRIB2ファイルから取り出す。Trueを指定すると読み込む。 - lalomima:データの緯度経度範囲を限定するときに与える範囲
4要素リスト [緯度の最大値,緯度の最小値,経度の最大値,経度の最小値] で与える。
GRIBから読みだす段階では全域、Datasetオブジェクトに段階でトリミングされる。 - timezone: 戻り値に標準時を追加したいときにタイムゾーンID('Asia/Tokyo'など)を指定
指定すると座標 awaretime が追加される。 - vorbose: デコード中のメンバー名と枚数を逐次表示するかどうかを指定
Trueを与えると表示する。何も指定しないと表示しない。
- grblist: 対象とするGRIBファイルへのパスのリスト
- 戻り値:xarray.DataSetオブジェクト
- 引数:
- awtime(xarr,tzid='UTC')
DataArray または Datasetオブジェクトの座標 time をUTCとみなし、その時刻に相当する aware な pandas.timestamp オブジェクトのリストを作成する関数- 引数:
- xarr: 時刻座標値を読み出す DataArray/Dataset オブジェクト
- tzid: タイムゾーンID('Asia/Tokyo’など)省略した場合はUTCと見なす
- 戻り値:指定したタイムゾーン表現の aware な pandas.timestamp オブジェクトのリスト
- 引数:
- strft_range(start='202001010000', end=None, format='%Y%m%d%H%M', iter=1, **step)
連続した日時の文字列のリストを作成する関数- 文字列の書式は format で与える
- 開始日時は、start にdatetimeオブジェクトか文字列(書式はformat)で指定する。
- 終了日時は、end にdatetimeオブジェクトか文字列(書式はformat)で指定する。
- itera に個数(デフォルトは1)を指定することで、end を省略することができる。
- end と itera が共に与えられた場合は、itera が無視される。
- 刻む時間間隔は「hours=1」(デフォルト)、「minutes=30」などの形で指定する。
- 戻り値:日時を表現する文字列のリスト
3 時刻の取り扱い¶
気象 × ○○ の分析をする際、○○ のデータはほぼ間違いなく日本標準時(JST)で整理されています。一方、気象庁の全てのGPVプロダクトは、時刻に中央協定時(UTC)を使用しているので、単純に並べてしまうと下図のように9時間ずれてしまいます。従って、気象庁のGPVプロダクト分析に使用する際には、異なる二つの時刻表記に折り合いを付けなければなりません。
図 アメダス「東京」における気温と、そこに最も近い格子点のGSM-GPV予測値を単純に並べたグラフ
どのように折り合いをつけるかはケースバイケースですが、ここでは、タイムゾーンを導入して 処理する方法を学びます。
3.1 aware な時刻¶
当たり前のことですが、「午前6時」と聞いたとき、私たちは全員、日本時間の午前6時と受け止めます。日常生活はこれで全く問題ありません。しかし、その一方で、国際線旅客機の発着時刻を確認するときは「日本時間の午前6時」という表現にして、注意深くコミュニケーションを図ることもあります。このような状況はコンピューティングの世界にも存在し、両者を区別したいとき、前者を naiveな時刻、後者を awareな時刻 と呼んでいます。もちろん、Pythonでも言語仕様としては両者を区別でき、どちらを使うこともできます。詳細は( https://docs.python.org/3/library/datetime.html )を参照してください。
言語としては使い分けられるものの、ライブラリの一部には使い分けられないものも存在します。そして、残念なことに、私たちがここで使っている xarray のオブジェクトは使い分けができない部類です。この問題をカバーするため、wxbcgribXの関数 getgpv には、 awareな時刻 の座標を付加する機能を持たせてあります。
それでは、この機能の使用方法を学びましょう。素材として、2023年6月1日00時UTCを初期時刻とするGSM-GPV を使用します。準備として、ライブラリーをインポートし、GRIBファイルへのパスを作成します。
# ライブラリのインポート
import pandas as pd # 時刻処理に便利なライブラリ
import wxbcgribx as wx
# ファイルへのパスのリスト作成
yyyymmdd = '20230601'
prdct = 'gsm'
fd = ['0000-0100','0101-0200','0201-0300','0301-0400','0401-0500','0501-0512',
'0515-0700','0703-0900','0903-1100']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD{xx}_grib2.bin'
for xx in fd]
GPV に aware な時刻座標を追加するには、関数 getgpv にキーワード引数 timezone を追加し タイムゾーンID ('Asia/Tokyo' などの文字列) を与えます。
ds = wx.getgpv(grb_paths, 'TMP', timezone='Asia/Tokyo') # 引数 timezone を使用
da = ds['TMP_2maboveground']
以下を実行して、概要フォームを表示させてください(概要フォームが表示されたら、見やすくするためにデータ本体の "お皿" をクリックし、表示を畳んでください)
da
<xarray.DataArray 'TMP_2maboveground' (time: 177, latitude: 301, longitude: 241)> array([[[301.50333, 301.4799 , 301.44864, ..., 301.1205 , 301.12833, 301.13614], [301.3705 , 301.33145, 301.2924 , ..., 301.1205 , 301.12833, 301.12833], [301.21426, 301.1752 , 301.15176, ..., 301.1127 , 301.1127 , 301.1205 ], ..., [286.13614, 286.08145, 286.07364, ..., 275.94083, 276.00333, 276.07364], [285.9799 , 285.90958, 285.90176, ..., 275.90176, 275.95645, 276.02676], [286.14395, 286.0502 , 285.9955 , ..., 275.8549 , 275.9174 , 275.9877 ]], [[301.5029 , 301.47165, 301.42477, ..., 301.21384, 301.22165, 301.22946], [301.36227, 301.3076 , 301.23727, ..., 301.20602, 301.22165, 301.22946], [301.20602, 301.1279 , 301.08102, ..., 301.1982 , 301.21384, 301.21384], ... [280.60638, 280.37982, 280.03607, ..., 278.02826, 278.0595 , 278.07513], [280.71576, 280.45795, 280.06732, ..., 277.95795, 277.96576, 277.97357], [280.69232, 280.41107, 279.97357, ..., 277.91107, 277.91888, 277.91888]], [[301.95477, 301.94696, 301.93915, ..., 301.05634, 301.0407 , 301.01727], [301.94696, 301.93915, 301.92352, ..., 301.07977, 301.05634, 301.0251 ], [301.93134, 301.93915, 301.92352, ..., 301.0954 , 301.07196, 301.0329 ], ..., [287.67352, 287.61102, 287.54852, ..., 278.4626 , 278.48602, 278.49384], [287.68915, 287.62665, 287.54852, ..., 278.45477, 278.4704 , 278.48602], [287.76727, 287.69696, 287.5954 , ..., 278.43134, 278.45477, 278.4704 ]]], dtype=float32) Coordinates: * latitude (latitude) float64 20.0 20.1 20.2 20.3 ... 49.7 49.8 49.9 50.0 * longitude (longitude) float64 120.0 120.1 120.2 120.4 ... 149.8 149.9 150.0 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-12 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-12T09:00:0... Attributes: short_name: TMP_2maboveground long_name: Temperature level: 2 m above ground units: K
Coordinates (座標) のセクションを着目してください。4行目にawaretime が存在します。この行の右端の "お皿" アイコンをクリックしてみてください。2023年6月1日00時UTCより9時間進んだ時刻が刻まれた配列が表示されます。そして、タイムゾーンID 'Azia/Tokyo' が付加しており、時刻はawareであることも確認できます。
次に、3行目 time に対して同様にして座標値を見てみましょう。こちらは、UTCの時刻が表示されていますが、これには準拠するタイムゾーンの表示がなく、naive な時刻だということが確認できます。
3.1.1 日本時間の折れ線グラフ¶
awaretime に保持される座標値の配列をオブジェクトから取り出して、グラフの x 軸データとして用いれば、時系列グラフの横軸を日本標準時(JST)で表示することができます。グラフをメソッド plot で描く場合であれば、キーワード引数 x に、座標名「'awaretime'」を代入します。
以下に、アメダス「舘野(つくば)」に最も近い格子点における気温予測値の折れ線グラフを、メソッドplotを用いて JST を横軸にして描くスクリプトを示します。
lat, lon = 36.05666667, 140.125 #「舘野(つくば)」
da.sel(latitude=lat, longitude=lon, method='nearest').plot(x='awaretime', figsize=(10,3))
[<matplotlib.lines.Line2D at 0x29106f6eef0>]
3.1.2 日本時間でのGPVの選択¶
GPVをDataArrayオブジェクトとして取り扱うことの大きなメリットの一つに、上の例でも使用する、メソッド sel の存在がありますが、大変残念なことに、メソッド sel で使用できる座標は、次元の構成に使われている座標に限定 されています。このため、例えば、日本標準時の2023年06月02日00時 の2次元データをGPVから取り出すときに、以下のようにすることはできません。
da.sel(awaretime='2023-06-02T00:00:00+0900') # これは不可
UTC 表現の時刻を使用するしかないのですが、aware な時刻は異なるタイムゾーンの時刻表現を簡単に作り出せるので、以下の手順を踏むことで、日本標準時の文字列やオブジェクトでメソッド sel が使用できます。
日本標準時の文字列 → 日本標準時のawareなオブジェクト → タイムゾーンがUTCに変更されたオブジェクト → UTCのnaiveなオブジェクト
以下に、2023年06月02日00時JSTにおける気温予測値をこの方法で抽出し、分布図とするスクリプトを示します。
jststr = '2023-06-02T00:00+0900' # 日本標準時の文字列
print(jststr)
awjst = pd.to_datetime(jststr) # 文字列をawareなオブジェクトに
print(awjst)
awutc = awjst.tz_convert('utc') # タイムゾーンをUTCに変更
print(awutc)
nvutc = awutc.tz_localize(None) # タイムゾーン情報を除去
print(nvutc)
da.sel(time=nvutc).plot(cmap='RdYlGn_r')
2023-06-02T00:00+0900 2023-06-02 00:00:00+09:00 2023-06-01 15:00:00+00:00 2023-06-01 15:00:00
<matplotlib.collections.QuadMesh at 0x29107172680>
この変換は機械的なものなので、ライブラリ wxbcgribX に、関数 nvutc として用意してあります。これを使えば、上のスクリプトは、以下のように簡略化できます。
jststr = '2023-06-02T00:00+0900' # 日本標準時の文字列
da.sel(time=wx.nvutc(jststr)).plot(cmap='RdYlGn_r')
3.2 日本標準時のデータとの突き合わせ¶
日本標準時で整理されているデータと GPV データとを aware な時刻を利用して突き合わせる方法、つまり、
[時刻],[日本標準時で整理されているデータ],[GPVのデータ]を、以下のような3列の表にまとめる方法を考えてみましょう。
時刻 | 突合対象のデータ | GPVのデータ |
---|---|---|
2023-06-02 00:00 | : | : |
2023-06-02 00:00 | : | : |
2023-06-02 00:00 | : | : |
: | : | : |
日本標準時で整理されているデータ として、今回はアメダスデータを使うことにします。GPVのデータ としては、3.1 章で使用した GSM-GPV を使用します(これだと 気象×気象 ですが許してください)。
3.2.1 pandas.DataFrame オブジェクト¶
表形式に整理されたデータをPythonで扱うには、ライブラリ pandas が提供するオブジェクト DataFrame を利用するのが最も効率的です。pandas については、3.1 章のPython スクリプトに "時刻処理に便利なライブラリ" とコメントを付けましたが、それは付随的な機能であり、本領は表形式のデータを処理することです。そして、オブジェクト DataFrame は、表形式のデータの入れ物として提供するオブジェクトです。表計算ソフトExcelでいうところのワークシートに相当するものです。
pandas についての完全な説明はホームページ (https://pandas.pydata.org/) にあります。とても広く使われているライブラリなので、ネット上にはこれ以外にも多くの解説記事が存在します。
今回の表を格納する DataFrame オブジェクトを df とすれば、 df は、もっとも簡単には、下に示すようにして作成できます。この例でわかるように、オブジェクト DataFrame は、埋め込むデータがまだない段階でも作ることができます。このような場合、表の内容には「NaN」(not a number) が埋め込まれます。
import pandas as pd
idx = ['2023-06-02 00:00','2023-06-02 01:00','2023-06-02 02:00']
cols = ['突合データ','GPV抽出データ']
df = pd.DataFrame(data=None, index=idx, columns=cols)
多くの場合、突き合わせるべきデータはCSVファイルで提供されます。ライブラリpandasが提供する関数 readcsv を使うと、CSVファイルから直接 Dataframe オブジェクトを生成することができます。これを使って、気象庁からダウンロードしたアメダスデータの DataFrame オブジェクトを作成してみましょう。
使用するサンプルは、アメダス「舘野(つくば)」における、2023年6月1日1時から6月16日00時までの時別気温のCSVファイルです。気象庁のホームページ( https://www.data.jma.go.jp/gmd/risk/obsdl/ )からダウンロードしたものです。フォルダ challenge4/jmadata に、data.csvとして置いています。内容は以下のようなものです。
CSVは、ファイルフォーマットの観点からはかなりいい加減なので、これに対応するため、メソッド read_csv には数えきれないほど多数のオプションがあります。上記のファイルを読むために、今回は以下を使用します;
- キーワード引数 encoding :文字コードを指定する
- キーワード引数 skiprows :データの上から5行目までを読みとばす
- キーワード引数 names :列見出しを付け直す
- キーワード引数 index_col :最初の列をインデックスに指定する
- キーワード引数 parse_dates :日付の文字列を 時刻と解釈する
以下を実行してください、CSVファイルのアメダスデータが、Dataframe オブジェクトの形式で df_xdata に格納されます。
df_xdata = pd.read_csv( "./jmadata/data.csv",
encoding='Shift_JIS',
skiprows=5,
names=('Ta_amd','品質番号','均質番号'),
index_col=0,
parse_dates=True
)
結果を確認してみましょう、以下を実行して、概要を表示させてください。
df_xdata
Ta_amd | 品質番号 | 均質番号 | |
---|---|---|---|
2023-06-01 01:00:00 | 12.7 | 8 | 1 |
2023-06-01 02:00:00 | 12.3 | 8 | 1 |
2023-06-01 03:00:00 | 11.8 | 8 | 1 |
2023-06-01 04:00:00 | 11.7 | 8 | 1 |
2023-06-01 05:00:00 | 11.5 | 8 | 1 |
... | ... | ... | ... |
2023-06-15 20:00:00 | 20.2 | 8 | 1 |
2023-06-15 21:00:00 | 20.2 | 8 | 1 |
2023-06-15 22:00:00 | 20.0 | 8 | 1 |
2023-06-15 23:00:00 | 20.0 | 8 | 1 |
2023-06-16 00:00:00 | 19.8 | 8 | 1 |
360 rows × 3 columns
2023年6月1日1時から1月16日0時までの気温データが保存されているのがわかります。
さて、品質番号と均質番号は今回は不要なので、列ごとを削除しましょう。それにはメソッド drop を使用します。全般的に言えることですが、Dataframe のメソッドは、列の挿入や削除など自身に変更を加えるものの場合、本当には変更せず、"変更するとこうなりますケド" という姿(view と呼びます)を返します。従って、本当に変更してしまう場合は、結果を自分自身に代入します。
df_xdata = df_xdata.drop(columns=['品質番号','均質番号'])
df_xdata
Ta_amd | |
---|---|
2023-06-01 01:00:00 | 12.7 |
2023-06-01 02:00:00 | 12.3 |
2023-06-01 03:00:00 | 11.8 |
2023-06-01 04:00:00 | 11.7 |
2023-06-01 05:00:00 | 11.5 |
... | ... |
2023-06-15 20:00:00 | 20.2 |
2023-06-15 21:00:00 | 20.2 |
2023-06-15 22:00:00 | 20.0 |
2023-06-15 23:00:00 | 20.0 |
2023-06-16 00:00:00 | 19.8 |
360 rows × 1 columns
これで、時刻と気温データだけのオブジェクトになりました。
3.2.2 データの突き合わせ¶
日本標準時のデータ同士の突き合わせであれば、3.2.1 のオブジェクトでOKですが、今回はUTCのデータと突き合わせるのでもうひと手間必要です。以下のようにして、時刻を awareな日本標準時に変更します(時刻が日本標準時表記であることを認識させます)。
df_xdata.index = df_xdata.index.tz_localize('Asia/Tokyo')
すると、naive だったインデックスの時刻 aware になり、表示も以下のようになります。
実際に変更してみましょう。以下を実行してください。
df_xdata.index = df_xdata.index.tz_localize('Asia/Tokyo') # お前は日本標準時だと言い聞かせる
これで突合するデータができました。確認したい人は以下を実行してください。
df_xdata
Ta_amd | |
---|---|
2023-06-01 01:00:00+09:00 | 12.7 |
2023-06-01 02:00:00+09:00 | 12.3 |
2023-06-01 03:00:00+09:00 | 11.8 |
2023-06-01 04:00:00+09:00 | 11.7 |
2023-06-01 05:00:00+09:00 | 11.5 |
... | ... |
2023-06-15 20:00:00+09:00 | 20.2 |
2023-06-15 21:00:00+09:00 | 20.2 |
2023-06-15 22:00:00+09:00 | 20.0 |
2023-06-15 23:00:00+09:00 | 20.0 |
2023-06-16 00:00:00+09:00 | 19.8 |
360 rows × 1 columns
次はGPVデータの取り出しです。
アメダス「舘野(つくば)」に最も近い格子点の気温を取り出します。ついでに、単位の変更(K → Celsius)もしておきましょう。
lat, lon = 36.05666667, 140.125 #「舘野(つくば)」
da_tateno = da.sel(latitude=lat, longitude=lon, method='nearest') - 273.15
da_tateno として取り出されたのは、時間の次元をもつ1次元の DataArray オブジェクトです。これから、 DataFrame オブジェクト df_gpv を作ります。index(行見出し)値に座標 awaretime を使用するのがミソです。
df_gpv = pd.DataFrame(data=da_tateno.data ,index=da_tateno['awaretime'].data, columns=['Ta_gpv'])
df_gpv
Ta_gpv | |
---|---|
2023-06-01 09:00:00+09:00 | 20.314270 |
2023-06-01 10:00:00+09:00 | 21.548218 |
2023-06-01 11:00:00+09:00 | 22.487213 |
2023-06-01 12:00:00+09:00 | 23.248871 |
2023-06-01 13:00:00+09:00 | 23.517914 |
... | ... |
2023-06-11 21:00:00+09:00 | 21.579315 |
2023-06-12 00:00:00+09:00 | 20.724731 |
2023-06-12 03:00:00+09:00 | 19.998657 |
2023-06-12 06:00:00+09:00 | 20.386078 |
2023-06-12 09:00:00+09:00 | 24.375092 |
177 rows × 1 columns
これで、GSM-GPV から取り出した気温の DataFrame ができました。
それでは両者を結合します。結合には、ライブラリ pandas が提供する関数 merge を使用して以下のようにします。
df = pd.merge(df_xdata, df_gpv, # 結合するDataFrame
how='outer', # どちらかにデータがあれば結果に含めるので'outer'
left_index=True, # 左のDataFrameのインデックスをキーとして使うからTrue
right_index=True # 右のDataFrameのインデックスもキーとして使うからTrue
)
df
Ta_amd | Ta_gpv | |
---|---|---|
2023-06-01 01:00:00+09:00 | 12.7 | NaN |
2023-06-01 02:00:00+09:00 | 12.3 | NaN |
2023-06-01 03:00:00+09:00 | 11.8 | NaN |
2023-06-01 04:00:00+09:00 | 11.7 | NaN |
2023-06-01 05:00:00+09:00 | 11.5 | NaN |
... | ... | ... |
2023-06-15 20:00:00+09:00 | 20.2 | NaN |
2023-06-15 21:00:00+09:00 | 20.2 | NaN |
2023-06-15 22:00:00+09:00 | 20.0 | NaN |
2023-06-15 23:00:00+09:00 | 20.0 | NaN |
2023-06-16 00:00:00+09:00 | 19.8 | NaN |
360 rows × 2 columns
これで完成しました。GPVは、設定した期間の最初の部分と最後の部分のデータが無いので、列全体が NaN であるかのように見えていますが、きちんとマージされているので安心してください。折れ線グラフにして確認してみましょう。pandas にも描画のメソッドがあるので、これを使います。メソッドの名前はxarrayと同じ、plot です。
df.plot(figsize=(10,3), marker='o')
<AxesSubplot: >
今回は、突き合わせるデータの時刻もGPVの時刻も両方 aware な日本標準時に揃えてマージしましたが、Pandas の関数 merge は、時刻がawareでさえあれば、仮にタイムゾーンが一致していなくても、時刻をキーとするマージを正しく行います。
3.2.3 タイムゾーンの変更¶
必要となる機会は多くないと思いますが、timestamp オブジェクトのタイムゾーンの変更の仕方を示しておきます。
現在、オブジェクト df は、日本標準時で整理されています。これを中央協定時で整理し直す場合には、インデックスにメソッド tz_convert を適用し、引数に文字列 'UTC' を与えます。
以下を実行して、時刻が中央協定時で表現されていることを確認してください。
dfUTC = df.copy() # dfのコピーを作りこれを弄る
dfUTC.index = df.index.tz_convert('UTC') # インデックスをUTC表示にしてくれと頼む
dfUTC
Ta_amd | Ta_gpv | |
---|---|---|
2023-05-31 16:00:00+00:00 | 12.7 | NaN |
2023-05-31 17:00:00+00:00 | 12.3 | NaN |
2023-05-31 18:00:00+00:00 | 11.8 | NaN |
2023-05-31 19:00:00+00:00 | 11.7 | NaN |
2023-05-31 20:00:00+00:00 | 11.5 | NaN |
... | ... | ... |
2023-06-15 11:00:00+00:00 | 20.2 | NaN |
2023-06-15 12:00:00+00:00 | 20.2 | NaN |
2023-06-15 13:00:00+00:00 | 20.0 | NaN |
2023-06-15 14:00:00+00:00 | 20.0 | NaN |
2023-06-15 15:00:00+00:00 | 19.8 | NaN |
360 rows × 2 columns
4 補間¶
第 4 章 を始めるにあたり、Python カーネルをリスタート してメモリをきれいにしておきましょう。 ツールバーにある輪のようなアイコンをクリックし、確認に対して [Restart] をクリックしてください。
その後に、下記を実行して必要となるライブラリーをインポートとしてください。
import numpy as np
import xarray as xr
import pandas as pd
import wxbcgribx as wx
補間とは、格子点が存在ない時空間位置における値を周囲の格子点値から推定することです。
気象値のGPVを扱う過程では、格子点間隔が大きなGPVからより小さな格子間隔のGPVを推定する場合によく使用されますが、補間は、あくまで任意の場所の値の推定ですから、細かな格子間隔のGPVから疎らに値を取り出して格子間隔が大きなGPVを作ることも、処理としては可能です。
DataArray オブジェクトには、データを補間するメソッドとして、interp と、interp_like が用意されており、これらを用いて GPV を簡単に補間することができます。
4.1 メソッド interp¶
メソッド interp は、DataArrayオブジェクトを補間する基本的なメソッドです。メソッド interp は次のようにして用います。
da_new = da.interp( 座標名 = 補間したい座標値の配列, method = "補間方法")
ここで、da は補間をしたい GPV を格納する DataArray オブジェクト、da_new は補間された GPV を格納する DataArray オブジェクトです。座標名 は、補間をする座標の名前(time, longitude, latitude のいずれか1つ)です。
それでは、全球数値予報モデルGPV (GSM-GPV) を素材に、これを実際に使ってみましょう。
まずは、サンプルのGRIBファイルから2023年6月1日0時UTCを初期値とする気温予報のGSM-GPVをあらためて読み込み、DataArrayオブジェクト da_gsm に格納しておきます。
yyyymmdd = '20230601' # 読み込むGPVの初期値の日付
prdct = 'gsm'
fd = ['0000-0100','0101-0200','0201-0300','0301-0400','0401-0500','0501-0512',
'0515-0700','0703-0900','0903-1100']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_GSM_GPV_Rjp_Gll0p1deg_Lsurf_FD{xx}_grib2.bin'
for xx in fd]
ds = wx.getgpv(grb_paths, 'TMP', timezone='Asia/Tokyo')
da_gsm = ds['TMP_2maboveground']
GSM-GPVの予報は、初期日から132時間(5.5日)先までについては1時間毎ですが、それから先は3時間毎であり一様ではありません。以下を実行すると、時刻の座標値が表示されるので、そのことを確認しましょう。
2023年06月06日12時00分(UTC)までは1時間刻みで、その次は3時間刻みとなっていることが確認できるでしょうか。
print(da_gsm.time.dt.strftime('%Y-%m-%dT%H:%M').data)
['2023-06-01T00:00' '2023-06-01T01:00' '2023-06-01T02:00' '2023-06-01T03:00' '2023-06-01T04:00' '2023-06-01T05:00' '2023-06-01T06:00' '2023-06-01T07:00' '2023-06-01T08:00' '2023-06-01T09:00' '2023-06-01T10:00' '2023-06-01T11:00' '2023-06-01T12:00' '2023-06-01T13:00' '2023-06-01T14:00' '2023-06-01T15:00' '2023-06-01T16:00' '2023-06-01T17:00' '2023-06-01T18:00' '2023-06-01T19:00' '2023-06-01T20:00' '2023-06-01T21:00' '2023-06-01T22:00' '2023-06-01T23:00' '2023-06-02T00:00' '2023-06-02T01:00' '2023-06-02T02:00' '2023-06-02T03:00' '2023-06-02T04:00' '2023-06-02T05:00' '2023-06-02T06:00' '2023-06-02T07:00' '2023-06-02T08:00' '2023-06-02T09:00' '2023-06-02T10:00' '2023-06-02T11:00' '2023-06-02T12:00' '2023-06-02T13:00' '2023-06-02T14:00' '2023-06-02T15:00' '2023-06-02T16:00' '2023-06-02T17:00' '2023-06-02T18:00' '2023-06-02T19:00' '2023-06-02T20:00' '2023-06-02T21:00' '2023-06-02T22:00' '2023-06-02T23:00' '2023-06-03T00:00' '2023-06-03T01:00' '2023-06-03T02:00' '2023-06-03T03:00' '2023-06-03T04:00' '2023-06-03T05:00' '2023-06-03T06:00' '2023-06-03T07:00' '2023-06-03T08:00' '2023-06-03T09:00' '2023-06-03T10:00' '2023-06-03T11:00' '2023-06-03T12:00' '2023-06-03T13:00' '2023-06-03T14:00' '2023-06-03T15:00' '2023-06-03T16:00' '2023-06-03T17:00' '2023-06-03T18:00' '2023-06-03T19:00' '2023-06-03T20:00' '2023-06-03T21:00' '2023-06-03T22:00' '2023-06-03T23:00' '2023-06-04T00:00' '2023-06-04T01:00' '2023-06-04T02:00' '2023-06-04T03:00' '2023-06-04T04:00' '2023-06-04T05:00' '2023-06-04T06:00' '2023-06-04T07:00' '2023-06-04T08:00' '2023-06-04T09:00' '2023-06-04T10:00' '2023-06-04T11:00' '2023-06-04T12:00' '2023-06-04T13:00' '2023-06-04T14:00' '2023-06-04T15:00' '2023-06-04T16:00' '2023-06-04T17:00' '2023-06-04T18:00' '2023-06-04T19:00' '2023-06-04T20:00' '2023-06-04T21:00' '2023-06-04T22:00' '2023-06-04T23:00' '2023-06-05T00:00' '2023-06-05T01:00' '2023-06-05T02:00' '2023-06-05T03:00' '2023-06-05T04:00' '2023-06-05T05:00' '2023-06-05T06:00' '2023-06-05T07:00' '2023-06-05T08:00' '2023-06-05T09:00' '2023-06-05T10:00' '2023-06-05T11:00' '2023-06-05T12:00' '2023-06-05T13:00' '2023-06-05T14:00' '2023-06-05T15:00' '2023-06-05T16:00' '2023-06-05T17:00' '2023-06-05T18:00' '2023-06-05T19:00' '2023-06-05T20:00' '2023-06-05T21:00' '2023-06-05T22:00' '2023-06-05T23:00' '2023-06-06T00:00' '2023-06-06T01:00' '2023-06-06T02:00' '2023-06-06T03:00' '2023-06-06T04:00' '2023-06-06T05:00' '2023-06-06T06:00' '2023-06-06T07:00' '2023-06-06T08:00' '2023-06-06T09:00' '2023-06-06T10:00' '2023-06-06T11:00' '2023-06-06T12:00' '2023-06-06T15:00' '2023-06-06T18:00' '2023-06-06T21:00' '2023-06-07T00:00' '2023-06-07T03:00' '2023-06-07T06:00' '2023-06-07T09:00' '2023-06-07T12:00' '2023-06-07T15:00' '2023-06-07T18:00' '2023-06-07T21:00' '2023-06-08T00:00' '2023-06-08T03:00' '2023-06-08T06:00' '2023-06-08T09:00' '2023-06-08T12:00' '2023-06-08T15:00' '2023-06-08T18:00' '2023-06-08T21:00' '2023-06-09T00:00' '2023-06-09T03:00' '2023-06-09T06:00' '2023-06-09T09:00' '2023-06-09T12:00' '2023-06-09T15:00' '2023-06-09T18:00' '2023-06-09T21:00' '2023-06-10T00:00' '2023-06-10T03:00' '2023-06-10T06:00' '2023-06-10T09:00' '2023-06-10T12:00' '2023-06-10T15:00' '2023-06-10T18:00' '2023-06-10T21:00' '2023-06-11T00:00' '2023-06-11T03:00' '2023-06-11T06:00' '2023-06-11T09:00' '2023-06-11T12:00' '2023-06-11T15:00' '2023-06-11T18:00' '2023-06-11T21:00' '2023-06-12T00:00']
北緯35度、東経137度の格子点値をグラフに表すと以下のようになります。
このデータを時間について補間し、全期間にわたって間隔が1時間である GPV を作りましょう。
そのためには、全期間にわたって間隔が1時間である時刻の配列を用意する必要があります。ここでは、2つの方法を紹介します。どちらでも大丈夫です。
# その1 ライブラリpandasの 関数date_rangeを使う方法
beg = pd.to_datetime(da_gsm.time[0].data) # GPVの最初の時刻を取り出しbegに代入
end = pd.to_datetime(da_gsm.time[-1].data) # GPVの最後の時刻を取り出しbegに代入
time1h = pd.date_range(beg, end, freq="H" )
time1h
DatetimeIndex(['2023-06-01 00:00:00', '2023-06-01 01:00:00', '2023-06-01 02:00:00', '2023-06-01 03:00:00', '2023-06-01 04:00:00', '2023-06-01 05:00:00', '2023-06-01 06:00:00', '2023-06-01 07:00:00', '2023-06-01 08:00:00', '2023-06-01 09:00:00', ... '2023-06-11 15:00:00', '2023-06-11 16:00:00', '2023-06-11 17:00:00', '2023-06-11 18:00:00', '2023-06-11 19:00:00', '2023-06-11 20:00:00', '2023-06-11 21:00:00', '2023-06-11 22:00:00', '2023-06-11 23:00:00', '2023-06-12 00:00:00'], dtype='datetime64[ns]', length=265, freq='H')
(結果のセルの左側にあるブルーの棒線をクリックすると表示を畳むことができます。邪魔なようなら畳んでください)
# その2 ライブラリnumpyの 関数arangeを使う方法
oh = np.timedelta64(1,'h')
time1h = np.arange(da_gsm.time[0].data, da_gsm.time[-1].data+oh, oh)
time1h
array(['2023-06-01T00:00:00.000000000', '2023-06-01T01:00:00.000000000', '2023-06-01T02:00:00.000000000', '2023-06-01T03:00:00.000000000', '2023-06-01T04:00:00.000000000', '2023-06-01T05:00:00.000000000', '2023-06-01T06:00:00.000000000', '2023-06-01T07:00:00.000000000', '2023-06-01T08:00:00.000000000', '2023-06-01T09:00:00.000000000', '2023-06-01T10:00:00.000000000', '2023-06-01T11:00:00.000000000', '2023-06-01T12:00:00.000000000', '2023-06-01T13:00:00.000000000', '2023-06-01T14:00:00.000000000', '2023-06-01T15:00:00.000000000', '2023-06-01T16:00:00.000000000', '2023-06-01T17:00:00.000000000', '2023-06-01T18:00:00.000000000', '2023-06-01T19:00:00.000000000', '2023-06-01T20:00:00.000000000', '2023-06-01T21:00:00.000000000', '2023-06-01T22:00:00.000000000', '2023-06-01T23:00:00.000000000', '2023-06-02T00:00:00.000000000', '2023-06-02T01:00:00.000000000', '2023-06-02T02:00:00.000000000', '2023-06-02T03:00:00.000000000', '2023-06-02T04:00:00.000000000', '2023-06-02T05:00:00.000000000', '2023-06-02T06:00:00.000000000', '2023-06-02T07:00:00.000000000', '2023-06-02T08:00:00.000000000', '2023-06-02T09:00:00.000000000', '2023-06-02T10:00:00.000000000', '2023-06-02T11:00:00.000000000', '2023-06-02T12:00:00.000000000', '2023-06-02T13:00:00.000000000', '2023-06-02T14:00:00.000000000', '2023-06-02T15:00:00.000000000', '2023-06-02T16:00:00.000000000', '2023-06-02T17:00:00.000000000', '2023-06-02T18:00:00.000000000', '2023-06-02T19:00:00.000000000', '2023-06-02T20:00:00.000000000', '2023-06-02T21:00:00.000000000', '2023-06-02T22:00:00.000000000', '2023-06-02T23:00:00.000000000', '2023-06-03T00:00:00.000000000', '2023-06-03T01:00:00.000000000', '2023-06-03T02:00:00.000000000', '2023-06-03T03:00:00.000000000', '2023-06-03T04:00:00.000000000', '2023-06-03T05:00:00.000000000', '2023-06-03T06:00:00.000000000', '2023-06-03T07:00:00.000000000', '2023-06-03T08:00:00.000000000', '2023-06-03T09:00:00.000000000', '2023-06-03T10:00:00.000000000', '2023-06-03T11:00:00.000000000', '2023-06-03T12:00:00.000000000', '2023-06-03T13:00:00.000000000', '2023-06-03T14:00:00.000000000', '2023-06-03T15:00:00.000000000', '2023-06-03T16:00:00.000000000', '2023-06-03T17:00:00.000000000', '2023-06-03T18:00:00.000000000', '2023-06-03T19:00:00.000000000', '2023-06-03T20:00:00.000000000', '2023-06-03T21:00:00.000000000', '2023-06-03T22:00:00.000000000', '2023-06-03T23:00:00.000000000', '2023-06-04T00:00:00.000000000', '2023-06-04T01:00:00.000000000', '2023-06-04T02:00:00.000000000', '2023-06-04T03:00:00.000000000', '2023-06-04T04:00:00.000000000', '2023-06-04T05:00:00.000000000', '2023-06-04T06:00:00.000000000', '2023-06-04T07:00:00.000000000', '2023-06-04T08:00:00.000000000', '2023-06-04T09:00:00.000000000', '2023-06-04T10:00:00.000000000', '2023-06-04T11:00:00.000000000', '2023-06-04T12:00:00.000000000', '2023-06-04T13:00:00.000000000', '2023-06-04T14:00:00.000000000', '2023-06-04T15:00:00.000000000', '2023-06-04T16:00:00.000000000', '2023-06-04T17:00:00.000000000', '2023-06-04T18:00:00.000000000', '2023-06-04T19:00:00.000000000', '2023-06-04T20:00:00.000000000', '2023-06-04T21:00:00.000000000', '2023-06-04T22:00:00.000000000', '2023-06-04T23:00:00.000000000', '2023-06-05T00:00:00.000000000', '2023-06-05T01:00:00.000000000', '2023-06-05T02:00:00.000000000', '2023-06-05T03:00:00.000000000', '2023-06-05T04:00:00.000000000', '2023-06-05T05:00:00.000000000', '2023-06-05T06:00:00.000000000', '2023-06-05T07:00:00.000000000', '2023-06-05T08:00:00.000000000', '2023-06-05T09:00:00.000000000', '2023-06-05T10:00:00.000000000', '2023-06-05T11:00:00.000000000', '2023-06-05T12:00:00.000000000', '2023-06-05T13:00:00.000000000', '2023-06-05T14:00:00.000000000', '2023-06-05T15:00:00.000000000', '2023-06-05T16:00:00.000000000', '2023-06-05T17:00:00.000000000', '2023-06-05T18:00:00.000000000', '2023-06-05T19:00:00.000000000', '2023-06-05T20:00:00.000000000', '2023-06-05T21:00:00.000000000', '2023-06-05T22:00:00.000000000', '2023-06-05T23:00:00.000000000', '2023-06-06T00:00:00.000000000', '2023-06-06T01:00:00.000000000', '2023-06-06T02:00:00.000000000', '2023-06-06T03:00:00.000000000', '2023-06-06T04:00:00.000000000', '2023-06-06T05:00:00.000000000', '2023-06-06T06:00:00.000000000', '2023-06-06T07:00:00.000000000', '2023-06-06T08:00:00.000000000', '2023-06-06T09:00:00.000000000', '2023-06-06T10:00:00.000000000', '2023-06-06T11:00:00.000000000', '2023-06-06T12:00:00.000000000', '2023-06-06T13:00:00.000000000', '2023-06-06T14:00:00.000000000', '2023-06-06T15:00:00.000000000', '2023-06-06T16:00:00.000000000', '2023-06-06T17:00:00.000000000', '2023-06-06T18:00:00.000000000', '2023-06-06T19:00:00.000000000', '2023-06-06T20:00:00.000000000', '2023-06-06T21:00:00.000000000', '2023-06-06T22:00:00.000000000', '2023-06-06T23:00:00.000000000', '2023-06-07T00:00:00.000000000', '2023-06-07T01:00:00.000000000', '2023-06-07T02:00:00.000000000', '2023-06-07T03:00:00.000000000', '2023-06-07T04:00:00.000000000', '2023-06-07T05:00:00.000000000', '2023-06-07T06:00:00.000000000', '2023-06-07T07:00:00.000000000', '2023-06-07T08:00:00.000000000', '2023-06-07T09:00:00.000000000', '2023-06-07T10:00:00.000000000', '2023-06-07T11:00:00.000000000', '2023-06-07T12:00:00.000000000', '2023-06-07T13:00:00.000000000', '2023-06-07T14:00:00.000000000', '2023-06-07T15:00:00.000000000', '2023-06-07T16:00:00.000000000', '2023-06-07T17:00:00.000000000', '2023-06-07T18:00:00.000000000', '2023-06-07T19:00:00.000000000', '2023-06-07T20:00:00.000000000', '2023-06-07T21:00:00.000000000', '2023-06-07T22:00:00.000000000', '2023-06-07T23:00:00.000000000', '2023-06-08T00:00:00.000000000', '2023-06-08T01:00:00.000000000', '2023-06-08T02:00:00.000000000', '2023-06-08T03:00:00.000000000', '2023-06-08T04:00:00.000000000', '2023-06-08T05:00:00.000000000', '2023-06-08T06:00:00.000000000', '2023-06-08T07:00:00.000000000', '2023-06-08T08:00:00.000000000', '2023-06-08T09:00:00.000000000', '2023-06-08T10:00:00.000000000', '2023-06-08T11:00:00.000000000', '2023-06-08T12:00:00.000000000', '2023-06-08T13:00:00.000000000', '2023-06-08T14:00:00.000000000', '2023-06-08T15:00:00.000000000', '2023-06-08T16:00:00.000000000', '2023-06-08T17:00:00.000000000', '2023-06-08T18:00:00.000000000', '2023-06-08T19:00:00.000000000', '2023-06-08T20:00:00.000000000', '2023-06-08T21:00:00.000000000', '2023-06-08T22:00:00.000000000', '2023-06-08T23:00:00.000000000', '2023-06-09T00:00:00.000000000', '2023-06-09T01:00:00.000000000', '2023-06-09T02:00:00.000000000', '2023-06-09T03:00:00.000000000', '2023-06-09T04:00:00.000000000', '2023-06-09T05:00:00.000000000', '2023-06-09T06:00:00.000000000', '2023-06-09T07:00:00.000000000', '2023-06-09T08:00:00.000000000', '2023-06-09T09:00:00.000000000', '2023-06-09T10:00:00.000000000', '2023-06-09T11:00:00.000000000', '2023-06-09T12:00:00.000000000', '2023-06-09T13:00:00.000000000', '2023-06-09T14:00:00.000000000', '2023-06-09T15:00:00.000000000', '2023-06-09T16:00:00.000000000', '2023-06-09T17:00:00.000000000', '2023-06-09T18:00:00.000000000', '2023-06-09T19:00:00.000000000', '2023-06-09T20:00:00.000000000', '2023-06-09T21:00:00.000000000', '2023-06-09T22:00:00.000000000', '2023-06-09T23:00:00.000000000', '2023-06-10T00:00:00.000000000', '2023-06-10T01:00:00.000000000', '2023-06-10T02:00:00.000000000', '2023-06-10T03:00:00.000000000', '2023-06-10T04:00:00.000000000', '2023-06-10T05:00:00.000000000', '2023-06-10T06:00:00.000000000', '2023-06-10T07:00:00.000000000', '2023-06-10T08:00:00.000000000', '2023-06-10T09:00:00.000000000', '2023-06-10T10:00:00.000000000', '2023-06-10T11:00:00.000000000', '2023-06-10T12:00:00.000000000', '2023-06-10T13:00:00.000000000', '2023-06-10T14:00:00.000000000', '2023-06-10T15:00:00.000000000', '2023-06-10T16:00:00.000000000', '2023-06-10T17:00:00.000000000', '2023-06-10T18:00:00.000000000', '2023-06-10T19:00:00.000000000', '2023-06-10T20:00:00.000000000', '2023-06-10T21:00:00.000000000', '2023-06-10T22:00:00.000000000', '2023-06-10T23:00:00.000000000', '2023-06-11T00:00:00.000000000', '2023-06-11T01:00:00.000000000', '2023-06-11T02:00:00.000000000', '2023-06-11T03:00:00.000000000', '2023-06-11T04:00:00.000000000', '2023-06-11T05:00:00.000000000', '2023-06-11T06:00:00.000000000', '2023-06-11T07:00:00.000000000', '2023-06-11T08:00:00.000000000', '2023-06-11T09:00:00.000000000', '2023-06-11T10:00:00.000000000', '2023-06-11T11:00:00.000000000', '2023-06-11T12:00:00.000000000', '2023-06-11T13:00:00.000000000', '2023-06-11T14:00:00.000000000', '2023-06-11T15:00:00.000000000', '2023-06-11T16:00:00.000000000', '2023-06-11T17:00:00.000000000', '2023-06-11T18:00:00.000000000', '2023-06-11T19:00:00.000000000', '2023-06-11T20:00:00.000000000', '2023-06-11T21:00:00.000000000', '2023-06-11T22:00:00.000000000', '2023-06-11T23:00:00.000000000', '2023-06-12T00:00:00.000000000'], dtype='datetime64[ns]')
これを使って、メソッド interp で時間についての補間をします。以下を実行してください。
# 補間
da1 = da_gsm.interp(time=time1h, method='cubic')
補間後のGPVを、折れ線グラフで見てみましょう。以下を実行してください。
lat,lon = 35.2, 137.0
time = da1.time
var = da1.sel(latitude=lat, longitude=lon)
wx.ts(time,var,
title = f'(N{lat:0.2f}, E{lon:0.2f})',
xlabel = f"Time (UTC)",
ylabel = f"{da1.attrs['long_name']} [{da1.attrs['units']}]",
llabel = 'interpolated',
)
全期間にわたって、1時間おきのデータとなっていることが確認されました。
今実行したのは時間についての補間ですが、同様の方法で、緯度、経度についても補完することができます。また、今回は、補間方法として、cubic を採用しましたが、それ以外に 'nearest' と 'linear' を指定できます。どの補間が適切かはケースバイケースで検討してください。
4.2 メソッド interp_like¶
メソッド interp_like は、既存の DataArray の座標とぴったりと合うように補間をするメソッドです。同じ処理はメソッド interp を複数回適用しても可能で、このメソッドは、近道、と理解してください。
メソッド interp_like は以下のようにして用います。
da_new = da.interp_like(da_ref)
ここで、da は補間をかけたい GPV の DataArray オブジェクト、da_ref は格子を一致させたい GPV の DataArray オブジェクトを示します。
それでは、空間範囲を日本の中部域、時間範囲を2023年6月1日0時UTCから6月04日06時UTCまでとする MSM-GPV を下敷きとして、次元と座標がこれに一致するように GSM-GPV を補完してみましょう。
緯度経度格子の違いを予め把握しておくため、二つのGPVの分布図を並べて示します。
まず、MSM-GPVを、領域を中部域に限って読み込んで DataArray オブジェクト da_msm に格納します。
yyyymmdd = '20230601'
prdct = 'msm'
fh = ['00-15', '16-33', '34-39', '40-51', '52-78']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_MSM_GPV_Rjp_Lsurf_FH{xx}_grib2.bin'
for xx in fh]
area = [34, 36, 135, 138] # GPVを読み込む領域
ds = wx.getgpv(grb_paths, 'TMP', timezone='Asia/Tokyo', lalomima=area)
da_msm = ds['TMP_1D5maboveground']
次に概要リストを表示させて、時刻、緯度、経度、座標を確認できるようにしておきましょう。データ本体は畳んでおいてください。
da_msm
<xarray.DataArray 'TMP_1D5maboveground' (time: 79, latitude: 39, longitude: 47)> array([[[294.25366, 295.1521 , 296.07397, ..., 295.6521 , 295.6443 , 295.63647], [293.99585, 295.67554, 296.04272, ..., 295.60522, 295.60522, 295.61304], [293.6521 , 294.23804, 296.32397, ..., 295.5349 , 295.54272, 295.54272], ..., [293.87866, 293.51147, 293.55835, ..., 287.07397, 288.81616, 291.24585], [293.5974 , 293.23022, 292.88647, ..., 287.85522, 288.42554, 290.51147], [293.26147, 292.7068 , 292.2693 , ..., 287.9568 , 288.2849 , 289.56616]], [[294.46075, 295.6795 , 297.242 , ..., 295.6717 , 295.69513, 295.68732], [294.16388, 296.21857, 296.96857, ..., 295.63263, 295.63263, 295.6092 ], [293.84357, 294.4217 , 296.95294, ..., 295.6092 , 295.58575, 295.5467 ], ... [293.25098, 293.18848, 293.18848, ..., 290.53223, 292.40723, 296.29785], [293.17285, 293.1416 , 293.09473, ..., 291.90723, 292.40723, 295.71973], [293.09473, 293.06348, 293.00098, ..., 292.2666 , 292.40723, 294.3291 ]], [[294.4373 , 295.8748 , 297.78104, ..., 295.8748 , 295.84354, 295.8279 ], [294.2498 , 296.5154 , 297.4998 , ..., 295.8279 , 295.8123 , 295.78104], [293.9529 , 294.9529 , 297.96854, ..., 295.78104, 295.7654 , 295.73416], ..., [293.2654 , 293.2029 , 293.15604, ..., 290.54666, 292.34354, 296.2029 ], [293.2498 , 293.15604, 293.10916, ..., 291.85916, 292.3123 , 295.5154 ], [293.21854, 293.1248 , 293.03104, ..., 292.2654 , 292.3748 , 294.21854]]], dtype=float32) Coordinates: * latitude (latitude) float64 34.05 34.1 34.15 34.2 ... 35.85 35.9 35.95 * longitude (longitude) float64 135.1 135.1 135.2 135.2 ... 137.8 137.9 137.9 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-04T06:00:00 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-04T15:00:0... Attributes: short_name: TMP_1D5maboveground long_name: Temperature level: 1.5 m above ground units: K
では、メソッド interp_like で da_gsm (GSM-GPV)を補間して da_msm と同じ時空間範囲の GPV を作り、da_gsmi に格納しましょう。
da_gsmi = da_gsm.interp_like(da_msm, method='linear')
以下を実行して概要フォームを表示させ、時空間範囲を確認してください。GSM-GPVの予報期間は132時間先までですが、補間後は78時間先で切られています。緯度経度についても、もともとの領域は日本列島がすっぽり入る範囲でしたが、領域を限定したMSM-GPVの範囲でトリミングされていることが分かります。
da_gsmi
<xarray.DataArray 'TMP_2maboveground' (time: 79, latitude: 39, longitude: 47)> array([[[293.25527954, 294.01504517, 294.39004517, ..., 294.28652954, 294.26113892, 294.24356079], [293.05020142, 293.80020142, 294.21426392, ..., 294.17910767, 294.15176392, 294.13613892], [292.89590454, 293.57754517, 294.00527954, ..., 294.10098267, 294.07363892, 294.05801392], ..., [292.69668579, 292.55020142, 292.39395142, ..., 289.30801392, 289.84317017, 290.12246704], [292.14785767, 292.01895142, 291.87051392, ..., 289.54629517, 290.12832642, 290.44082642], [291.61465454, 291.50332642, 291.37832642, ..., 289.68692017, 290.28067017, 290.56973267]], [[295.05172729, 295.72555542, 295.91500854, ..., 294.50094604, 294.47555542, 294.46774292], [294.93258667, 295.62008667, 295.88571167, ..., 294.34274292, 294.31539917, 294.30758667], [294.77047729, 295.40524292, 295.73922729, ..., 294.19430542, 294.16696167, 294.15719604], ... [293.42828369, 293.37164307, 293.24664307, ..., 293.71929932, 294.32476807, 294.35406494], [293.29742432, 293.23883057, 293.15679932, ..., 293.69976807, 294.44195557, 294.61773682], [293.28570557, 293.22711182, 293.15679932, ..., 293.57672119, 294.40679932, 294.65289307]], [[295.83178711, 296.46264648, 296.75366211, ..., 295.95678711, 295.91186523, 295.86889648], [295.89624023, 296.48999023, 296.76342773, ..., 295.88452148, 295.83374023, 295.79077148], [295.98413086, 296.50561523, 296.76342773, ..., 295.79077148, 295.73999023, 295.70483398], ..., [293.74975586, 293.72827148, 293.61499023, ..., 293.50170898, 294.13842773, 294.21264648], [293.58374023, 293.56811523, 293.50561523, ..., 293.57983398, 294.35717773, 294.55639648], [293.52905273, 293.51342773, 293.46850586, ..., 293.43725586, 294.29077148, 294.54272461]]]) Coordinates: * latitude (latitude) float64 34.05 34.1 34.15 34.2 ... 35.85 35.9 35.95 * longitude (longitude) float64 135.1 135.1 135.2 135.2 ... 137.8 137.9 137.9 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-04T06:00:00 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-04T15:00:0... Attributes: short_name: TMP_2maboveground long_name: Temperature level: 2 m above ground units: K
それでは、6月1日3時UTCにおけるこのデータの分布図を表示させてみましょう。
da_gsmi.sel(time='2023-06-01T03:00').plot(cmap='RdYlGn_r',vmin=288, vmax=300)
<matplotlib.collections.QuadMesh at 0x29109c89750>
データの形式は da_msm と同じになりましたが、値についてはだいぶ異なることが分かります。
5 縮約¶
縮約とは、定められた格子領域内の格子点値を統計してその領域内の代表値を求めることです。代表値としては、平均値や最大値、最小値、標準偏差などが一般的です。GPVの処理においては、時間に関しては、時別のGPVから日別や月別のGPVを作る処理が代表的です。空間(緯度経度)に関しては、格子点間隔が小さい GPV を一定の時空間領域で集計して、より大きな格子領域を単位とする分布データを作成する処理が代表的です。 縮約は大変重要な処理ですが、一度縮約をしてしまうと、それはもはや GPV(格子点値) ではなく、メッシュデータ となることをよく認識しておいてください。
5.1 時間の縮約¶
時間方向の縮約は、データを所定の期間毎にグルーピングする処理と、夫々のグループに対する集計を組み合わせて行います。DataArrayには、時間のグルーピングのためのメソッド resample が用意され、集計には、メソッド mean や、max、min 等が用意されているので、これらを連結すれば時間方向の縮約が行えます。
それでは、時別のGPVから日平均値のGPVを作成する処理を例に、メソッド resample とメソッド mean の使い方を学びましょう。集計が正しくとれているかを容易に確認できるよう、数値予報モデルのGPVではなく、下の図のような作為的な1次元GPVを素材に使用します。
ところで、気象データを日集計するとき、気象庁は、日の範囲を、日界(00JST)が過ぎて最初のデータから次の日界を過ぎない最後のデータまでと決めています。従って、毎正時のデータを日集計する場合は、01JSTのデータから翌日00JSTのデータを集計します。
上の図の時別データは、この定義で日平均したときに、6月2日が1、6月3日が2、6月4日が3となるようになっています。
では、まずは、このようなデータの DataArray オブジェクト da を作り出しましょう。以下を実行してください
time = pd.date_range('2023-06-01T00', '2023-06-04T15:00', freq="H" )
data = [0]*16
for i in range(3):
data = data+[(i+1)]*24
da = xr.DataArray(data, coords={'time':time},dims=['time'])
tzid = 'Asia/Tokyo'
awaretime = list(pd.Series(time).dt.tz_localize('UTC').dt.tz_convert(tzid))
da = da.assign_coords(awaretime=('time', awaretime))
da.awaretime.attrs = {'long_name':'time', 'units':tzid}
da.plot(marker='o',figsize=(10,3))
[<matplotlib.lines.Line2D at 0x291099ab790>]
このデータを日平均します。 メソッド resample は多数のキーワード引数を持ちますが、気象庁の方法で UTC のデータを JST の日界でグルーピングするときには、下のように使用します。
- time 分割の単位(文字列)
日集計なので「1日」を示す文字列 '1D' を与えます。 - offset 集計を開始する時刻をずらす量
JSTの日界はUTCの日界より9時間早いので、文字列 '-9H'を与えます。 - closed 分割単位のどちらを「閉」とするか('right'または'left')
「閉」とは、境界を範囲に含める時に使用する数学の言葉です。今回は、区間の開始の時刻は含めず終了の時刻を含めるので、「閉」は右側(数直線の右側)となります。よって文字列'right' を指定します。 - label グループ化した区間に付与する時刻('right'または'left')
区間平均するので、特定の時刻の値ではありませんが、識別のために時刻のラベルが必要です。区間の最初か最後の時刻のみが指定できます。最初を指定する 'left' を与えましょう。
なお、今回用意したDataArrayオブジェクトには、座標 awaretime を付加しているので、これに基づいてグルーピングすれば集計がきわめて単純になりそうですが、残念ながらその作戦は使えません。メソッド resample の仕様により、グルーピングに際し aware な時刻は UTC に変換されてしまうので、UTC である座標 time を使用したのと同じ結果にしかなりません。
グルーピングしたものに対する集計は、集計のためのメソッド mean を連結すればそれで大丈夫です。
以下に、時別データを日本時間を日界として日平均するスクリプトを示します。実行してください。
dad = da.resample(time='1D', offset='-9H', closed='right', label='left').mean()
概要リストを表示させ、処理結果を見てみましょう。次を実行してください。
dad
<xarray.DataArray (time: 4)> array([0., 1., 2., 3.]) Coordinates: * time (time) datetime64[ns] 2023-05-31T15:00:00 ... 2023-06-03T15:00:00
日平均値が正しく計算できていることが確認できました。しかし、その一方で、縮約により、aware な時刻の座標が消失してしまいました。メソッド resample は、縮約に使用した時刻座標は再構成しますが、時間の次元に結びついた他の座標の構築まではしてくれないので、残念ですが自分で付け直さなければなりません。
DataArrayオブジェクトに新しい座標を追加するにはメソッド assign_coords を使用します。このメソッドは、自身に座標を追加するのではなく、座標が追加された新しいDataArrayオブジェクトを作るメソッドなので、既存のオブジェクトを変更するには、代入しなおす必要があります。また、UTCの時刻座標から指定したタイムゾーンの時刻の座標値を作るには、wxbcgribXの関数 awtime を用います。
tzid = 'Asia/Tokyo'
at = wx.awtime(dad,tzid)
dad = dad.assign_coords(awaretime=('time', at))
dad.awaretime.attrs = {'long_name':'time', 'units':tzid}
dad
<xarray.DataArray (time: 4)> array([0., 1., 2., 3.]) Coordinates: * time (time) datetime64[ns] 2023-05-31T15:00:00 ... 2023-06-03T15:00:00 awaretime (time) object 2023-06-01T00:00:00+09:00 ... 2023-06-04T00:00:0...
一般に、日界をシフトしてグルーピングする場合、最初と最後のグループにはデータの欠落が生じるため有効ではありません。よって、これらを除去する必要があります。それには以下のようにします。
dad = dad[1:-1] # 配列の最初と最後を除く部分で上書き
dad
<xarray.DataArray (time: 2)> array([1., 2.]) Coordinates: * time (time) datetime64[ns] 2023-06-01T15:00:00 2023-06-02T15:00:00 awaretime (time) object 2023-06-02T00:00:00+09:00 2023-06-03T00:00:00+09:00
以上見てきた通り、時別や3時間毎のGPVから日平均値のGPVを作成するには、気象庁の日界でのグルーピング、平均計算、日本標準時座標の付与、無効部分のトリミングが必要です。
参考に、4.1.1 章で用意した GSM-GPV の DataArrayオブジェクト da_gsm に対してこれらの処理を行い、日平均値のGPV da_gsm_d を計算するスクリプトを下に示しておきます。各処理がどこで実行されているか、もう皆さんにはお分かりですね。
da_gsm_d = da_gsm.resample(time='1D',
offset='-9H',
closed='right',
label='left'
).mean()[1:-1]
tzid = 'Asia/Tokyo'
da_gsm_d = da_gsm_d.assign_coords(awaretime=('time', wx.awtime(da_gsm_d,tzid)))
da_gsm_d.awaretime.attrs = {'long_name':'time', 'units':tzid}
da_gsm_d
<xarray.DataArray 'TMP_2maboveground' (time: 10, latitude: 301, longitude: 241)> array([[[301.8403 , 301.87222, 301.9028 , ..., 301.04572, 301.03854, 301.0301 ], [301.8185 , 301.84586, 301.87515, ..., 301.034 , 301.02652, 301.0187 ], [301.79117, 301.81525, 301.84323, ..., 301.02292, 301.0151 , 301.00763], ..., [287.7355 , 287.5438 , 287.28238, ..., 276.75082, 276.77426, 276.79868], [287.72964, 287.49106, 287.16714, ..., 276.64404, 276.65543, 276.67105], [287.83316, 287.56332, 287.21793, ..., 276.57828, 276.58057, 276.59164]], [[301.87967, 301.90536, 301.93436, ..., 300.92166, 300.90994, 300.8956 ], [301.8702 , 301.89105, 301.9119 , ..., 300.90897, 300.8956 , 300.88065], [301.86533, 301.8829 , 301.8979 , ..., 300.8979 , 300.8826 , 300.86826], ... [286.5559 , 286.46216, 286.33618, ..., 279.5979 , 279.6223 , 279.64185], [286.50806, 286.4231 , 286.29712, ..., 279.49048, 279.51587, 279.54224], [286.66724, 286.53345, 286.36255, ..., 279.4026 , 279.42993, 279.45728]], [[301.5842 , 301.60373, 301.6086 , ..., 300.8635 , 300.87814, 300.89377], [301.56662, 301.58224, 301.58517, ..., 300.86838, 300.88498, 300.9055 ], [301.5549 , 301.5637 , 301.5627 , ..., 300.87912, 300.89767, 300.91916], ..., [286.10568, 285.88107, 285.59103, ..., 278.16818, 278.17502, 278.17795], [286.09494, 285.88986, 285.61154, ..., 278.11252, 278.12814, 278.13596], [286.29025, 286.05295, 285.73166, ..., 278.07736, 278.093 , 278.09982]]], dtype=float32) Coordinates: * latitude (latitude) float64 20.0 20.1 20.2 20.3 ... 49.7 49.8 49.9 50.0 * longitude (longitude) float64 120.0 120.1 120.2 120.4 ... 149.8 149.9 150.0 * time (time) datetime64[ns] 2023-06-01T15:00:00 ... 2023-06-10T15:00:00 awaretime (time) object 2023-06-02T00:00:00+09:00 ... 2023-06-11T00:00:0... Attributes: short_name: TMP_2maboveground long_name: Temperature level: 2 m above ground units: K
5.2 空間の縮約について¶
メソッド resample は時間に特化したグルーピングのメソッドです。空間に対しては、メソッド groupe_by を使用しますが、残念なことに、このメソッドは、現時点では、特定の1つの次元でしかデータをグルーピングすることができません。二次元データを一筆書きの要領でむりやり一次元にして行う方法もあるようですが、とても難解です。このため、計算時間に制約がない限りは、現時点では、メソッド sel で縮約しようとする領域(Grid Cell)内のGPVを取り出して集計するという計算をループで繰り返す方が分かりやすく確実と思われます。
6 結合¶
異なる二つのGPVの結合のしかたは何種も考えられますが、ここでは、気象庁GPVの分析において使用頻度が高い3種類の結合について説明します。
6.1 特定の次元に沿った結合¶
MSM-GPV は、78時間先までの一つながりの数値予報モデルGPVですが、予報時間で分割された複数のGRIBファイルで提供されるので、それぞれのファイルから個別にGPVを読み出した場合はそれらを繋ぎ直す必要があります。
このケースでは、それぞれのGPVは、次元(時間、緯度、経度)自体は共通していて、時間の座標範囲だけが異なっています。 このようなGPVの結合には、Xarrayの関数 concat を用い、以下のようにします。
gpv_new = xr.concat( [gpv_1, gpv_2], 'time' )
ここで、gpv_1、gpv_2 は結合しようとするGPV、gpv_new は結合されたGPVのDataArrayオブジェクトです。そして、文字列「time」は、時間の次元に付けられている名前です。
ライブラリ wxbcgribx の関数 getgpv は、分割されたGRIBファイルを一括して読み込む機能を持っていますが、敢えてひとつづつ読み込んで手動で連結し、この方法を習得しましょう。まず、初期時刻~15時間先までのGPV、gpv00_15 を作ります。
# 初期時刻~15時間先までのGPV
grb_dir = f'./jmadata/msm/2023/202306'
grb_path = f'{grb_dir}/Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin'
ds = wx.getgpv(grb_path, 'TMP', timezone='Asia/Tokyo')
gpv00_15 = ds['TMP_1D5maboveground']
北緯35度、東経137度の格子点値の時系列グラフは以下の通りです。
gpv00_15.sel(latitude=35.0, longitude=137.0).plot(figsize=(5,3))
[<matplotlib.lines.Line2D at 0x29103593130>]
同様にして、16時間先~33時間先までのGPV、gpv16_33 を作ります。
# 16時間先~33時間先までのGPV
grb_path = f'{grb_dir}/Z__C_RJTD_20230601000000_MSM_GPV_Rjp_Lsurf_FH16-33_grib2.bin'
ds = wx.getgpv(grb_path, 'TMP', timezone='Asia/Tokyo')
gpv16_33 = ds['TMP_1D5maboveground']
gpv16_33.sel(latitude=35.0, longitude=137.0).plot(figsize=(5,3))
[<matplotlib.lines.Line2D at 0x29109af7220>]
用意ができたところで、両者を結合します。
gpv00_33 = xr.concat( [gpv00_15, gpv16_33], 'time' )
では、北緯35度、東経137度の格子点値の時系列グラフを描いて結合を確認しましょう。以下を実行してください。
gpv00_33.sel(latitude=35.0, longitude=137.0).plot(figsize=(10,3))
[<matplotlib.lines.Line2D at 0x2910a2f1d80>]
時間方向に連結され、連続したGPVとなっていることが確認できました。
6.2 上書きの結合¶
気象庁のGPVの多くは、新しい初期値に基づいて数時間に1回更新されています。これを利用するにあたり、「新しい予報が出されたらそれに更新するが、それより前の気象値については、以前に出された予報値をそのまま使いたい」というようなケースが時としてあります。このようなケースでは、それぞれのGPVは、次元(時間、緯度、経度)自体は共通していて、時間の座標範囲がオーバーラップしつつ異なっています。
また、広い範囲のGPVがタイルのように緯度経度領域で分割されて提供されていて、必要とする緯度経度領域を得るために取得した小分けのGPVを空間的に繋げる必要が生じることもあります。このケースでは、次元(時間、緯度、経度)自体は共通していて、緯度範囲が異なっています。
これら両方のケースにおいて、DataArrayオブジェクトのメソッド combine_first でGPVを結合することができます。combine_first はメソッドなので以下のように使用します。
gpv_new = gpv_2.combine_first(gpv_1)
ここで、gpv_1 は結合しようとするGPV、gpv_2 はこれに上書きをするGPVのDataArrayオブジェクトです。結合された gpv_new は、gpv_2 と、これがカバーしていない時空間部分のgpv_1が結合されたGPVです。gpv_new の座標は、gpv_1 も gpv_2 も納まるように自動的に拡張されます。
4.2 章で作成した、MSM-GPVの座標系に補間したGSM-GPV、da_gsmi と、それとずれた領域のMSM-GPVをこの方法で結合してみましょう。 4.2 章で読みだした領域は、北緯34~36度、東経135~138度でしたが、それより北西に1.0度ずれた領域のMSM-GPVを改めてGRIBファイルから読み込み、da_msm2 とします。
yyyymmdd = '20230601'
prdct = 'msm'
fh = ['00-15', '16-33', '34-39', '40-51', '52-78']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_MSM_GPV_Rjp_Lsurf_FH{xx}_grib2.bin'
for xx in fh]
area = [35.0, 37.0, 136.0, 139.0] # GPVを読み込む領域
# 34.0, 36.0, 135.0, 138.0 4.2章で読んだ領域
ds = wx.getgpv(grb_paths, 'TMP', timezone='Asia/Tokyo', lalomima=area)
da_msm2 = ds['TMP_1D5maboveground']
4.2 章と同じ時刻における分布図を描画します。
da_msm2.sel(time='2023-06-01T03:00').plot(cmap='RdYlGn_r',vmin=288, vmax=300)
<matplotlib.collections.QuadMesh at 0x2910a383e20>
da_gsmi と da_msm2 を、メソッド combine_first で結合します。
da_new = da_msm2.combine_first(da_gsmi)
それでは、4.2 章と同じ時刻における分布図を描いてみましょう。どちらのデータが優先されているか確認してください。
da_new.sel(time='2023-06-01T03:00').plot(cmap='RdYlGn_r',vmin=288, vmax=300)
<matplotlib.collections.QuadMesh at 0x2910b632dd0>
この例は、空間に関する結合ですが、時間について適用すれば、「新しい予報が出されたらそれに更新するが、それより前の気象値については、以前に出された予報値をそのまま使う」という結合を実現できます。
6.3 次元を追加する結合¶
これまで、地上気象予測のGPVだけを対象にしてきましたが、指定等圧面のGPVに目を向けてみましょう。MSM-GPVの上空のGPVをサンプルのGRIBファイルから読み込んで、概要フォームを表示します。
yyyymmdd = '20230601'
prdct = 'msm'
fh = ['00-15', '18-33', '36-39', '42-51', '54-78']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_MSM_GPV_Rjp_L-pall_FH{xx}_grib2.bin'
for xx in fh]
ds = wx.getgpv(grb_paths, ['HGT', 'RH', 'TMP', 'UGRD', 'VGRD', 'VVEL'], timezone='Asia/Tokyo')
ds
<xarray.Dataset> Dimensions: (latitude: 253, longitude: 241, time: 27) Coordinates: * latitude (latitude) float64 22.4 22.5 22.6 22.7 ... 47.3 47.4 47.5 47.6 * longitude (longitude) float64 120.0 120.1 120.2 ... 149.8 149.9 150.0 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-04T06:00:00 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-04T15:00... Data variables: (12/92) HGT_1000mb (time, latitude, longitude) float32 1.054 -1.946 ... 38.41 UGRD_1000mb (time, latitude, longitude) float32 4.439 4.549 ... 4.717 5.186 VGRD_1000mb (time, latitude, longitude) float32 -13.09 -13.51 ... -2.124 TMP_1000mb (time, latitude, longitude) float32 301.0 301.0 ... 277.6 277.6 VVEL_1000mb (time, latitude, longitude) float32 0.0 0.0 ... 0.01822 0.02994 RH_1000mb (time, latitude, longitude) float32 81.22 81.04 ... 90.11 90.52 ... ... VVEL_150mb (time, latitude, longitude) float32 0.0 0.0 ... 0.000834 HGT_100mb (time, latitude, longitude) float32 1.673e+04 ... 1.625e+04 UGRD_100mb (time, latitude, longitude) float32 -8.205 -8.267 ... 13.3 VGRD_100mb (time, latitude, longitude) float32 -3.854 -3.744 ... 6.574 TMP_100mb (time, latitude, longitude) float32 197.5 197.5 ... 223.9 223.8 VVEL_100mb (time, latitude, longitude) float32 0.0 0.0 ... -0.02619 Attributes: Conventions: COARDS History: created by wgrib2 GRIB2_grid_template: 0
Data variables のセクションが畳まれているのでクリックして展開してください。きわめて多くの GPV が保持されていることが分かります。そして、よく見ると、GPVが気圧毎に整理されていることがその理由であることが分かります。
特定気圧面の分析であればこれでも不自由ありませんが、例えば、気温の鉛直分布を分析する際などには、時刻、気圧、緯度、経度の4次元 GPV としてまとめたほうが効率的です。このような結合を「次元を追加する結合」と表現しました。このような結合は、次に示す3つのステップを踏んで実行します。
- 次元の追加
- 次元への座標値の付与
- 追加した次元に沿った結合
それでは、まず手始めに、1000 hP の気温のGPVに、気圧(p)の次元を追加しましょう。現在の次元:時刻(time)、緯度(latitude)、経度(longitude) の2番目に気圧(p)を追加挿入します。これにはメソッド expand_dims を用いて以下のようにします。
da = ds['TMP_1000mb'].expand_dims('p',axis=1)
da
<xarray.DataArray 'TMP_1000mb' (time: 27, p: 1, latitude: 253, longitude: 241)> array([[[[300.991 , 301.0457 , 301.01443, ..., 299.5613 , 299.5535 , 299.53787], [300.88943, 300.94412, 300.866 , ..., 299.5535 , 299.5457 , 299.52225], [300.78787, 300.84256, 300.81912, ..., 299.53787, 299.52225, 299.50662], ..., [288.50662, 288.491 , 287.991 , ..., 278.65506, 278.63162, 278.59256], [288.03787, 288.10037, 288.07693, ..., 278.6082 , 278.59256, 278.5457 ], [287.9207 , 287.75662, 287.6863 , ..., 278.4988 , 278.491 , 278.46756]]], [[[300.91122, 300.94247, 300.98154, ..., 299.80185, 299.79404, 299.7706 ], [300.80966, 300.8409 , 300.9503 , ..., 299.78622, 299.7706 , 299.7628 ], [300.72372, 300.73154, 301.49716, ..., 299.7784 , 299.7628 , ... 277.52908], [294.99002, 295.20096, 295.31815, ..., 277.56033, 277.5447 , 277.4822 ], [294.80252, 294.77127, 294.77908, ..., 277.52127, 277.50565, 277.45877]]], [[[300.88187, 300.88968, 300.94437, ..., 299.56937, 299.60062, 299.62405], [300.7178 , 300.74905, 300.80374, ..., 299.53812, 299.56155, 299.60062], [300.5928 , 300.60843, 300.8428 , ..., 299.49124, 299.50687, 299.53812], ..., [296.96 , 297.04593, 296.63968, ..., 277.70218, 277.68655, 277.63968], [296.41312, 296.55374, 296.61624, ..., 277.66312, 277.6475 , 277.60062], [296.31937, 296.22562, 296.21 , ..., 277.61624, 277.60843, 277.585 ]]]], dtype=float32) Coordinates: * latitude (latitude) float64 22.4 22.5 22.6 22.7 ... 47.3 47.4 47.5 47.6 * longitude (longitude) float64 120.0 120.1 120.2 120.4 ... 149.8 149.9 150.0 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-04T06:00:00 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-04T15:00:0... Dimensions without coordinates: p Attributes: short_name: TMP_1000mb long_name: Temperature level: 1000 mb units: K
概要フォームの最上段を見てください。意図した通り、次元 p が2番目追加挿入されました。しかし、Coordinates のセクションを見ると、ここに p は存在しません。まだ座標が与えられていないためです。
それでは次に、座標を与えましょう。これには、メソッド assign_coords を用いて以下のようにします。ここでは、名前と単位も設定しています。
da = da.assign_coords({'p': [1000]})
da.name = 'TMP'
da['p'].attrs = {'long_name':'atmospheric pressure', 'units':'hPa'}
da.attrs['short_name'] = da.name
da.attrs['level'] = ''
da
<xarray.DataArray 'TMP' (time: 27, p: 1, latitude: 253, longitude: 241)> array([[[[300.991 , 301.0457 , 301.01443, ..., 299.5613 , 299.5535 , 299.53787], [300.88943, 300.94412, 300.866 , ..., 299.5535 , 299.5457 , 299.52225], [300.78787, 300.84256, 300.81912, ..., 299.53787, 299.52225, 299.50662], ..., [288.50662, 288.491 , 287.991 , ..., 278.65506, 278.63162, 278.59256], [288.03787, 288.10037, 288.07693, ..., 278.6082 , 278.59256, 278.5457 ], [287.9207 , 287.75662, 287.6863 , ..., 278.4988 , 278.491 , 278.46756]]], [[[300.91122, 300.94247, 300.98154, ..., 299.80185, 299.79404, 299.7706 ], [300.80966, 300.8409 , 300.9503 , ..., 299.78622, 299.7706 , 299.7628 ], [300.72372, 300.73154, 301.49716, ..., 299.7784 , 299.7628 , ... 277.52908], [294.99002, 295.20096, 295.31815, ..., 277.56033, 277.5447 , 277.4822 ], [294.80252, 294.77127, 294.77908, ..., 277.52127, 277.50565, 277.45877]]], [[[300.88187, 300.88968, 300.94437, ..., 299.56937, 299.60062, 299.62405], [300.7178 , 300.74905, 300.80374, ..., 299.53812, 299.56155, 299.60062], [300.5928 , 300.60843, 300.8428 , ..., 299.49124, 299.50687, 299.53812], ..., [296.96 , 297.04593, 296.63968, ..., 277.70218, 277.68655, 277.63968], [296.41312, 296.55374, 296.61624, ..., 277.66312, 277.6475 , 277.60062], [296.31937, 296.22562, 296.21 , ..., 277.61624, 277.60843, 277.585 ]]]], dtype=float32) Coordinates: * latitude (latitude) float64 22.4 22.5 22.6 22.7 ... 47.3 47.4 47.5 47.6 * longitude (longitude) float64 120.0 120.1 120.2 120.4 ... 149.8 149.9 150.0 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-04T06:00:00 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-04T15:00:0... * p (p) int32 1000 Attributes: short_name: TMP long_name: Temperature level: units: K
Coordinates のセクションを見てください。次元 p に座標値 1000 が付与され、名称、単位も設定されました。
ここまでくれば、あとは、6.1 章で学習した、「特定の次元に沿った結合」をすれば完成です。今回は全部で15個のGPVを結合するので、ループを組むことにします。
# ループでの実行
levels = [1000,975,950,900,850,800,700,600,500,400,300,250,200,150,100]
for level in levels[1:]: #1000hPについては既に作ってあるので除く
dc = ds[f'TMP_{level}mb'].expand_dims('p',axis=1)
dc = dc.assign_coords({'p': [level]})
da = xr.concat([da,dc], 'p')
さあ、うまく結合できたでしょうか。概要フォームを表示して結果を確認してみましょう。
da
<xarray.DataArray 'TMP' (time: 27, p: 15, latitude: 253, longitude: 241)> array([[[[300.991 , 301.0457 , 301.01443, ..., 299.5613 , 299.5535 , 299.53787], [300.88943, 300.94412, 300.866 , ..., 299.5535 , 299.5457 , 299.52225], [300.78787, 300.84256, 300.81912, ..., 299.53787, 299.52225, 299.50662], ..., [288.50662, 288.491 , 287.991 , ..., 278.65506, 278.63162, 278.59256], [288.03787, 288.10037, 288.07693, ..., 278.6082 , 278.59256, 278.5457 ], [287.9207 , 287.75662, 287.6863 , ..., 278.4988 , 278.491 , 278.46756]], [[298.80154, 298.86404, 298.86404, ..., 297.4656 , 297.44217, 297.43436], [298.69998, 298.7703 , 298.72342, ..., 297.44998, 297.43436, 297.41873], [298.59842, 298.67654, 298.69998, ..., 297.43436, 297.41092, 297.4031 ], ... [221.07849, 221.0863 , 221.09412, ..., 225.03943, 225.07068, 225.09412], [221.04724, 221.06287, 221.07849, ..., 225.09412, 225.09412, 225.10193], [221.0238 , 221.03943, 221.06287, ..., 225.05505, 225.0238 , 225.01599]], [[194.13713, 194.16057, 194.19963, ..., 200.00432, 199.97307, 199.96526], [194.16838, 194.1762 , 194.21526, ..., 200.09026, 200.02776, 200.01213], [194.23088, 194.22307, 194.20744, ..., 200.19182, 200.09807, 200.0512 ], ..., [218.78557, 218.81682, 218.84807, ..., 223.44963, 223.44182, 223.434 ], [218.81682, 218.85588, 218.90276, ..., 223.65276, 223.684 , 223.69963], [218.8715 , 218.90276, 218.95744, ..., 223.84026, 223.85588, 223.84026]]]], dtype=float32) Coordinates: * latitude (latitude) float64 22.4 22.5 22.6 22.7 ... 47.3 47.4 47.5 47.6 * longitude (longitude) float64 120.0 120.1 120.2 120.4 ... 149.8 149.9 150.0 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-04T06:00:00 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-04T15:00:0... * p (p) int32 1000 975 950 900 850 800 ... 400 300 250 200 150 100 Attributes: short_name: TMP long_name: Temperature level: units: K
うまく結合できたようです。それでは、北緯35度、東経137度の格子点における気温の鉛直分布を可視化してみましょう。
縦軸を特殊な形で付けるので、メソッド plot ではなく、ライブラリ matplotlib を用いて個々のアイテムを一つ一つ設定します。
from matplotlib import pyplot as plt
from matplotlib import ticker as mt
gpv = da.sel(time=slice('2023-06-01T15:00','2023-06-02T15:00'), latitude=35.0, longitude=137.0)
lbl = [f'{str(xx)[:13]}Z' for xx in gpv.time.data]
ytic = [1000,900,800,700,600,500,400,300,200,100]
fig, ax = plt.subplots(figsize=(5,5),ncols=1,nrows=1)
ax.plot(gpv.data.T, gpv.p.data, label=lbl)
ax.set_yscale('log')
ax.set_yticks(ytic)
ax.get_yaxis().set_major_formatter(mt.ScalarFormatter())
ax.invert_yaxis()
ax.set_ylabel(f'{gpv.p.long_name} ({gpv.p.units})')
ax.set_xlabel(f'{gpv.long_name} [{gpv.attrs["units"]}]')
ax.set_title(f'(N{gpv.latitude.data:0.2f}, E{gpv.longitude.data:0.2f})')
ax.legend()
<matplotlib.legend.Legend at 0x2910b8717e0>
7 アンサンブル数値予報モデル GPV プロダクト¶
アンサンブル(ensemble)とは「集団」という意味を持ちます。アンサンブル数値予報とは、ある一つの気象予測モデルを、少しづつ微妙に変えた多数の初期気象状態で実行した結果(メンバーと呼びます)の集合のことです。
気象予測は、数値予報モデルと予測をスタートする時点での大気の状態(初期値)の両方があって初めて行うことができますが、ある時点における全地球上の大気状態を知ることは不可能なので、モデルに与える初期値には、様々な観測値で得られた断片的な情報から推定したものを用いるしかありません。このため、帰結としての予測値もまた実際とは必ず乖離します。しかも、気象のメカニズムは大変複雑なので、その乖離は千差万別です。
とはいえ、初期値が少しぐらい違っても予測結果が概ね同じようなものとなる場合、その予測は当たりそうだと期待できます。逆に、初期条件がほんの少し違うだけで予測結果がばらばらになる場合は、予測の信頼性は低いということが推察できます。このように、アンサンブルメンバーを多数用意すれば、それらを統計することで予報の信頼性や予報値の範囲などを推定することが可能になります。アンサンブル数値予報モデルGPVはこのような理由から作成されています。
7.1 気象庁のアンサンブルプロダクト¶
現在、気象庁はそれぞれに特徴を持った3つのアンサンブル数値予報システムを運用しており、それらに基づき以下9種類のアンサンブル数値予報モデルGPVを作成し、気象業務支援センターから提供しています。
- メソアンサンブル予報システム
資料:配信資料に関する仕様No.13101( https://www.data.jma.go.jp/suishin/shiyou/pdf/no13101 )- メソアンサンブル数値予報モデルGPV
- 全球アンサンブル予報システム
資料:配信資料に関する仕様No.12802( https://www.data.jma.go.jp/suishin/shiyou/pdf/no12802 )- 週間アンサンブル数値予報モデルGPV(日本域)
- 週間アンサンブル数値予報モデルGPV(全球域)
- 2週間アンサンブル数値予報モデルGPV(日本域)
- 2週間アンサンブル数値予報モデルGPV(全球域)
- 1か月アンサンブル数値予報モデルGPV(日本域)
- 1か月アンサンブル数値予報モデルGPV(全球域)
- 台風アンサンブル数値予報モデルGPV(日本域)
- 季節アンサンブル予報システム
資料:配信資料に関する仕様No.20114( https://www.data.jma.go.jp/suishin/shiyou/pdf/no20114 )- 6か月アンサンブル数値予報モデルGPV(全球域)
下図に、上記プロダクトのうち、日本周辺域 を対象とし 定常的 に配信される5種類のプロダクトについて、地上気象要素の予報期間とグリッド間隔の関係を整理したものを示します。
M1:メソアンサンブル数値予報モデルGPV、G1:週間アンサンブル数値予報モデルGPV(日本域)、G2:2週間アンサンブル数値予報モデルGPV(日本域)、G3:1か月アンサンブル数値予報モデルGPV(日本域)、S1:6か月アンサンブル数値予報モデルGPV(全球域)
7.2 アンサンブル平均/スプレッド¶
第 7 章を学習するにあたり、もう一度、カーネルをリスタートしましょう。ツールバーの輪矢印アイコンをクリックし、確認に対して [Restart] をクリックしてください。
# ライブラリーのインポート
import numpy as np
import xarray as xr
import pandas as pd
import wxbcgribx as wx
アンサンブル数値予報モデルGPV に対する分析で最も基本的なものはメンバーの統計です。中でもよく計算されるのは平均と標準偏差で、それぞれ、アンサンブル平均、アンサンブルスプレッドと呼ばれます。
それでは、メソアンサンブル予報システムのプロダクト、メソアンサンブル数値予報モデルGPV のデータを素材に、これらの計算方法を学びます。
使用するのは、初期値が2023年6月1日00UTC のもので、GRIBファイルは jmadata/meps/2023/202306/ に、保存されています。予報期間に基づき3分割になっています。
Z__C_RJTD_20230601000000_MEPS_GPV_Rjp_Lsurf_FH00-15_grib2.bin
Z__C_RJTD_20230601000000_MEPS_GPV_Rjp_Lsurf_FH18-33_grib2.bin
Z__C_RJTD_20230601000000_MEPS_GPV_Rjp_Lsurf_FH36-39_grib2.bin
以下を実行すると、jmadata/meps に保存されるサンプルのGRIBファイルからデータを取り出し、DataSetオブジェクトにまとめて概要フォームを表示します。データの次元に、他のGPVには見られない次元 member が存在することを確認してください。
# GRIBファイルへのパスのリストを作成
yyyymmdd = '20230601'
prdct = 'meps'
fd = ['00-15', '18-33', '36-39']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_MEPS_GPV_Rjp_Lsurf_FH{xx}_grib2.bin'
for xx in fd]
# DataSetオブジェクトとして読み込み
ds = wx.getgpv(grb_paths, ['APCP', 'DSWRF', 'PRMSL', 'TMP', 'UGRD', 'VGRD'],timezone='Asia/Tokyo')
# 概要フォームの表示
ds
<xarray.Dataset> Dimensions: (latitude: 505, longitude: 481, time: 14, member: 21) 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] 2023-06-01 ... 2023-06-02T15:0... * member (member) <U14 'ENS=hi-res ctl' 'ENS=+1' ... 'ENS=-10' awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-... Data variables: PRMSL_meansealevel (member, time, latitude, longitude) float32 1e+05 ..... UGRD_10maboveground (member, time, latitude, longitude) float32 4.583 ..... VGRD_10maboveground (member, time, latitude, longitude) float32 -13.51 .... TMP_1D5maboveground (member, time, latitude, longitude) float32 301.2 ..... APCP_surface (member, time, latitude, longitude) float32 nan ... ... DSWRF_surface (member, time, latitude, longitude) float32 nan ... 0.0 Attributes: Conventions: COARDS History: created by wgrib2 GRIB2_grid_template: 0
ではさらに、名古屋地方気象台に最も近い格子点における地上気温をこれから取り出します。
lat, lon = 35.1667, 136.965 # 名古屋地方気象台の緯度経度
gpv = ds['TMP_1D5maboveground'].sel(latitude=lat, longitude=lon, method='nearest')
gpv
<xarray.DataArray 'TMP_1D5maboveground' (member: 21, time: 14)> array([[295.0896 , 296.79352, 296.23358, 294.74805, 292.9697 , 292.90015, 293.19348, 293.74637, 296.52814, 296.47235, 296.58392, 293.203 , 292.4676 , 291.8237 ], [295.11102, 296.65582, 296.06918, 294.67532, 292.63626, 292.60825, 292.87036, 293.1708 , 296.00955, 296.3969 , 296.81256, 294.6732 , 293.2252 , 292.5018 ], [295.06064, 296.07883, 295.51517, 294.4176 , 292.99554, 292.71176, 293.35788, 294.13568, 296.30704, 296.8293 , 296.98157, 293.84094, 292.63123, 292.1224 ], [294.8644 , 296.91516, 295.457 , 294.43533, 292.86234, 292.2871 , 293.05737, 293.75003, 296.4486 , 296.3727 , 296.61417, 296.66064, 294.90256, 293.77625], [295.30704, 297.52606, 297.27225, 295.39423, 293.56497, 293.06198, 293.16452, 293.4872 , 296.1985 , 296.34918, 295.72696, 294.004 , 292.27393, 291.48074], [295.0015 , 297.149 , 296.63647, 293.9147 , 293.25888, 293.31924, 293.91998, 294.48935, 296.10822, 296.58905, 296.55432, 296.1316 , 296.54456, 293.80905], [295.16418, 296.89685, 296.1231 , 292.93777, 292.18 , 292.0579 , 291.56036, 291.63596, 293.9129 , 296.6176 , 296.25635, 294.32214, ... 291.30646, 291.44 , 292.71875, 293.06876, 293.30972, 292.89023, 292.41397, 291.84595], [295.2439 , 297.76868, 296.45032, 293.31287, 293.01025, 293.57825, 293.479 , 294.42337, 296.7591 , 296.0409 , 292.68652, 292.25858, 291.4712 , 291.08987], [294.92755, 296.87872, 296.64206, 294.7372 , 292.88193, 292.484 , 292.72348, 293.4066 , 295.45584, 296.17963, 296.44165, 296.5737 , 296.4077 , 295.21387], [295.0934 , 297.09436, 296.15445, 294.7468 , 293.32007, 292.82968, 293.244 , 293.6235 , 296.17282, 296.7136 , 296.89398, 296.87567, 294.55267, 293.51465], [295.07812, 296.73468, 296.2501 , 294.8379 , 292.61447, 292.68567, 292.79202, 293.63235, 296.18344, 296.11523, 296.49234, 292.926 , 291.89554, 291.63385], [295.15457, 297.40176, 296.91602, 294.8079 , 292.60458, 292.41406, 292.93283, 294.06543, 296.68927, 296.86383, 293.7718 , 292.53836, 291.44098, 291.06808], [295.01862, 296.94186, 296.1076 , 294.52103, 292.90314, 293.54337, 293.58298, 294.0557 , 296.42444, 296.2921 , 296.12048, 296.1574 , 296.27768, 295.58887]], dtype=float32) Coordinates: latitude float64 35.15 longitude float64 136.9 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-02T15:00:00 * member (member) <U14 'ENS=hi-res ctl' 'ENS=+1' ... 'ENS=+10' 'ENS=-10' awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-03T00:00:0... Attributes: short_name: TMP_1D5maboveground long_name: Temperature level: 1.5 m above ground units: K
可視化してみましょう。21本の折れ線からなるこのこのグラフの場合、DataArrayオブジェクトのメソッド plot ではどうしても見にくいものとなるので、ライブラリ Matplotlib の 機能を使ってグラフを作ります。
import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(12,4),ncols=1,nrows=1)
ax.plot(gpv.awaretime.data, gpv.data.T, label=gpv.member.data) # 配列の形を整合させるために転置(.T)が必要
ax.set_ylabel(f'{gpv.long_name} [{gpv.attrs["units"]}]')
ax.set_xlabel(f'{gpv.awaretime.long_name} ({gpv.awaretime.units})')
ax.set_title(f'{gpv.name} (N{gpv.latitude.data:0.2f}, E{gpv.longitude.data:0.2f})')
ax.legend(bbox_to_anchor=(0.5, -0.15), loc='upper center', borderaxespad=0, fontsize=8, ncol=9)
<matplotlib.legend.Legend at 0x2910a3122c0>
それではアンサンブル平均を取ります。それには、DataArrayオブジェクトのメソッド mean を使用します。引数には、平均を取る次元を指定します。今回はアンサンブルメンバーに対して平均を取るので、'member' を指定します。
gpv.mean('member')
<xarray.DataArray 'TMP_1D5maboveground' (time: 14)> array([295.08505, 296.97516, 296.243 , 294.5146 , 292.89426, 292.77798, 292.9975 , 293.58765, 295.74927, 296.2895 , 296.04816, 294.64175, 293.43106, 292.44354], dtype=float32) Coordinates: latitude float64 35.15 longitude float64 136.9 * time (time) datetime64[ns] 2023-06-01 ... 2023-06-02T15:00:00 awaretime (time) object 2023-06-01T09:00:00+09:00 ... 2023-06-03T00:00:0...
21メンバーの気温が平均され、1次元の配列となりました。
可視化してみましょう。気合を入れて、メンバーとアンサンブル平均の両方を一つのグラフにします。
fig, ax = plt.subplots(figsize=(12,4),ncols=1,nrows=1)
ax.plot(gpv.awaretime.data, gpv.data.T, label=gpv.member.data) # 配列の形を整合させるために転置(.T)が必要
ax.plot(gpv.awaretime.data, gpv.mean('member'), # アンサンブル平均を同じ軸空間に重ね書き
label='ensamble mean',
color='black', linewidth=4, marker="o"
)
ax.set_ylabel(f'{gpv.long_name} [{gpv.attrs["units"]}]')
ax.set_xlabel(f'{gpv.awaretime.long_name} ({gpv.awaretime.units})')
ax.set_title(f'{gpv.name} (N{gpv.latitude.data:0.2f}, E{gpv.longitude.data:0.2f})')
ax.legend(bbox_to_anchor=(0.5, -0.15), loc='upper center', borderaxespad=0, fontsize=8, ncol=9)
<matplotlib.legend.Legend at 0x2910ced3b20>
アンサンブルスプレッドもを可視化しましょう。一般に標準偏差や総和は、メンバーの値からは大きく異なるので、描画のテクニックとして、Y軸を別に設けなければなりません。
fig, ax = plt.subplots(figsize=(12,4),ncols=1,nrows=1)
ax.plot(gpv.awaretime.data, gpv.data.T, label=gpv.member.data)
ax.set_ylabel(f'{gpv.long_name} [{gpv.attrs["units"]}]')
ax.set_xlabel(f'{gpv.awaretime.long_name} ({gpv.awaretime.units})')
ax.set_title(f'{gpv.name} (N{gpv.latitude.data:0.2f}, E{gpv.longitude.data:0.2f})')
ax.legend(bbox_to_anchor=(0.5, -0.15), loc='upper center', borderaxespad=0, fontsize=8, ncol=9)
ax2 = ax.twinx() # 第2軸をとるにはこうします。
ax2.plot(gpv.awaretime.data, gpv.std(dim='member'),
label='ensamble splead',
color='black', linewidth=4, marker="o"
)
ax2.set_ylabel(f'Spread of {gpv.long_name} [{gpv.attrs["units"]}]')
ax2.legend(loc='best')
<matplotlib.legend.Legend at 0x2910d1da950>
上の処理では、特定地点の時系列データを取り出したうえでアンサンブルに対する統計を取りましたが、必ずしもこの順序である必要はありません。まず最初にアンサンブルに対する統計を取って3次元のGPVにし、その後に特定地点の時系列データを取り出しても構いません。
DataArrayオブジェクトの次元を統計するメソッドは、このほかに、max (最大値)、min (最小値)、var (分散)、median (中央値)、sum (総和)などがあります。
7.3 アンサンブルGPVの結合¶
7.3.1 メンバーの連結¶
下の図を見てください。右側は、「週間アンサンブル数値予報モデルGPV(日本域)」で、左側は「2週間アンサンブル数値予報モデルGPV(日本域)」から特定の格子点における気温予測を取り出してグラフ化したものです。これらのプロダクトは、同じ全球アンサンブルシステムの同一の初期値からの計算結果であり、元来はひとつながりのものです。言ってしまえば、両者は1つのものを2つに分けて違う名前でばら売りしているものなので、両者を接続して利用したほうが便利なことがあります。
同じ初期日の週間アンサンブル数値予報モデルGPV(日本域) (左)と 2週間アンサンブル数値予報モデルGPV(日本域) (右)の例
分割されているひとつながりのアンサンブルプロダクトを結合するには2つの方法があります。第1の方法は、連結すべきGPVのGRIBファイルを、関数 getgpv で一度に読み込む方法です。関数 getgpv は、読み込みを指定された複数のGPVの時間軸を見て、これらが重なっていなければ、メンバー名が同じもの同士を時間方向に連結するように作られています。関数に与えるリストを連結するだけなの簡単です。通常はこの方法を用いればよいでしょう。
第2の方法は、6 章で学習した、関数 concat で結合する方法です。何らかの理由で、すでに二つのオブジェクトになっている場合に使用してください。今のケースでは、緯度、経度、メンバー名が同じ2つのデータを時刻方向に連結するので、関数を以下のようにして用います。
gpv_new = xr.concat( [gpv_1, gpv_2], 'time' )
ここで、 gpv_1 と gpv_2 は結合する DataArray オブジェクト、gpv_new は結合されたオブジェクトです
7.3.2 メンバーの合成¶
アンサンブル予報の分析においてはメンバーについての統計が重要ですが、そのためには統計に耐えるメンバー数が必要です。十分なメンバー数は分析によって異なると思いますが、メソアンサンブル数値予報モデルGPV(21)や1か月アンサンブル数値予報モデルGPV(25)、6か月アンサンブル数値予報モデルGPV(5)のメンバー数は明らかに不十分です(下表参照)。これらに対しては、初期日時が異なる複数の予報を一まとめにしてメンバー数を増やしてから分析をする必要があります。
プロダクト | 予報期間 | 時間ステップ | 予報頻度 | メンバー数 | 地上面格子サイズ |
---|---|---|---|---|---|
メソアンサンブル数値予報モデルGPV | 39時間 | 3時間 | 4回/日 | 21 | 0.05度×0.0625度 |
週間アンサンブル数値予報モデルGPV(日本域) | 264時間 | 3時間 | 2回/日 | 51 | 0.375度×0.375度 |
2週間アンサンブル数値予報モデルGPV(日本域) | 267~432時間 | 3時間 | 1回/日 | 51 | 0.375度×0.375度 |
1か月アンサンブル数値予報モデルGPV(日本域) | 435~816時間 | 3時間 | 1回/7日 | 火曜日25,水曜日25 | 0.375度×0.375度 |
6か月アンサンブル数値予報モデルGPV(全球域) | 240日 | 1日 | 1回/日 | 5 | 1.25度× 1.25度 |
それでは、1か月アンサンブル数値予報モデルGPV(日本域)を素材とし、火曜日と水曜日に発表された2日分のプロダクトを結合して50メンバーのアンサンブル予報データを合成してみしましょう。
これにも二通りの方法があります。第1の方法は、7.3.1 章と同じ方法で、火曜日の予報のGRIBファイルと水曜日の予報のGRIBファイルを、一つのファイルリストにして一度に読む方法です。関数 getgpv は、読み込みを指定された複数のGPVの時間軸を見て、これらが重なっていなければ同じ名前のメンバー同士を時間方向に連結しますが、重なりがある場合は別なメンバーとして並べるように作られています。今回のケースでは、二つのアンサンブルGPVは期間が1日違うだけで大部分は重なっているので、メンバーは連結されずに別々に扱われます。
それでは、第1の方法でメンバーを合成します。まず、火曜日の予報のGRIBファイルへのパスと、水曜日の予報のそれとを一緒のリストにします。
# 第1の方法
# 火曜日の予報のGRIBファイルへのパスと、水曜日の予報のそれとを一緒のリストにする
yyyymmdd = '20230530'
yyyymmdds = wx.strft_range(start=yyyymmdd, format='%Y%m%d', iter=2, days=1)
prdct = 'geps/g3'
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}120000_EPSG_GPV_Rjp_Gll0p375deg_Lsurf_FD1803-3400_grib2.bin'
for yyyymmdd in yyyymmdds]
grb_paths
['./jmadata/geps/g3/2023/202305/Z__C_RJTD_20230530120000_EPSG_GPV_Rjp_Gll0p375deg_Lsurf_FD1803-3400_grib2.bin', './jmadata/geps/g3/2023/202305/Z__C_RJTD_20230531120000_EPSG_GPV_Rjp_Gll0p375deg_Lsurf_FD1803-3400_grib2.bin']
# 第1の方法
# GPVを一度に読み込む
ds = wx.getgpv(grb_paths, 'TMP', timezone='Asia/Tokyo')
ds
<xarray.Dataset> Dimensions: (latitude: 83, longitude: 83, time: 136, member: 50) Coordinates: * latitude (latitude) float64 19.5 19.88 20.25 ... 49.5 49.88 50.25 * longitude (longitude) float64 119.6 120.0 120.4 ... 150.0 150.4 * time (time) datetime64[ns] 2023-06-17T15:00:00 ... 2023-07-... * member (member) <U15 'ENS=low-res ctl' 'ENS=+1' ... 'ENS=-12' awaretime (time) object 2023-06-18T00:00:00+09:00 ... 2023-07-04... Data variables: TMP_2maboveground (member, time, latitude, longitude) float32 301.1 ... ... Attributes: Conventions: COARDS History: created by wgrib2 GRIB2_grid_template: 0
DataSetオブジェクトの概要フォームには、メンバー数が50と示されています。メンバーが別なものとして合成されました。
Coordinats のセクション、member の項目の右端のお皿アイコンをクリックして、メンバー名を表示させてください。それぞれの GPV で使われていたメンバー名がそのまま使われています。アンサンブル数値予報モデル GPV の個々のメンバーは、メソッド sel にメンバー名を指定することで選び出せますが、この方法で合成されたGPVに適用すると、同じ名前のメンバーが常に2つ選びだされ、区別することはできません。ただし、データ自体は独立なので、インデックス番号がわかっていれば区別することは可能です。関数 getgpv は初期時刻がより早いものにより小さいメンバーのインデックスを付けるので、識別子指定で取り出された二つのメンバーのうち、インデックスが0のものが火曜日、1のものが水曜日のものです。
また、このようにしてメンバー数を拡張したGPVでは、予報期間が重なっていない部分については、データが使えないことに注意してください。
以上を確認するため、合成した1か月アンサンブル予報GPVから特定格子点の気温予測を取り出し、さらに、メンバー名が「ENS=low-res ctl」のものを取り出して、折れ線グラフにしてみましょう。
lat, lon = 35.1667, 136.965 # 名古屋地方気象台の緯度経度
gpv = ds['TMP_2maboveground'].sel(latitude=lat, longitude=lon, method='nearest').sel(member='ENS=low-res ctl')
gpv.plot.line(x='time', figsize=(12,4))
gpv
<xarray.DataArray 'TMP_2maboveground' (member: 2, time: 136)> array([[290.04294, 289.31305, 289.825 , 294.71313, 297.9852 , 299.32175, 297.44168, 293.35687, 291.4228 , 289.61658, 290.7325 , 295.6193 , 297.8225 , 297.9764 , 296.26535, 294.15665, 293.62253, 292.50964, 292.641 , 296.69226, 298.47787, 298.41724, 295.9886 , 294.62006, 294.37402, 293.9212 , 293.8212 , 297.22025, 298.9618 , 299.0916 , 296.781 , 294.22452, 293.13254, 292.7237 , 293.13068, 295.6802 , 296.30334, 296.60757, 295.9524 , 293.9651 , 292.81906, 292.1818 , 292.55823, 296.92892, 299.938 , 300.48813, 299.3222 , 295.7342 , 293.77576, 293.5261 , 294.39777, 297.7011 , 298.1187 , 296.64758, 294.2007 , 291.83542, 291.06693, 290.65814, 290.83707, 290.98962, 292.11948, 294.47308, 294.20648, 293.15643, 292.9085 , 292.96704, 293.14792, 296.1407 , 298.4708 , 299.11584, 296.49106, 292.54395, 290.47876, 289.41318, 289.92532, 294.83374, 297.88516, 298.65732, 296.9463 , 293.45398, 292.16437, 291.37094, 291.99045, 294.09668, 295.4175 , 295.2473 , 294.08453, 292.3226 , 291.49042, 291.19858, 291.27237, 293.283 , 294.842 , 296.36322, 295.6581 , 293.93338, 293.3153 , 292.84402, 293.14362, 294.81982, 294.98584, 295.69577, 294.68747, 293.47592, 293.08136, 292.60144, 292.49124, 294.9089 , 295.44498, 295.98035, 294.50992, 293.59976, 293.35474, 293.51642, 294.00638, 295.72208, 299.29572, 297.60782, 297.068 , 295.74332, ... 296.11826, 297.4969 , 297.58786, 297.14606, 295.93274, 295.37943, 295.3383 , 295.09512, 294.74457, 295.05457, 295.01486, 295.82538, 296.01154, 294.6325 , 294.34027, 293.8954 , 293.8661 , 294.44662, 294.4503 , 294.86075, 294.69443, 294.51013, 294.30682, 293.79825, 294.01303, 295.56354, 298.41428, 299.2057 , 298.35947, 295.62335, 294.90518, 294.15427, 294.73056, 297.90637, 298.0101 , 298.17938, 297.2651 , 296.50177, 296.08356, 295.5933 , 295.90753, 298.27612, 299.3902 , 299.59854, 298.38477, 297.2344 , 296.5504 , 294.93304, 293.34045, 296.28522, 299.65012, 300.35925, 298.361 , 294.068 , 291.90366, 290.13046, 290.3598 , 296.11685, 299.14813, 297.85873, 296.65256, 294.96274, 293.51117, 292.2894 , 292.19064, 294.53427, 297.8818 , 300.3011 , 299.6792 , 296.5864 , 295.7861 , 294.8624 , 295.39352, 300.44147, 302.41098, 302.82336, 301.32556, 297.84973, 295.9389 , 295.12836, 296.0739 , 300.69806, 301.9912 , 302.06094, 300.16858, 296.89093, 294.85428, 293.69778, 293.79 , 300.0881 , 302.45297, 302.6479 , 300.12964, 296.39075, 294.68234, 294.13083, 294.52948, 300.26877, 302.95233, 303.17838, 300.35632, 297.70496, 296.17172, 294.88953, 294.9279 , 300.12463, 302.69232, 303.06354, 300.80423, 297.48502, 295.57303, 294.50937, 294.89084, 299.72034, 302.8696 , 303.5492 , 301.90976, 298.51944]], dtype=float32) Coordinates: latitude float64 35.25 longitude float64 136.9 * time (time) datetime64[ns] 2023-06-17T15:00:00 ... 2023-07-04T12:00:00 * member (member) <U15 'ENS=low-res ctl' 'ENS=low-res ctl' awaretime (time) object 2023-06-18T00:00:00+09:00 ... 2023-07-04T21:00:0... Attributes: short_name: TMP_2maboveground long_name: Temperature level: 2 m above ground units: K
メンバーを合成する第2の方法も、関数 concat で行う方法です。ちょっと分かりにくいかもしれませんが、このとき、指定する次元は member です。
gpv_new = xr.concat( [gpv_1, gpv_2], 'member' )
ここで、 gpv_1 と gpv_2 は結合する DataArray オブジェクト、gpv_new は結合されたオブジェクトです。
何らかの必要からメンバー名をユニークな識別子で間違いなく管理したい場合は、結合前にメソッド assign_coords ( https://docs.xarray.dev/en/stable/generated/xarray.DataArray.assign_coords.html )を使用し、重複しないメンバー名を与え直してから結合するとよいでしょう。
7.4 降水量には気を付けましょう¶
気象庁GPVデータ分析チャレンジ!入門 においても説明しましたが、気象庁のGPVプロダクトにおける降水量データに対しては、注意すべき点が3つあります。これらはアンサンブル数値予報モデルGPVについても言えることですので注意してください。
- 空間に関しては格子点値だが、時間に関して瞬時値ではなく時間ステップ内で集計された"メッシュ値"である
- 時間ステップ内の集計値であることから、初期時刻の値が存在しない
- 「どのような集計か?」がプロダクトによって異なる
アンサンブルGPVであることを活かす分析としては、予定される降水の降り出し時刻の分布や特定地点で予想される降水量積算値の分布が考えられます。一方、いくつかのケースでは、降水量の単純なアンサンブル平均は、全く意味をなさないことがあります。
以下に、サンプルとして収録する、メソアンサンブル数値予報モデルGPV、週間アンサンブル数値予報モデルGPV(日本域)、6か月アンサンブル数値予報モデルGPVの降水量を、極簡易な図にするスクリプトを掲載するので、グラフを確認し、各自の関心に照らしたアンサンブルならではの利用法や、しても意味がないことなどを考察してください。
# メソアンサンブル数値予報モデル
yyyymmdd = '20230601'
prdct = 'meps'
fd = ['00-15', '18-33', '36-39']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_MEPS_GPV_Rjp_Lsurf_FH{xx}_grib2.bin'
for xx in fd]
ds = wx.getgpv(grb_paths, 'APCP', timezone='Asia/Tokyo')
gpv = ds['APCP_surface'].sel(latitude=35.1667, longitude=136.965, method='nearest')
wx.ts(gpv.time.data, gpv.data )
# 週間アンサンブル数値予報モデルGPV
yyyymmdd = '20230531'
prdct = 'geps/g1'
fd = ['0000-0512', '0515-1100']
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = [f'{grb_dir}/Z__C_RJTD_{yyyymmdd}120000_EPSG_GPV_Rjp_Gll0p375deg_Lsurf_FD{xx}_grib2.bin'
for xx in fd]
ds = wx.getgpv(grb_paths, 'APCP', timezone='Asia/Tokyo')
gpv = ds['APCP_surface'].sel(latitude=35.1667, longitude=136.965, method='nearest')
wx.ts(gpv.time.data, gpv.data)
# 6か月アンサンブル数値予報モデルGPV
yyyymmdd = '20230601'
ele = 'Prr_Emb'
prdct = 'cps3'
grb_dir = f'./jmadata/{prdct}/{yyyymmdd[0:4]}/{yyyymmdd[0:6]}'
grb_paths = f'{grb_dir}/Z__C_RJTD_{yyyymmdd}000000_EPSC_MGPV_Rgl_Gll1p25deg_Lsurf_{ele}_grib2.bin'
ds = wx.getgpv(grb_paths, 'var discipline=0 center=34 local_table=1 parmcat=1 parm=210', timezone='Asia/Tokyo')
gpv = ds['var0_1_210_surface'].sel(latitude=35.1667, longitude=136.965, method='nearest')
wx.ts(gpv.awaretime, gpv.data)
おわりに¶
この研修では、気象庁が対外的に作成しているGPVデータを分析する際に必要となるテクニックの基礎を学習しました。
GPVデータは、Pyhtonのライブラリ Xarray が提供する Dataset および DataArray
で欠落なく情報を保持することができ、多くの処理が、メソッドや関数により簡単に実行できることが理解いただけたかと思います。Xarray には、今回は説明しなかったメソッドや関数がまだ沢山あるので、機会を捉えて習得していただきたいと思います。aware な時間を次元に持てない、グルーピングを1次元でしか行えないなど、残念なところも見られますが、開発コミュニティは大変活発なので、近い将来にこれらも解決するかもしれません。
気象庁GPVが準拠する中央協定時と日本標準時を aware な時刻でつなげる方法については、大袈裟と感じた方もいるかもしれません、時刻データを9時間遅らせればそれでよいようですが、データを他者と共有する際、それは、コーラのボトルに別な飲み物を入れて渡すのと同じ危うさを持っています。したがって、折り合いの付け方についてはよく注意して選択してください。
学習を通して、GRIBファイルからの気象庁GPVの読み込みは、関数 getgpv で簡単に行えるように感じたか方もいるかもしれませんが、getgpv は、気象庁の全てのGPVプロダクトに対して動作確認をしているわけではないので、プロダクトによっては wgrib2 にまで戻ってデーコード法を検討する必要があるかもしれません。そのようなときは、wgrib2 がもつ本来の使い方と多様なオプションを思い出していただければと思います。
最後に、この研修では、「分析チャレンジ!」と言いながら、データのハンドリング法の学習に終始し「分析」といえるようなものは行わなかったので、その点不満の残る方もいたかもしれません。「分析」はデータを自在に扱うスキルと明確な目的と対象のもとで行うものなので、高度かつ個別的にならざるを得ず、基礎編の範疇に入れられませんでした。ご了承ください。
著作権について¶
(C) 2024 WXBC
<利用条件>
本書は、本書に記載した要件・技術・方式に関する内容が変更されないこと、および出典を明示いただくことを前提に、無償でその全部または一部を複製、翻案、翻訳、転記、引用、公衆送信等して利用できます。なお、全体または一部を複製、翻案、翻訳された場合は、本書にある著作権表示および利用条件を明示してください。
<免責事項>
本書の著作権者は、本書の記載内容に関して、その正確性、商品性、利用目的への適合性等に関して保証するものではなく、特許権、著作権、その他の権利を侵害していないことを保証するものでもありません。本書の利用により生じた損害について、本書の著作権者は、法律上のいかなる責任も負いません。