PythonでGRIB2形式のGPVデータをデコードし、気象要素の変換計算などを経て、欲しいデータを取得する方法を解説します。
やり方はいくつかありますが、ここではPythonモジュールwxbcgrib
とwxparams
を使う方法を紹介します。
またwxbcgrib
内部で使われている、wgrib2
ソフトウェアとnetCDF4
モジュールについても紹介します。
MacのターミナルやWindowsのコマンドプロンプトで実行したコマンドをPythonで実行するには、subprocess
モジュールを使います。
ディレクトリやファイルの操作をするのに便利なos
モジュールとともに、まずはモジュールの読み込みをしましょう。
以下のセルを実行して下さい。
import os
import subprocess
subprocess
を使ってコマンドラインのコマンドを実行するには、以下のようにコードを書きます。
まずwgrib2
のバージョン情報を出力してみましょう。
コマンドラインで実行した時と同様に、wgrib2のバージョン情報が出力されればOKです。
#--subprocessでコマンドを実行
# Windowsの方は、wgrib2のパスを指定して実行(e.g. c:/wgrib2/wgrib2)
# Macの方は、wgrib2へのパスが通ってない場合はフルパスにして実行(e.g. /usr/local/bin/wgrib2)
rc = subprocess.run("wgrib2 --version",
shell=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True)
# 標準出力や標準エラーを出力
for line in rc.stderr.splitlines():
print(line)
for line in rc.stdout.splitlines():
print(line)
v0.2.0.8 2/2019 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, Sergey Varlamov
では次にweather
フォルダにある、数値予報モデルMSMのデータ(2020/09/05 12zベース、FT=0〜15)の概要(目録)を表示させてみましょう。
Z__C_RJTD_20200905120000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin
ダブルクォーテーションの中にコマンドを記述します。コマンドラインで実行した場合と同じ出力が得られればOKです。
### 練習問題(1)
# ダブルクォーテーション""の間にコマンドを入力して、実行して下さい。
#--subprocessでコマンドを実行
rc = subprocess.run("wgrib2 weather/Z__C_RJTD_20200905120000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin",
shell=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True)
# 標準出力や標準エラーを出力
for line in rc.stderr.splitlines():
print(line)
for line in rc.stdout.splitlines():
print(line)
1.1:0:d=2020090512:PRMSL:mean sea level:anl: 1.2:0:d=2020090512:PRES:surface:anl: 1.3:0:d=2020090512:UGRD:10 m above ground:anl: 1.4:0:d=2020090512:VGRD:10 m above ground:anl: 1.5:0:d=2020090512:TMP:1.5 m above ground:anl: 1.6:0:d=2020090512:RH:1.5 m above ground:anl: 1.7:0:d=2020090512:LCDC:surface:anl: 1.8:0:d=2020090512:MCDC:surface:anl: 1.9:0:d=2020090512:HCDC:surface:anl: 1.10:0:d=2020090512:TCDC:surface:anl: 1.11:0:d=2020090512:PRMSL:mean sea level:1 hour fcst: 1.12:0:d=2020090512:PRES:surface:1 hour fcst: 1.13:0:d=2020090512:UGRD:10 m above ground:1 hour fcst: 1.14:0:d=2020090512:VGRD:10 m above ground:1 hour fcst: 1.15:0:d=2020090512:TMP:1.5 m above ground:1 hour fcst: 1.16:0:d=2020090512:RH:1.5 m above ground:1 hour fcst: 1.17:0:d=2020090512:LCDC:surface:1 hour fcst: 1.18:0:d=2020090512:MCDC:surface:1 hour fcst: 1.19:0:d=2020090512:HCDC:surface:1 hour fcst: 1.20:0:d=2020090512:TCDC:surface:1 hour fcst: 1.21:0:d=2020090512:APCP:surface:0-1 hour acc fcst: 1.22:0:d=2020090512:DSWRF:surface:0-1 hour ave fcst: 1.23:0:d=2020090512:PRMSL:mean sea level:2 hour fcst: 1.24:0:d=2020090512:PRES:surface:2 hour fcst: 1.25:0:d=2020090512:UGRD:10 m above ground:2 hour fcst: 1.26:0:d=2020090512:VGRD:10 m above ground:2 hour fcst: 1.27:0:d=2020090512:TMP:1.5 m above ground:2 hour fcst: 1.28:0:d=2020090512:RH:1.5 m above ground:2 hour fcst: 1.29:0:d=2020090512:LCDC:surface:2 hour fcst: 1.30:0:d=2020090512:MCDC:surface:2 hour fcst: 1.31:0:d=2020090512:HCDC:surface:2 hour fcst: 1.32:0:d=2020090512:TCDC:surface:2 hour fcst: 1.33:0:d=2020090512:APCP:surface:1-2 hour acc fcst: 1.34:0:d=2020090512:DSWRF:surface:1-2 hour ave fcst: 1.35:0:d=2020090512:PRMSL:mean sea level:3 hour fcst: 1.36:0:d=2020090512:PRES:surface:3 hour fcst: 1.37:0:d=2020090512:UGRD:10 m above ground:3 hour fcst: 1.38:0:d=2020090512:VGRD:10 m above ground:3 hour fcst: 1.39:0:d=2020090512:TMP:1.5 m above ground:3 hour fcst: 1.40:0:d=2020090512:RH:1.5 m above ground:3 hour fcst: 1.41:0:d=2020090512:LCDC:surface:3 hour fcst: 1.42:0:d=2020090512:MCDC:surface:3 hour fcst: 1.43:0:d=2020090512:HCDC:surface:3 hour fcst: 1.44:0:d=2020090512:TCDC:surface:3 hour fcst: 1.45:0:d=2020090512:APCP:surface:2-3 hour acc fcst: 1.46:0:d=2020090512:DSWRF:surface:2-3 hour ave fcst: 1.47:0:d=2020090512:PRMSL:mean sea level:4 hour fcst: 1.48:0:d=2020090512:PRES:surface:4 hour fcst: 1.49:0:d=2020090512:UGRD:10 m above ground:4 hour fcst: 1.50:0:d=2020090512:VGRD:10 m above ground:4 hour fcst: 1.51:0:d=2020090512:TMP:1.5 m above ground:4 hour fcst: 1.52:0:d=2020090512:RH:1.5 m above ground:4 hour fcst: 1.53:0:d=2020090512:LCDC:surface:4 hour fcst: 1.54:0:d=2020090512:MCDC:surface:4 hour fcst: 1.55:0:d=2020090512:HCDC:surface:4 hour fcst: 1.56:0:d=2020090512:TCDC:surface:4 hour fcst: 1.57:0:d=2020090512:APCP:surface:3-4 hour acc fcst: 1.58:0:d=2020090512:DSWRF:surface:3-4 hour ave fcst: 1.59:0:d=2020090512:PRMSL:mean sea level:5 hour fcst: 1.60:0:d=2020090512:PRES:surface:5 hour fcst: 1.61:0:d=2020090512:UGRD:10 m above ground:5 hour fcst: 1.62:0:d=2020090512:VGRD:10 m above ground:5 hour fcst: 1.63:0:d=2020090512:TMP:1.5 m above ground:5 hour fcst: 1.64:0:d=2020090512:RH:1.5 m above ground:5 hour fcst: 1.65:0:d=2020090512:LCDC:surface:5 hour fcst: 1.66:0:d=2020090512:MCDC:surface:5 hour fcst: 1.67:0:d=2020090512:HCDC:surface:5 hour fcst: 1.68:0:d=2020090512:TCDC:surface:5 hour fcst: 1.69:0:d=2020090512:APCP:surface:4-5 hour acc fcst: 1.70:0:d=2020090512:DSWRF:surface:4-5 hour ave fcst: 1.71:0:d=2020090512:PRMSL:mean sea level:6 hour fcst: 1.72:0:d=2020090512:PRES:surface:6 hour fcst: 1.73:0:d=2020090512:UGRD:10 m above ground:6 hour fcst: 1.74:0:d=2020090512:VGRD:10 m above ground:6 hour fcst: 1.75:0:d=2020090512:TMP:1.5 m above ground:6 hour fcst: 1.76:0:d=2020090512:RH:1.5 m above ground:6 hour fcst: 1.77:0:d=2020090512:LCDC:surface:6 hour fcst: 1.78:0:d=2020090512:MCDC:surface:6 hour fcst: 1.79:0:d=2020090512:HCDC:surface:6 hour fcst: 1.80:0:d=2020090512:TCDC:surface:6 hour fcst: 1.81:0:d=2020090512:APCP:surface:5-6 hour acc fcst: 1.82:0:d=2020090512:DSWRF:surface:5-6 hour ave fcst: 1.83:0:d=2020090512:PRMSL:mean sea level:7 hour fcst: 1.84:0:d=2020090512:PRES:surface:7 hour fcst: 1.85:0:d=2020090512:UGRD:10 m above ground:7 hour fcst: 1.86:0:d=2020090512:VGRD:10 m above ground:7 hour fcst: 1.87:0:d=2020090512:TMP:1.5 m above ground:7 hour fcst: 1.88:0:d=2020090512:RH:1.5 m above ground:7 hour fcst: 1.89:0:d=2020090512:LCDC:surface:7 hour fcst: 1.90:0:d=2020090512:MCDC:surface:7 hour fcst: 1.91:0:d=2020090512:HCDC:surface:7 hour fcst: 1.92:0:d=2020090512:TCDC:surface:7 hour fcst: 1.93:0:d=2020090512:APCP:surface:6-7 hour acc fcst: 1.94:0:d=2020090512:DSWRF:surface:6-7 hour ave fcst: 1.95:0:d=2020090512:PRMSL:mean sea level:8 hour fcst: 1.96:0:d=2020090512:PRES:surface:8 hour fcst: 1.97:0:d=2020090512:UGRD:10 m above ground:8 hour fcst: 1.98:0:d=2020090512:VGRD:10 m above ground:8 hour fcst: 1.99:0:d=2020090512:TMP:1.5 m above ground:8 hour fcst: 1.100:0:d=2020090512:RH:1.5 m above ground:8 hour fcst: 1.101:0:d=2020090512:LCDC:surface:8 hour fcst: 1.102:0:d=2020090512:MCDC:surface:8 hour fcst: 1.103:0:d=2020090512:HCDC:surface:8 hour fcst: 1.104:0:d=2020090512:TCDC:surface:8 hour fcst: 1.105:0:d=2020090512:APCP:surface:7-8 hour acc fcst: 1.106:0:d=2020090512:DSWRF:surface:7-8 hour ave fcst: 1.107:0:d=2020090512:PRMSL:mean sea level:9 hour fcst: 1.108:0:d=2020090512:PRES:surface:9 hour fcst: 1.109:0:d=2020090512:UGRD:10 m above ground:9 hour fcst: 1.110:0:d=2020090512:VGRD:10 m above ground:9 hour fcst: 1.111:0:d=2020090512:TMP:1.5 m above ground:9 hour fcst: 1.112:0:d=2020090512:RH:1.5 m above ground:9 hour fcst: 1.113:0:d=2020090512:LCDC:surface:9 hour fcst: 1.114:0:d=2020090512:MCDC:surface:9 hour fcst: 1.115:0:d=2020090512:HCDC:surface:9 hour fcst: 1.116:0:d=2020090512:TCDC:surface:9 hour fcst: 1.117:0:d=2020090512:APCP:surface:8-9 hour acc fcst: 1.118:0:d=2020090512:DSWRF:surface:8-9 hour ave fcst: 1.119:0:d=2020090512:PRMSL:mean sea level:10 hour fcst: 1.120:0:d=2020090512:PRES:surface:10 hour fcst: 1.121:0:d=2020090512:UGRD:10 m above ground:10 hour fcst: 1.122:0:d=2020090512:VGRD:10 m above ground:10 hour fcst: 1.123:0:d=2020090512:TMP:1.5 m above ground:10 hour fcst: 1.124:0:d=2020090512:RH:1.5 m above ground:10 hour fcst: 1.125:0:d=2020090512:LCDC:surface:10 hour fcst: 1.126:0:d=2020090512:MCDC:surface:10 hour fcst: 1.127:0:d=2020090512:HCDC:surface:10 hour fcst: 1.128:0:d=2020090512:TCDC:surface:10 hour fcst: 1.129:0:d=2020090512:APCP:surface:9-10 hour acc fcst: 1.130:0:d=2020090512:DSWRF:surface:9-10 hour ave fcst: 1.131:0:d=2020090512:PRMSL:mean sea level:11 hour fcst: 1.132:0:d=2020090512:PRES:surface:11 hour fcst: 1.133:0:d=2020090512:UGRD:10 m above ground:11 hour fcst: 1.134:0:d=2020090512:VGRD:10 m above ground:11 hour fcst: 1.135:0:d=2020090512:TMP:1.5 m above ground:11 hour fcst: 1.136:0:d=2020090512:RH:1.5 m above ground:11 hour fcst: 1.137:0:d=2020090512:LCDC:surface:11 hour fcst: 1.138:0:d=2020090512:MCDC:surface:11 hour fcst: 1.139:0:d=2020090512:HCDC:surface:11 hour fcst: 1.140:0:d=2020090512:TCDC:surface:11 hour fcst: 1.141:0:d=2020090512:APCP:surface:10-11 hour acc fcst: 1.142:0:d=2020090512:DSWRF:surface:10-11 hour ave fcst: 1.143:0:d=2020090512:PRMSL:mean sea level:12 hour fcst: 1.144:0:d=2020090512:PRES:surface:12 hour fcst: 1.145:0:d=2020090512:UGRD:10 m above ground:12 hour fcst: 1.146:0:d=2020090512:VGRD:10 m above ground:12 hour fcst: 1.147:0:d=2020090512:TMP:1.5 m above ground:12 hour fcst: 1.148:0:d=2020090512:RH:1.5 m above ground:12 hour fcst: 1.149:0:d=2020090512:LCDC:surface:12 hour fcst: 1.150:0:d=2020090512:MCDC:surface:12 hour fcst: 1.151:0:d=2020090512:HCDC:surface:12 hour fcst: 1.152:0:d=2020090512:TCDC:surface:12 hour fcst: 1.153:0:d=2020090512:APCP:surface:11-12 hour acc fcst: 1.154:0:d=2020090512:DSWRF:surface:11-12 hour ave fcst: 1.155:0:d=2020090512:PRMSL:mean sea level:13 hour fcst: 1.156:0:d=2020090512:PRES:surface:13 hour fcst: 1.157:0:d=2020090512:UGRD:10 m above ground:13 hour fcst: 1.158:0:d=2020090512:VGRD:10 m above ground:13 hour fcst: 1.159:0:d=2020090512:TMP:1.5 m above ground:13 hour fcst: 1.160:0:d=2020090512:RH:1.5 m above ground:13 hour fcst: 1.161:0:d=2020090512:LCDC:surface:13 hour fcst: 1.162:0:d=2020090512:MCDC:surface:13 hour fcst: 1.163:0:d=2020090512:HCDC:surface:13 hour fcst: 1.164:0:d=2020090512:TCDC:surface:13 hour fcst: 1.165:0:d=2020090512:APCP:surface:12-13 hour acc fcst: 1.166:0:d=2020090512:DSWRF:surface:12-13 hour ave fcst: 1.167:0:d=2020090512:PRMSL:mean sea level:14 hour fcst: 1.168:0:d=2020090512:PRES:surface:14 hour fcst: 1.169:0:d=2020090512:UGRD:10 m above ground:14 hour fcst: 1.170:0:d=2020090512:VGRD:10 m above ground:14 hour fcst: 1.171:0:d=2020090512:TMP:1.5 m above ground:14 hour fcst: 1.172:0:d=2020090512:RH:1.5 m above ground:14 hour fcst: 1.173:0:d=2020090512:LCDC:surface:14 hour fcst: 1.174:0:d=2020090512:MCDC:surface:14 hour fcst: 1.175:0:d=2020090512:HCDC:surface:14 hour fcst: 1.176:0:d=2020090512:TCDC:surface:14 hour fcst: 1.177:0:d=2020090512:APCP:surface:13-14 hour acc fcst: 1.178:0:d=2020090512:DSWRF:surface:13-14 hour ave fcst: 1.179:0:d=2020090512:PRMSL:mean sea level:15 hour fcst: 1.180:0:d=2020090512:PRES:surface:15 hour fcst: 1.181:0:d=2020090512:UGRD:10 m above ground:15 hour fcst: 1.182:0:d=2020090512:VGRD:10 m above ground:15 hour fcst: 1.183:0:d=2020090512:TMP:1.5 m above ground:15 hour fcst: 1.184:0:d=2020090512:RH:1.5 m above ground:15 hour fcst: 1.185:0:d=2020090512:LCDC:surface:15 hour fcst: 1.186:0:d=2020090512:MCDC:surface:15 hour fcst: 1.187:0:d=2020090512:HCDC:surface:15 hour fcst: 1.188:0:d=2020090512:TCDC:surface:15 hour fcst: 1.189:0:d=2020090512:APCP:surface:14-15 hour acc fcst: 1.190:0:d=2020090512:DSWRF:surface:14-15 hour ave fcst:
それでは次は、wgrib2
のmatch
オプションを使って183番目の要素を、V
オプションを使ってその詳細情報を出力してみましょう。
オプションが増えるとコマンドが長くなりますので、下記のように適宜変数を活用するとコードが見やすくなります。
ここでf文字列
と呼ばれる機能を使い、変数からコマンドの文字列を構成しています。
#--subprocessでコマンドを実行
wgrib2 = "wgrib2"
gpvdir = "weather"
gpvname = "Z__C_RJTD_20200905120000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
gpvpath = os.path.join(gpvdir, gpvname)
option = "-match '1.183:'"
# f文字列を使ってコマンドを作成
print("実行するコマンド")
print(f"{wgrib2} {gpvpath} {option} -V")
print("")
rc = subprocess.run(f"{wgrib2} {gpvpath} {option} -V",
shell=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True)
# 標準出力や標準エラーを出力
for line in rc.stderr.splitlines():
print(line)
for line in rc.stdout.splitlines():
print(line)
実行するコマンド wgrib2 weather/Z__C_RJTD_20200905120000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin -match '1.183:' -V 1.183:0:vt=2020090603:1.5 m above ground:15 hour fcst:TMP Temperature [K]: ndata=242905:undef=0:mean=298.983:min=283.563:max=306.165 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
だいぶ慣れてきましたでしょうか?
ではいよいよGRIB2形式のMSMデータを、NetCDF形式に変換しましょう。
wgrib2
にnetcdf
オプションを、以下のように与えて実行します。
wgrib2 <GRIB2ファイル> -netcdf <出力するNetCDFファイル名>
ここではGRIB2ファイル名.nc
という命名則でNetCDFファイルを作るようにします。
変数を活用して、ダブルクォーテーションの中にコマンドを記述し、セルを実行して下さい。
### 練習問題(2)
# ダブルクォーテーション""の間にコマンドを入力して、実行して下さい。
#--subprocessでコマンドを実行
# 変数wgrib2, 変数gpvpathは上で定義済み
ncname = gpvname + ".nc"
ncpath = os.path.join(".", ncname)
rc = subprocess.run(f"{wgrib2} {gpvpath} -netcdf {ncpath}",
shell=True,
stdout=subprocess.PIPE,
stderr=subprocess.PIPE,
universal_newlines=True)
# 標準出力や標準エラーを出力
for line in rc.stderr.splitlines():
print(line)
for line in rc.stdout.splitlines():
print(line)
1.1:0:d=2020090512:PRMSL:mean sea level:anl: 1.2:0:d=2020090512:PRES:surface:anl: 1.3:0:d=2020090512:UGRD:10 m above ground:anl: 1.4:0:d=2020090512:VGRD:10 m above ground:anl: 1.5:0:d=2020090512:TMP:1.5 m above ground:anl: 1.6:0:d=2020090512:RH:1.5 m above ground:anl: 1.7:0:d=2020090512:LCDC:surface:anl: 1.8:0:d=2020090512:MCDC:surface:anl: 1.9:0:d=2020090512:HCDC:surface:anl: 1.10:0:d=2020090512:TCDC:surface:anl: 1.11:0:d=2020090512:PRMSL:mean sea level:1 hour fcst: 1.12:0:d=2020090512:PRES:surface:1 hour fcst: 1.13:0:d=2020090512:UGRD:10 m above ground:1 hour fcst: 1.14:0:d=2020090512:VGRD:10 m above ground:1 hour fcst: 1.15:0:d=2020090512:TMP:1.5 m above ground:1 hour fcst: 1.16:0:d=2020090512:RH:1.5 m above ground:1 hour fcst: 1.17:0:d=2020090512:LCDC:surface:1 hour fcst: 1.18:0:d=2020090512:MCDC:surface:1 hour fcst: 1.19:0:d=2020090512:HCDC:surface:1 hour fcst: 1.20:0:d=2020090512:TCDC:surface:1 hour fcst: 1.21:0:d=2020090512:APCP:surface:0-1 hour acc fcst: 1.22:0:d=2020090512:DSWRF:surface:0-1 hour ave fcst: 1.23:0:d=2020090512:PRMSL:mean sea level:2 hour fcst: 1.24:0:d=2020090512:PRES:surface:2 hour fcst: 1.25:0:d=2020090512:UGRD:10 m above ground:2 hour fcst: 1.26:0:d=2020090512:VGRD:10 m above ground:2 hour fcst: 1.27:0:d=2020090512:TMP:1.5 m above ground:2 hour fcst: 1.28:0:d=2020090512:RH:1.5 m above ground:2 hour fcst: 1.29:0:d=2020090512:LCDC:surface:2 hour fcst: 1.30:0:d=2020090512:MCDC:surface:2 hour fcst: 1.31:0:d=2020090512:HCDC:surface:2 hour fcst: 1.32:0:d=2020090512:TCDC:surface:2 hour fcst: 1.33:0:d=2020090512:APCP:surface:1-2 hour acc fcst: 1.34:0:d=2020090512:DSWRF:surface:1-2 hour ave fcst: 1.35:0:d=2020090512:PRMSL:mean sea level:3 hour fcst: 1.36:0:d=2020090512:PRES:surface:3 hour fcst: 1.37:0:d=2020090512:UGRD:10 m above ground:3 hour fcst: 1.38:0:d=2020090512:VGRD:10 m above ground:3 hour fcst: 1.39:0:d=2020090512:TMP:1.5 m above ground:3 hour fcst: 1.40:0:d=2020090512:RH:1.5 m above ground:3 hour fcst: 1.41:0:d=2020090512:LCDC:surface:3 hour fcst: 1.42:0:d=2020090512:MCDC:surface:3 hour fcst: 1.43:0:d=2020090512:HCDC:surface:3 hour fcst: 1.44:0:d=2020090512:TCDC:surface:3 hour fcst: 1.45:0:d=2020090512:APCP:surface:2-3 hour acc fcst: 1.46:0:d=2020090512:DSWRF:surface:2-3 hour ave fcst: 1.47:0:d=2020090512:PRMSL:mean sea level:4 hour fcst: 1.48:0:d=2020090512:PRES:surface:4 hour fcst: 1.49:0:d=2020090512:UGRD:10 m above ground:4 hour fcst: 1.50:0:d=2020090512:VGRD:10 m above ground:4 hour fcst: 1.51:0:d=2020090512:TMP:1.5 m above ground:4 hour fcst: 1.52:0:d=2020090512:RH:1.5 m above ground:4 hour fcst: 1.53:0:d=2020090512:LCDC:surface:4 hour fcst: 1.54:0:d=2020090512:MCDC:surface:4 hour fcst: 1.55:0:d=2020090512:HCDC:surface:4 hour fcst: 1.56:0:d=2020090512:TCDC:surface:4 hour fcst: 1.57:0:d=2020090512:APCP:surface:3-4 hour acc fcst: 1.58:0:d=2020090512:DSWRF:surface:3-4 hour ave fcst: 1.59:0:d=2020090512:PRMSL:mean sea level:5 hour fcst: 1.60:0:d=2020090512:PRES:surface:5 hour fcst: 1.61:0:d=2020090512:UGRD:10 m above ground:5 hour fcst: 1.62:0:d=2020090512:VGRD:10 m above ground:5 hour fcst: 1.63:0:d=2020090512:TMP:1.5 m above ground:5 hour fcst: 1.64:0:d=2020090512:RH:1.5 m above ground:5 hour fcst: 1.65:0:d=2020090512:LCDC:surface:5 hour fcst: 1.66:0:d=2020090512:MCDC:surface:5 hour fcst: 1.67:0:d=2020090512:HCDC:surface:5 hour fcst: 1.68:0:d=2020090512:TCDC:surface:5 hour fcst: 1.69:0:d=2020090512:APCP:surface:4-5 hour acc fcst: 1.70:0:d=2020090512:DSWRF:surface:4-5 hour ave fcst: 1.71:0:d=2020090512:PRMSL:mean sea level:6 hour fcst: 1.72:0:d=2020090512:PRES:surface:6 hour fcst: 1.73:0:d=2020090512:UGRD:10 m above ground:6 hour fcst: 1.74:0:d=2020090512:VGRD:10 m above ground:6 hour fcst: 1.75:0:d=2020090512:TMP:1.5 m above ground:6 hour fcst: 1.76:0:d=2020090512:RH:1.5 m above ground:6 hour fcst: 1.77:0:d=2020090512:LCDC:surface:6 hour fcst: 1.78:0:d=2020090512:MCDC:surface:6 hour fcst: 1.79:0:d=2020090512:HCDC:surface:6 hour fcst: 1.80:0:d=2020090512:TCDC:surface:6 hour fcst: 1.81:0:d=2020090512:APCP:surface:5-6 hour acc fcst: 1.82:0:d=2020090512:DSWRF:surface:5-6 hour ave fcst: 1.83:0:d=2020090512:PRMSL:mean sea level:7 hour fcst: 1.84:0:d=2020090512:PRES:surface:7 hour fcst: 1.85:0:d=2020090512:UGRD:10 m above ground:7 hour fcst: 1.86:0:d=2020090512:VGRD:10 m above ground:7 hour fcst: 1.87:0:d=2020090512:TMP:1.5 m above ground:7 hour fcst: 1.88:0:d=2020090512:RH:1.5 m above ground:7 hour fcst: 1.89:0:d=2020090512:LCDC:surface:7 hour fcst: 1.90:0:d=2020090512:MCDC:surface:7 hour fcst: 1.91:0:d=2020090512:HCDC:surface:7 hour fcst: 1.92:0:d=2020090512:TCDC:surface:7 hour fcst: 1.93:0:d=2020090512:APCP:surface:6-7 hour acc fcst: 1.94:0:d=2020090512:DSWRF:surface:6-7 hour ave fcst: 1.95:0:d=2020090512:PRMSL:mean sea level:8 hour fcst: 1.96:0:d=2020090512:PRES:surface:8 hour fcst: 1.97:0:d=2020090512:UGRD:10 m above ground:8 hour fcst: 1.98:0:d=2020090512:VGRD:10 m above ground:8 hour fcst: 1.99:0:d=2020090512:TMP:1.5 m above ground:8 hour fcst: 1.100:0:d=2020090512:RH:1.5 m above ground:8 hour fcst: 1.101:0:d=2020090512:LCDC:surface:8 hour fcst: 1.102:0:d=2020090512:MCDC:surface:8 hour fcst: 1.103:0:d=2020090512:HCDC:surface:8 hour fcst: 1.104:0:d=2020090512:TCDC:surface:8 hour fcst: 1.105:0:d=2020090512:APCP:surface:7-8 hour acc fcst: 1.106:0:d=2020090512:DSWRF:surface:7-8 hour ave fcst: 1.107:0:d=2020090512:PRMSL:mean sea level:9 hour fcst: 1.108:0:d=2020090512:PRES:surface:9 hour fcst: 1.109:0:d=2020090512:UGRD:10 m above ground:9 hour fcst: 1.110:0:d=2020090512:VGRD:10 m above ground:9 hour fcst: 1.111:0:d=2020090512:TMP:1.5 m above ground:9 hour fcst: 1.112:0:d=2020090512:RH:1.5 m above ground:9 hour fcst: 1.113:0:d=2020090512:LCDC:surface:9 hour fcst: 1.114:0:d=2020090512:MCDC:surface:9 hour fcst: 1.115:0:d=2020090512:HCDC:surface:9 hour fcst: 1.116:0:d=2020090512:TCDC:surface:9 hour fcst: 1.117:0:d=2020090512:APCP:surface:8-9 hour acc fcst: 1.118:0:d=2020090512:DSWRF:surface:8-9 hour ave fcst: 1.119:0:d=2020090512:PRMSL:mean sea level:10 hour fcst: 1.120:0:d=2020090512:PRES:surface:10 hour fcst: 1.121:0:d=2020090512:UGRD:10 m above ground:10 hour fcst: 1.122:0:d=2020090512:VGRD:10 m above ground:10 hour fcst: 1.123:0:d=2020090512:TMP:1.5 m above ground:10 hour fcst: 1.124:0:d=2020090512:RH:1.5 m above ground:10 hour fcst: 1.125:0:d=2020090512:LCDC:surface:10 hour fcst: 1.126:0:d=2020090512:MCDC:surface:10 hour fcst: 1.127:0:d=2020090512:HCDC:surface:10 hour fcst: 1.128:0:d=2020090512:TCDC:surface:10 hour fcst: 1.129:0:d=2020090512:APCP:surface:9-10 hour acc fcst: 1.130:0:d=2020090512:DSWRF:surface:9-10 hour ave fcst: 1.131:0:d=2020090512:PRMSL:mean sea level:11 hour fcst: 1.132:0:d=2020090512:PRES:surface:11 hour fcst: 1.133:0:d=2020090512:UGRD:10 m above ground:11 hour fcst: 1.134:0:d=2020090512:VGRD:10 m above ground:11 hour fcst: 1.135:0:d=2020090512:TMP:1.5 m above ground:11 hour fcst: 1.136:0:d=2020090512:RH:1.5 m above ground:11 hour fcst: 1.137:0:d=2020090512:LCDC:surface:11 hour fcst: 1.138:0:d=2020090512:MCDC:surface:11 hour fcst: 1.139:0:d=2020090512:HCDC:surface:11 hour fcst: 1.140:0:d=2020090512:TCDC:surface:11 hour fcst: 1.141:0:d=2020090512:APCP:surface:10-11 hour acc fcst: 1.142:0:d=2020090512:DSWRF:surface:10-11 hour ave fcst: 1.143:0:d=2020090512:PRMSL:mean sea level:12 hour fcst: 1.144:0:d=2020090512:PRES:surface:12 hour fcst: 1.145:0:d=2020090512:UGRD:10 m above ground:12 hour fcst: 1.146:0:d=2020090512:VGRD:10 m above ground:12 hour fcst: 1.147:0:d=2020090512:TMP:1.5 m above ground:12 hour fcst: 1.148:0:d=2020090512:RH:1.5 m above ground:12 hour fcst: 1.149:0:d=2020090512:LCDC:surface:12 hour fcst: 1.150:0:d=2020090512:MCDC:surface:12 hour fcst: 1.151:0:d=2020090512:HCDC:surface:12 hour fcst: 1.152:0:d=2020090512:TCDC:surface:12 hour fcst: 1.153:0:d=2020090512:APCP:surface:11-12 hour acc fcst: 1.154:0:d=2020090512:DSWRF:surface:11-12 hour ave fcst: 1.155:0:d=2020090512:PRMSL:mean sea level:13 hour fcst: 1.156:0:d=2020090512:PRES:surface:13 hour fcst: 1.157:0:d=2020090512:UGRD:10 m above ground:13 hour fcst: 1.158:0:d=2020090512:VGRD:10 m above ground:13 hour fcst: 1.159:0:d=2020090512:TMP:1.5 m above ground:13 hour fcst: 1.160:0:d=2020090512:RH:1.5 m above ground:13 hour fcst: 1.161:0:d=2020090512:LCDC:surface:13 hour fcst: 1.162:0:d=2020090512:MCDC:surface:13 hour fcst: 1.163:0:d=2020090512:HCDC:surface:13 hour fcst: 1.164:0:d=2020090512:TCDC:surface:13 hour fcst: 1.165:0:d=2020090512:APCP:surface:12-13 hour acc fcst: 1.166:0:d=2020090512:DSWRF:surface:12-13 hour ave fcst: 1.167:0:d=2020090512:PRMSL:mean sea level:14 hour fcst: 1.168:0:d=2020090512:PRES:surface:14 hour fcst: 1.169:0:d=2020090512:UGRD:10 m above ground:14 hour fcst: 1.170:0:d=2020090512:VGRD:10 m above ground:14 hour fcst: 1.171:0:d=2020090512:TMP:1.5 m above ground:14 hour fcst: 1.172:0:d=2020090512:RH:1.5 m above ground:14 hour fcst: 1.173:0:d=2020090512:LCDC:surface:14 hour fcst: 1.174:0:d=2020090512:MCDC:surface:14 hour fcst: 1.175:0:d=2020090512:HCDC:surface:14 hour fcst: 1.176:0:d=2020090512:TCDC:surface:14 hour fcst: 1.177:0:d=2020090512:APCP:surface:13-14 hour acc fcst: 1.178:0:d=2020090512:DSWRF:surface:13-14 hour ave fcst: 1.179:0:d=2020090512:PRMSL:mean sea level:15 hour fcst: 1.180:0:d=2020090512:PRES:surface:15 hour fcst: 1.181:0:d=2020090512:UGRD:10 m above ground:15 hour fcst: 1.182:0:d=2020090512:VGRD:10 m above ground:15 hour fcst: 1.183:0:d=2020090512:TMP:1.5 m above ground:15 hour fcst: 1.184:0:d=2020090512:RH:1.5 m above ground:15 hour fcst: 1.185:0:d=2020090512:LCDC:surface:15 hour fcst: 1.186:0:d=2020090512:MCDC:surface:15 hour fcst: 1.187:0:d=2020090512:HCDC:surface:15 hour fcst: 1.188:0:d=2020090512:TCDC:surface:15 hour fcst: 1.189:0:d=2020090512:APCP:surface:14-15 hour acc fcst: 1.190:0:d=2020090512:DSWRF:surface:14-15 hour ave fcst:
このipynbファイルがあるフォルダを開き、NetCDFファイルが作成されていることを確認して下さい。
ここでファイルサイズを見て下さい。元々のGRIB2が約70MB、NetCDFが約180MBほどではないでしょうか。
一般にはGRIB2の方が圧縮性が高いため、ファイルサイズも小さくなっています。
以上で、Python上でwgrib2
を実行し、GRIB2のGPVデータをNetCDFへ変換することができました。
目次へ戻る
NetCDF形式のファイルの読み込みには、その名もnetCDF4
モジュールを利用します。
また読み込んだデータを扱う際にはnumpy
も利用します。
まずはモジュールのimportから行います。以下のセルを実行して下さい。
import netCDF4
import numpy as np
netCDF4
のDataset
クラスを使ってNetCDFファイルを読み込みます。
以下のセルを実行するだけで、読み込みは完了です。
#--NetCDFファイルの読み込み
# 第1引数はファイルへのパス
# → 第1章で定義した変数を使い、第1章で作成したNetCDFファイルを読み込む
# 第2引数は読み込み専用オプションを表す
ds = netCDF4.Dataset(ncpath, "r")
変数ds
には、Datasetクラスにより作成されたオブジェクトが入っています。
NetCDFファイルから読み込んだデータの要素を見てみましょう。以下のセルを実行して下さい。
# データの要素名を出力
ds.variables.keys()
dict_keys(['latitude', 'longitude', 'time', 'PRMSL_meansealevel', 'PRES_surface', 'UGRD_10maboveground', 'VGRD_10maboveground', 'TMP_1D5maboveground', 'RH_1D5maboveground', 'LCDC_surface', 'MCDC_surface', 'HCDC_surface', 'TCDC_surface', 'APCP_surface', 'DSWRF_surface'])
緯度経度、予想対象時刻、各種気象要素が入っているようです。
気象要素は『要素名_高度情報』という命名則になっています。
ではこの中から気温を選択してみましょう。以下のセルを実行して下さい。
ds.variables['TMP_1D5maboveground']
<class 'netCDF4._netCDF4.Variable'> float32 TMP_1D5maboveground(time, latitude, longitude) _FillValue: 9.999e+20 short_name: TMP_1D5maboveground long_name: Temperature level: 1.5 m above ground units: K unlimited dimensions: time current shape = (16, 505, 481) filling on
何やら気温のデータに関するヘッダー情報のようなものが表示されました。
これらの情報は、以下のようにしてもう少し細かく分解できます。
# 気象要素名と単位を表示
print("データ名:", ds['TMP_1D5maboveground'].long_name)
print("単位:", ds['TMP_1D5maboveground'].units)
データ名: Temperature 単位: K
気温の値そのものには、以下のようにしてアクセスできます。
# 気温の値を表示
ds['TMP_1D5maboveground'][:]
masked_array( data=[[[302.0578 , 302.03436, 302.01874, ..., 302.08905, 302.04218, 302.0031 ], [302.02655, 301.9875 , 301.9406 , ..., 302.0656 , 302.01874, 301.9875 ], [301.9953 , 301.9406 , 301.8781 , ..., 302.05 , 302.01093, 301.97968], ..., [283.0031 , 282.925 , 282.9328 , ..., 286.76874, 286.76874, 286.76093], [282.88593, 282.78436, 282.76093, ..., 286.8078 , 286.8 , 286.79218], [283.02655, 282.8703 , 282.8156 , ..., 286.83124, 286.82343, 286.82343]], [[301.89423, 301.80048, 301.6911 , ..., 302.04266, 302.0114 , 301.98798], [301.86298, 301.7458 , 301.6208 , ..., 302.03485, 301.9958 , 301.97235], [301.8239 , 301.6911 , 301.53485, ..., 302.01923, 301.98798, 301.95673], ..., [282.53485, 282.34735, 282.33954, ..., 286.83954, 286.8239 , 286.8239 ], [282.4411 , 282.2458 , 282.1989 , ..., 286.8708 , 286.86298, 286.86298], [282.6208 , 282.4411 , 282.3083 , ..., 286.89423, 286.8864 , 286.8864 ]], [[301.6412 , 301.46152, 301.28964, ..., 302.0162 , 301.98495, 301.9537 ], [301.59433, 301.39902, 301.1959 , ..., 302.00058, 301.96933, 301.9459 ], [301.55527, 301.33652, 301.11777, ..., 301.99277, 301.96152, 301.93027], ..., [282.27402, 281.97714, 281.90683, ..., 286.96152, 286.9459 , 286.93027], [282.33652, 281.99277, 281.85214, ..., 286.99277, 286.96933, 286.9537 ], [282.4459 , 282.2584 , 281.99277, ..., 287.00058, 286.97714, 286.97714]], ..., [[300.76288, 300.70038, 300.6457 , ..., 302.03632, 301.99725, 301.966 ], [300.72382, 300.6535 , 300.59882, ..., 302.04413, 302.00507, 301.97382], [300.67694, 300.60663, 300.55194, ..., 302.05194, 302.01288, 301.98163], ..., [287.5832 , 287.5832 , 287.69257, ..., 287.3332 , 287.36444, 287.41913], [287.4582 , 287.30194, 287.43475, ..., 287.41132, 287.42694, 287.466 ], [287.4035 , 287.2707 , 287.31757, ..., 287.48944, 287.48163, 287.50507]], [[300.9903 , 300.95905, 300.9356 , ..., 302.0528 , 302.01373, 301.98248], [300.96686, 300.9356 , 300.91217, ..., 302.0606 , 302.02155, 301.9903 ], [300.95905, 300.91998, 300.89655, ..., 302.06842, 302.02936, 301.9981 ], ..., [289.1153 , 289.20123, 289.41998, ..., 287.3106 , 287.34186, 287.38873], [288.88092, 288.7481 , 288.95123, ..., 287.38873, 287.40436, 287.44342], [288.77936, 288.63873, 288.7403 , ..., 287.47467, 287.46686, 287.4903 ]], [[301.36017, 301.3133 , 301.3133 , ..., 302.05548, 302.01642, 301.97736], [301.3758 , 301.33673, 301.33673, ..., 302.0633 , 302.01642, 301.98517], [301.3836 , 301.35236, 301.3758 , ..., 302.0711 , 302.02423, 301.98517], ..., [290.26642, 290.36798, 290.61798, ..., 287.29767, 287.32892, 287.3758 ], [290.01642, 289.89142, 290.11798, ..., 287.36017, 287.39142, 287.4383 ], [289.89142, 289.7586 , 289.8758 , ..., 287.42267, 287.4383 , 287.47736]]], mask=False, fill_value=1e+20, dtype=float32)
見るとmasked_array
という特殊なフォーマットで値が格納されていることがわかります。
masked_array
とは、欠損値を含む場合に欠損箇所の情報を持たせることができるデータ型ですmasked_array
が有効です以下のようにすると、通常のnumpy.ndarray型で値を取り出すことができます。
# 気温の値を変数tempに格納
temp = ds['TMP_1D5maboveground'][:].data
print("データサイズ:", temp.shape)
temp
データサイズ: (16, 505, 481)
array([[[302.0578 , 302.03436, 302.01874, ..., 302.08905, 302.04218, 302.0031 ], [302.02655, 301.9875 , 301.9406 , ..., 302.0656 , 302.01874, 301.9875 ], [301.9953 , 301.9406 , 301.8781 , ..., 302.05 , 302.01093, 301.97968], ..., [283.0031 , 282.925 , 282.9328 , ..., 286.76874, 286.76874, 286.76093], [282.88593, 282.78436, 282.76093, ..., 286.8078 , 286.8 , 286.79218], [283.02655, 282.8703 , 282.8156 , ..., 286.83124, 286.82343, 286.82343]], [[301.89423, 301.80048, 301.6911 , ..., 302.04266, 302.0114 , 301.98798], [301.86298, 301.7458 , 301.6208 , ..., 302.03485, 301.9958 , 301.97235], [301.8239 , 301.6911 , 301.53485, ..., 302.01923, 301.98798, 301.95673], ..., [282.53485, 282.34735, 282.33954, ..., 286.83954, 286.8239 , 286.8239 ], [282.4411 , 282.2458 , 282.1989 , ..., 286.8708 , 286.86298, 286.86298], [282.6208 , 282.4411 , 282.3083 , ..., 286.89423, 286.8864 , 286.8864 ]], [[301.6412 , 301.46152, 301.28964, ..., 302.0162 , 301.98495, 301.9537 ], [301.59433, 301.39902, 301.1959 , ..., 302.00058, 301.96933, 301.9459 ], [301.55527, 301.33652, 301.11777, ..., 301.99277, 301.96152, 301.93027], ..., [282.27402, 281.97714, 281.90683, ..., 286.96152, 286.9459 , 286.93027], [282.33652, 281.99277, 281.85214, ..., 286.99277, 286.96933, 286.9537 ], [282.4459 , 282.2584 , 281.99277, ..., 287.00058, 286.97714, 286.97714]], ..., [[300.76288, 300.70038, 300.6457 , ..., 302.03632, 301.99725, 301.966 ], [300.72382, 300.6535 , 300.59882, ..., 302.04413, 302.00507, 301.97382], [300.67694, 300.60663, 300.55194, ..., 302.05194, 302.01288, 301.98163], ..., [287.5832 , 287.5832 , 287.69257, ..., 287.3332 , 287.36444, 287.41913], [287.4582 , 287.30194, 287.43475, ..., 287.41132, 287.42694, 287.466 ], [287.4035 , 287.2707 , 287.31757, ..., 287.48944, 287.48163, 287.50507]], [[300.9903 , 300.95905, 300.9356 , ..., 302.0528 , 302.01373, 301.98248], [300.96686, 300.9356 , 300.91217, ..., 302.0606 , 302.02155, 301.9903 ], [300.95905, 300.91998, 300.89655, ..., 302.06842, 302.02936, 301.9981 ], ..., [289.1153 , 289.20123, 289.41998, ..., 287.3106 , 287.34186, 287.38873], [288.88092, 288.7481 , 288.95123, ..., 287.38873, 287.40436, 287.44342], [288.77936, 288.63873, 288.7403 , ..., 287.47467, 287.46686, 287.4903 ]], [[301.36017, 301.3133 , 301.3133 , ..., 302.05548, 302.01642, 301.97736], [301.3758 , 301.33673, 301.33673, ..., 302.0633 , 302.01642, 301.98517], [301.3836 , 301.35236, 301.3758 , ..., 302.0711 , 302.02423, 301.98517], ..., [290.26642, 290.36798, 290.61798, ..., 287.29767, 287.32892, 287.3758 ], [290.01642, 289.89142, 290.11798, ..., 287.36017, 287.39142, 287.4383 ], [289.89142, 289.7586 , 289.8758 , ..., 287.42267, 287.4383 , 287.47736]]], dtype=float32)
気温の値は、そのオーダーからも絶対温度[K]であることがわかります。
またデータサイズは3次元になっており、(FT, Latitude, Longitude)
の数を表します。
これで気象要素を代表して気温の値を取り出すことができました。
では次は地上気圧のデータを取得してみましょう。
### 練習問題(3)
# 地上気圧のデータをnumpy.ndarray型で取得して下さい
pres = ds.variables['PRES_surface'][:].data
pres
array([[[100515.46 , 100515.46 , 100515.46 , ..., 101227.96 , 101227.96 , 101227.96 ], [100515.46 , 100515.46 , 100515.46 , ..., 101227.96 , 101227.96 , 101227.96 ], [100515.46 , 100515.46 , 100515.46 , ..., 101227.96 , 101227.96 , 101227.96 ], ..., [ 87502.96 , 87427.96 , 87565.46 , ..., 101927.96 , 101927.96 , 101927.96 ], [ 87340.46 , 87065.46 , 87165.46 , ..., 101927.96 , 101927.96 , 101927.96 ], [ 87402.96 , 87077.96 , 87102.96 , ..., 101927.96 , 101927.96 , 101927.96 ]], [[100562.16 , 100562.16 , 100562.16 , ..., 101199.66 , 101199.66 , 101199.66 ], [100562.16 , 100562.16 , 100562.16 , ..., 101199.66 , 101199.66 , 101199.66 ], [100562.16 , 100562.16 , 100562.16 , ..., 101212.16 , 101212.16 , 101212.16 ], ..., [ 87524.66 , 87449.66 , 87587.16 , ..., 101899.66 , 101912.16 , 101912.16 ], [ 87362.16 , 87087.16 , 87187.16 , ..., 101899.66 , 101912.16 , 101912.16 ], [ 87424.66 , 87112.16 , 87137.16 , ..., 101899.66 , 101899.66 , 101912.16 ]], [[100555.836, 100555.836, 100543.336, ..., 101143.336, 101155.836, 101155.836], [100555.836, 100555.836, 100543.336, ..., 101155.836, 101155.836, 101155.836], [100555.836, 100555.836, 100555.836, ..., 101155.836, 101155.836, 101155.836], ..., [ 87543.336, 87468.336, 87605.836, ..., 101905.836, 101905.836, 101905.836], [ 87380.836, 87105.836, 87205.836, ..., 101893.336, 101905.836, 101905.836], [ 87443.336, 87118.336, 87143.336, ..., 101893.336, 101893.336, 101905.836]], ..., [[100581.59 , 100569.09 , 100569.09 , ..., 101269.09 , 101269.09 , 101269.09 ], [100581.59 , 100569.09 , 100569.09 , ..., 101269.09 , 101269.09 , 101269.09 ], [100581.59 , 100569.09 , 100569.09 , ..., 101269.09 , 101269.09 , 101269.09 ], ..., [ 87794.09 , 87719.09 , 87856.59 , ..., 101806.59 , 101819.09 , 101819.09 ], [ 87631.59 , 87369.09 , 87456.59 , ..., 101806.59 , 101806.59 , 101819.09 ], [ 87706.59 , 87381.59 , 87406.59 , ..., 101806.59 , 101806.59 , 101806.59 ]], [[100562.49 , 100562.49 , 100549.99 , ..., 101224.99 , 101224.99 , 101237.49 ], [100562.49 , 100562.49 , 100549.99 , ..., 101237.49 , 101237.49 , 101237.49 ], [100562.49 , 100562.49 , 100549.99 , ..., 101237.49 , 101237.49 , 101237.49 ], ..., [ 87812.49 , 87737.49 , 87862.49 , ..., 101787.49 , 101787.49 , 101799.99 ], [ 87649.99 , 87374.99 , 87474.99 , ..., 101787.49 , 101787.49 , 101787.49 ], [ 87724.99 , 87399.99 , 87424.99 , ..., 101787.49 , 101787.49 , 101787.49 ]], [[100535.56 , 100523.06 , 100510.56 , ..., 101185.56 , 101185.56 , 101185.56 ], [100535.56 , 100523.06 , 100510.56 , ..., 101185.56 , 101185.56 , 101185.56 ], [100535.56 , 100523.06 , 100510.56 , ..., 101185.56 , 101185.56 , 101185.56 ], ..., [ 87810.56 , 87735.56 , 87873.06 , ..., 101760.56 , 101773.06 , 101773.06 ], [ 87648.06 , 87385.56 , 87473.06 , ..., 101760.56 , 101760.56 , 101760.56 ], [ 87723.06 , 87398.06 , 87423.06 , ..., 101760.56 , 101760.56 , 101760.56 ]]], dtype=float32)
以上で、netCDF4
モジュールを使ってNetCDF形式のGPVデータを読み込むことができました。
目次へ戻る
第1章・第2章ではwgrib2
を使ってGRIB2形式のGPVをNetCDFに変換し、netCDF4
モジュールを使ってそれを読み込んで処理する、という流れを説明してきましたが、これらの処理を一括で行えるPythonモジュールがあったら便利ですね。
そこで気象ビジネス推進コンソーシアムでは独自にPythonモジュールwxbcgrib
を作成し、それらの機能をまとめて実装しました。
wxbcgrib
を使えば、上でやってきた処理を1行で完了させることができます。
それではまずwxbcgrib
のインポートから行います。
wxbcgrib
はこのipynbファイルと同じフォルダにwxbcgrib.py
を置くことでインポートすることができます。
このipynbファイルと同じフォルダにwxbcgrib.py
があることを確認して、以下のセルを実行して下さい。
import wxbcgrib as wg
それでは早速、GRIB2形式のMSM-GPVファイルを、直接wxbcgrib
を使って読み込んでみましょう。
読み込みには、wxbcgrib
のDS_from_grnc
関数を使います。
wg.DS_from_grnc(GRIB2ファイルへのパス, 取り出したい気象要素名)
それでは以下のセルを実行して下さい。
#--GRIB2データを読み込む
# 読み込む気象要素は気温とする
gpvpath = "weather/Z__C_RJTD_20200905120000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin"
var = wg.DS_from_grnc(gpvpath,"TMP_1D5maboveground")
これで気温の値を読み込むことができました。簡単ですね!
なお今回は気温の要素名がTMP_1D5maboveground
と知っていたので指定できますが、知らなかったらどうしましょうか?
その時は少々強引ですが、要素名を指定せずに実行すると、エラーメッセージとして有効な要素名を表示してくれます。
(wxbcgrib.py
製作者のドキュメントより)
wg.DS_from_grnc(gpvpath,"")
!! ラベルは以下から選択してください: TCDC_surface MCDC_surface HCDC_surface PRES_surface UGRD_10maboveground DSWRF_surface RH_1D5maboveground PRMSL_meansealevel APCP_surface TMP_1D5maboveground LCDC_surface VGRD_10maboveground
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-16-844688de3732> in <module> ----> 1 wg.DS_from_grnc(gpvpath,"") ~/Desktop/気象データ演習/wxbcgrib.py in DS_from_grnc(paths, label, lalomima) 163 for path in paths: 164 #print(path) --> 165 msgs = msgs + msg_from_a_grnc(path,label, lalomima) 166 msgs.sort(key = lambda x: (x.tini,x.fcst)) 167 DS = DataSet(msgs) ~/Desktop/気象データ演習/wxbcgrib.py in msg_from_a_grnc(path, label, lalomima) 151 for lbl in nclabels: 152 print(" "*4,lbl) --> 153 raise ValueError 154 nc.Dataset.close(ncds) 155 return msgs ValueError:
さて変数var
には、netCDF4.Dataset
クラスと同じで、オブジェクトが入っています。
まず要素名と単位を表示してみましょう。
# 要素名と単位
var.name, var.unit
('Temperature', 'K')
次に予想対象時刻(validtime)を見てみましょう。
# 予想対象時刻
var.time
array([datetime.datetime(2020, 9, 5, 21, 0), datetime.datetime(2020, 9, 5, 22, 0), datetime.datetime(2020, 9, 5, 23, 0), datetime.datetime(2020, 9, 6, 0, 0), datetime.datetime(2020, 9, 6, 1, 0), datetime.datetime(2020, 9, 6, 2, 0), datetime.datetime(2020, 9, 6, 3, 0), datetime.datetime(2020, 9, 6, 4, 0), datetime.datetime(2020, 9, 6, 5, 0), datetime.datetime(2020, 9, 6, 6, 0), datetime.datetime(2020, 9, 6, 7, 0), datetime.datetime(2020, 9, 6, 8, 0), datetime.datetime(2020, 9, 6, 9, 0), datetime.datetime(2020, 9, 6, 10, 0), datetime.datetime(2020, 9, 6, 11, 0), datetime.datetime(2020, 9, 6, 12, 0)], dtype=object)
validtime、はdatetime型という日時を表す特別なフォーマットで格納されています。
例えば、一番上はFT=0で、2020/09/05 21:00(JST)
であることを表します。
気温のデータそのものは、以下のようにしてアクセスできます。
# データ本体
print("データサイズ:", var.data.shape)
var.data
データサイズ: (16, 505, 481)
array([[[302.0578 , 302.03436, 302.01874, ..., 302.08905, 302.04218, 302.0031 ], [302.02655, 301.9875 , 301.9406 , ..., 302.0656 , 302.01874, 301.9875 ], [301.9953 , 301.9406 , 301.8781 , ..., 302.05 , 302.01093, 301.97968], ..., [283.0031 , 282.925 , 282.9328 , ..., 286.76874, 286.76874, 286.76093], [282.88593, 282.78436, 282.76093, ..., 286.8078 , 286.8 , 286.79218], [283.02655, 282.8703 , 282.8156 , ..., 286.83124, 286.82343, 286.82343]], [[301.89423, 301.80048, 301.6911 , ..., 302.04266, 302.0114 , 301.98798], [301.86298, 301.7458 , 301.6208 , ..., 302.03485, 301.9958 , 301.97235], [301.8239 , 301.6911 , 301.53485, ..., 302.01923, 301.98798, 301.95673], ..., [282.53485, 282.34735, 282.33954, ..., 286.83954, 286.8239 , 286.8239 ], [282.4411 , 282.2458 , 282.1989 , ..., 286.8708 , 286.86298, 286.86298], [282.6208 , 282.4411 , 282.3083 , ..., 286.89423, 286.8864 , 286.8864 ]], [[301.6412 , 301.46152, 301.28964, ..., 302.0162 , 301.98495, 301.9537 ], [301.59433, 301.39902, 301.1959 , ..., 302.00058, 301.96933, 301.9459 ], [301.55527, 301.33652, 301.11777, ..., 301.99277, 301.96152, 301.93027], ..., [282.27402, 281.97714, 281.90683, ..., 286.96152, 286.9459 , 286.93027], [282.33652, 281.99277, 281.85214, ..., 286.99277, 286.96933, 286.9537 ], [282.4459 , 282.2584 , 281.99277, ..., 287.00058, 286.97714, 286.97714]], ..., [[300.76288, 300.70038, 300.6457 , ..., 302.03632, 301.99725, 301.966 ], [300.72382, 300.6535 , 300.59882, ..., 302.04413, 302.00507, 301.97382], [300.67694, 300.60663, 300.55194, ..., 302.05194, 302.01288, 301.98163], ..., [287.5832 , 287.5832 , 287.69257, ..., 287.3332 , 287.36444, 287.41913], [287.4582 , 287.30194, 287.43475, ..., 287.41132, 287.42694, 287.466 ], [287.4035 , 287.2707 , 287.31757, ..., 287.48944, 287.48163, 287.50507]], [[300.9903 , 300.95905, 300.9356 , ..., 302.0528 , 302.01373, 301.98248], [300.96686, 300.9356 , 300.91217, ..., 302.0606 , 302.02155, 301.9903 ], [300.95905, 300.91998, 300.89655, ..., 302.06842, 302.02936, 301.9981 ], ..., [289.1153 , 289.20123, 289.41998, ..., 287.3106 , 287.34186, 287.38873], [288.88092, 288.7481 , 288.95123, ..., 287.38873, 287.40436, 287.44342], [288.77936, 288.63873, 288.7403 , ..., 287.47467, 287.46686, 287.4903 ]], [[301.36017, 301.3133 , 301.3133 , ..., 302.05548, 302.01642, 301.97736], [301.3758 , 301.33673, 301.33673, ..., 302.0633 , 302.01642, 301.98517], [301.3836 , 301.35236, 301.3758 , ..., 302.0711 , 302.02423, 301.98517], ..., [290.26642, 290.36798, 290.61798, ..., 287.29767, 287.32892, 287.3758 ], [290.01642, 289.89142, 290.11798, ..., 287.36017, 287.39142, 287.4383 ], [289.89142, 289.7586 , 289.8758 , ..., 287.42267, 287.4383 , 287.47736]]], dtype=float32)
データサイズは3次元になっており、(FT, Latitude, Longitude)
の数を表します。
緯度経度情報もあります。以下のセルを実行して下さい。
# 緯度
var.lat
# 経度を表示したい場合は以下をコメントアウトして実行
# var.lon
array([22.4 , 22.45, 22.5 , 22.55, 22.6 , 22.65, 22.7 , 22.75, 22.8 , 22.85, 22.9 , 22.95, 23. , 23.05, 23.1 , 23.15, 23.2 , 23.25, 23.3 , 23.35, 23.4 , 23.45, 23.5 , 23.55, 23.6 , 23.65, 23.7 , 23.75, 23.8 , 23.85, 23.9 , 23.95, 24. , 24.05, 24.1 , 24.15, 24.2 , 24.25, 24.3 , 24.35, 24.4 , 24.45, 24.5 , 24.55, 24.6 , 24.65, 24.7 , 24.75, 24.8 , 24.85, 24.9 , 24.95, 25. , 25.05, 25.1 , 25.15, 25.2 , 25.25, 25.3 , 25.35, 25.4 , 25.45, 25.5 , 25.55, 25.6 , 25.65, 25.7 , 25.75, 25.8 , 25.85, 25.9 , 25.95, 26. , 26.05, 26.1 , 26.15, 26.2 , 26.25, 26.3 , 26.35, 26.4 , 26.45, 26.5 , 26.55, 26.6 , 26.65, 26.7 , 26.75, 26.8 , 26.85, 26.9 , 26.95, 27. , 27.05, 27.1 , 27.15, 27.2 , 27.25, 27.3 , 27.35, 27.4 , 27.45, 27.5 , 27.55, 27.6 , 27.65, 27.7 , 27.75, 27.8 , 27.85, 27.9 , 27.95, 28. , 28.05, 28.1 , 28.15, 28.2 , 28.25, 28.3 , 28.35, 28.4 , 28.45, 28.5 , 28.55, 28.6 , 28.65, 28.7 , 28.75, 28.8 , 28.85, 28.9 , 28.95, 29. , 29.05, 29.1 , 29.15, 29.2 , 29.25, 29.3 , 29.35, 29.4 , 29.45, 29.5 , 29.55, 29.6 , 29.65, 29.7 , 29.75, 29.8 , 29.85, 29.9 , 29.95, 30. , 30.05, 30.1 , 30.15, 30.2 , 30.25, 30.3 , 30.35, 30.4 , 30.45, 30.5 , 30.55, 30.6 , 30.65, 30.7 , 30.75, 30.8 , 30.85, 30.9 , 30.95, 31. , 31.05, 31.1 , 31.15, 31.2 , 31.25, 31.3 , 31.35, 31.4 , 31.45, 31.5 , 31.55, 31.6 , 31.65, 31.7 , 31.75, 31.8 , 31.85, 31.9 , 31.95, 32. , 32.05, 32.1 , 32.15, 32.2 , 32.25, 32.3 , 32.35, 32.4 , 32.45, 32.5 , 32.55, 32.6 , 32.65, 32.7 , 32.75, 32.8 , 32.85, 32.9 , 32.95, 33. , 33.05, 33.1 , 33.15, 33.2 , 33.25, 33.3 , 33.35, 33.4 , 33.45, 33.5 , 33.55, 33.6 , 33.65, 33.7 , 33.75, 33.8 , 33.85, 33.9 , 33.95, 34. , 34.05, 34.1 , 34.15, 34.2 , 34.25, 34.3 , 34.35, 34.4 , 34.45, 34.5 , 34.55, 34.6 , 34.65, 34.7 , 34.75, 34.8 , 34.85, 34.9 , 34.95, 35. , 35.05, 35.1 , 35.15, 35.2 , 35.25, 35.3 , 35.35, 35.4 , 35.45, 35.5 , 35.55, 35.6 , 35.65, 35.7 , 35.75, 35.8 , 35.85, 35.9 , 35.95, 36. , 36.05, 36.1 , 36.15, 36.2 , 36.25, 36.3 , 36.35, 36.4 , 36.45, 36.5 , 36.55, 36.6 , 36.65, 36.7 , 36.75, 36.8 , 36.85, 36.9 , 36.95, 37. , 37.05, 37.1 , 37.15, 37.2 , 37.25, 37.3 , 37.35, 37.4 , 37.45, 37.5 , 37.55, 37.6 , 37.65, 37.7 , 37.75, 37.8 , 37.85, 37.9 , 37.95, 38. , 38.05, 38.1 , 38.15, 38.2 , 38.25, 38.3 , 38.35, 38.4 , 38.45, 38.5 , 38.55, 38.6 , 38.65, 38.7 , 38.75, 38.8 , 38.85, 38.9 , 38.95, 39. , 39.05, 39.1 , 39.15, 39.2 , 39.25, 39.3 , 39.35, 39.4 , 39.45, 39.5 , 39.55, 39.6 , 39.65, 39.7 , 39.75, 39.8 , 39.85, 39.9 , 39.95, 40. , 40.05, 40.1 , 40.15, 40.2 , 40.25, 40.3 , 40.35, 40.4 , 40.45, 40.5 , 40.55, 40.6 , 40.65, 40.7 , 40.75, 40.8 , 40.85, 40.9 , 40.95, 41. , 41.05, 41.1 , 41.15, 41.2 , 41.25, 41.3 , 41.35, 41.4 , 41.45, 41.5 , 41.55, 41.6 , 41.65, 41.7 , 41.75, 41.8 , 41.85, 41.9 , 41.95, 42. , 42.05, 42.1 , 42.15, 42.2 , 42.25, 42.3 , 42.35, 42.4 , 42.45, 42.5 , 42.55, 42.6 , 42.65, 42.7 , 42.75, 42.8 , 42.85, 42.9 , 42.95, 43. , 43.05, 43.1 , 43.15, 43.2 , 43.25, 43.3 , 43.35, 43.4 , 43.45, 43.5 , 43.55, 43.6 , 43.65, 43.7 , 43.75, 43.8 , 43.85, 43.9 , 43.95, 44. , 44.05, 44.1 , 44.15, 44.2 , 44.25, 44.3 , 44.35, 44.4 , 44.45, 44.5 , 44.55, 44.6 , 44.65, 44.7 , 44.75, 44.8 , 44.85, 44.9 , 44.95, 45. , 45.05, 45.1 , 45.15, 45.2 , 45.25, 45.3 , 45.35, 45.4 , 45.45, 45.5 , 45.55, 45.6 , 45.65, 45.7 , 45.75, 45.8 , 45.85, 45.9 , 45.95, 46. , 46.05, 46.1 , 46.15, 46.2 , 46.25, 46.3 , 46.35, 46.4 , 46.45, 46.5 , 46.55, 46.6 , 46.65, 46.7 , 46.75, 46.8 , 46.85, 46.9 , 46.95, 47. , 47.05, 47.1 , 47.15, 47.2 , 47.25, 47.3 , 47.35, 47.4 , 47.45, 47.5 , 47.55, 47.6 ])
wxbcgribでは、緯度経度を指定して時系列データを取り出すことができます。
tsメソッドを使って、たった1行のコードでできてしまいます。引数として緯度・経度の順番で与えます。
変数名.ts(緯度, 経度)
今回は北緯35度・東経140度の値を取り出してみましょう。以下のセルを実行して下さい。
# 北緯35度・東経140度の時系列データを取得
ts_temp = var.ts(35.0, 140.0)
ts_temp
array([300.13593, 300.09735, 300.31308, 300.9329 , 301.23544, 301.71756, 301.7408 , 301.71555, 301.74518, 301.79312, 301.6907 , 301.31134, 300.7464 , 300.56757, 301.1153 , 301.3836 ], dtype=float32)
こんなに簡単にできてしまうなんて、便利ですね!
これで、第1章・第2章とで行った様々な処理を、wxbcgribで簡単に実行することができました。
それでは次は、北緯26.2度・東経125.75度の地上気圧の時系列データを、wxbcgrib
で取得してみましょう。
### 練習問題(4)
# 北緯26.2度・東経125.75度の地上気圧の時系列データを、wxbcgribで取得する
var = wg.DS_from_grnc(gpvpath,"PRES_surface")
ts_pres = var.ts(26.2, 125.75)
ts_pres
array([ 99977.96 , 100037.16 , 99968.336, 99906.9 , 99836.15 , 99770.95 , 99718.2 , 99668.4 , 99639.945, 99644.78 , 99638.1 , 99659.12 , 99654.36 , 99669.09 , 99662.49 , 99648.06 ], dtype=float32)
以上で、wxbcgribを使ってGRIB2形式のGPVデータをハンドリングすることができました。
目次へ戻る
第1章〜第3章でGPVデータをPythonで読み込むことができました。
すると次に必要になる処理として、読み込んだ気象要素の単位変換や、別の気象要素への変換があります。
これらの計算を簡単に行う方法として、Weather Data Science 加藤芳樹・史葉のオリジナルPythonモジュールwxparams
を紹介します。
wxparamsは、気象要素の変換計算や単位変換の計算を行う関数群から成り立っています。
ではwxparams
と、解説で使うnumpy
・pandas
をインポートします。以下のセルを実行して下さい。
import wxparams as wx
import numpy as np
import pandas as pd
まずはwxparams
を簡単に使ってみます。
気温20度・相対湿度75%の時の露点温度を計算してみましょう。
wxparams
はインプットとしてnumpy.ndarray
型を想定していますので、Numpy配列の形にして引数を渡します。
また露点温度の計算にはRH_to_Td
関数を使います。
以下のセルを実行して下さい。
# numpy.ndarray型でインプット
TEMP = np.array([20.0])
RH = np.array([75.0])
Td = wx.RH_to_Td(TEMP, RH)
Td
array([15.4380164])
戻り値もnumpy.ndarray
型になっています。
またpandas.Series
型でもOKです。まず以下のセルを実行して、データフレームを作ります。
# pandas.DataFrameを作成
df = pd.DataFrame({"TEMP":TEMP, "RH":RH})
df
TEMP | RH | |
---|---|---|
0 | 20.0 | 75.0 |
引数としてTEMP列とRH列を、以下のようにRH_to_Td
関数へ渡すこともできます。
df["Td"] = wx.RH_to_Td(df["TEMP"], df["RH"])
df
TEMP | RH | Td | |
---|---|---|---|
0 | 20.0 | 75.0 | 15.438016 |
簡単に概要を説明しましたので、練習用データフレームでwxparams
を使ってみましょう。
まず以下のセルで練習用データを定義します。
# 気温[K], 相対湿度[%], 気圧[Pa]
df = pd.DataFrame({"気温":[270.65, 292.45, 303.95],
"相対湿度":[85.2, 98.3, 48.1],
"気圧":[50000, 92500, 101325],
})
df
気温 | 相対湿度 | 気圧 | |
---|---|---|---|
0 | 270.65 | 85.2 | 50000 |
1 | 292.45 | 98.3 | 92500 |
2 | 303.95 | 48.1 | 101325 |
データの内容は、気温[K], 相対湿度[%], 気圧[Pa]です。
では露点温度を計算してみましょう。注意点として、RH_to_Td
関数へは気温を摂氏温度[C]で渡す必要があります。
### 練習問題(5)
# 露点温度[C]を計算する
df["気温"] = df["気温"] - 273.15
df["露点温度"] = wx.RH_to_Td(df["気温"], df["相対湿度"])
df
気温 | 相対湿度 | 気圧 | 露点温度 | |
---|---|---|---|---|
0 | -2.5 | 85.2 | 50000 | -4.642878 |
1 | 19.3 | 98.3 | 92500 | 19.025067 |
2 | 30.8 | 48.1 | 101325 | 18.571997 |
では次にTheta_e
関数で相当温位を計算してみましょう。各引数の単位に注意して下さい。
### 練習問題(6)
# 相当温位[K]を計算する
df["相当温位"] = wx.Theta_e(df["気温"], df["露点温度"], df["気圧"]/100)
df
気温 | 相対湿度 | 気圧 | 露点温度 | 相当温位 | |
---|---|---|---|---|---|
0 | -2.5 | 85.2 | 50000 | -4.642878 | 348.423863 |
1 | 19.3 | 98.3 | 92500 | 19.025067 | 343.375966 |
2 | 30.8 | 48.1 | 101325 | 18.571997 | 342.789851 |
ここまでは、テーブルデータ(pandas.DataFrame)の形で気象データを持っている場合を練習してきました。
次は、GPVから取り出したデータを使い、GPVには含まれていない気象要素を計算してみましょう。
ここでは計算する気象要素として『風向・風速』を扱ってみます。
第3章で使ったGPVデータをwxbcgrib
を使って読み込み、風のUV成分(東西風・南北風)を入手しましょう。
### 練習問題(7)
#--モデルデータから風向・風速を求める
# まずは風のUV成分を入手
U_wind = wg.DS_from_grnc(gpvpath,"UGRD_10maboveground")
V_wind = wg.DS_from_grnc(gpvpath,"VGRD_10maboveground")
読み込んだデータのデータサイズ(shape)と、単位を確認しましょう。
# データサイズをprint
print( U_wind.data.shape )
# 気象要素の単位をprint
print( U_wind.unit )
(16, 505, 481) m/s
読み込んだデータは3次元のnumpy.ndarray型ですが、wxparams
は多次元のNumpy配列を想定しています。
よって、多次元配列のまま引数として与えることができます。
それではUV_to_SpdDir
関数を用いて風向・風速を求めてみましょう。
注意点として、戻り値は2つあり、風速・風向の順番でreturnされてきます。
また、風速のデータサイズ(shape)と、FT=0における値を出力してみましょう。
### 練習問題(8)
# 風向・風速を計算
Wspd, Wdir = wx.UV_to_SpdDir(U_wind.data, V_wind.data)
# 風速のデータサイズを出力
print( Wspd.shape )
# FT=0の値を出力
Wspd[0,:,:]
(16, 505, 481)
array([[9.358936 , 9.180615 , 8.881036 , ..., 6.524155 , 6.6033683, 6.7736897], [9.375007 , 9.285649 , 9.040425 , ..., 6.5907264, 6.6962523, 6.863365 ], [9.420724 , 9.390047 , 9.168822 , ..., 6.6833477, 6.847068 , 6.9843736], ..., [1.5186125, 1.4945667, 1.4728137, ..., 5.7653747, 5.9132986, 6.1082764], [1.7121074, 1.7593732, 1.6912637, ..., 5.7870584, 5.94015 , 6.135818 ], [1.8670738, 1.9227737, 1.9341333, ..., 5.8708215, 5.996276 , 6.107845 ]], dtype=float32)
以上でwxparams
を使って気象要素の変換計算と、単位変換の計算をすることができました。
目次へ戻る