GrADS　MSM気圧面データ　取得＆描画　バイナリファイル編

SHARE

MSMデータを用いて、GrADSで描画したいと思います。

MSMデータとは、なんぞやという方は、こちらを参照してください。

GrADS　MSMデータを取得＆描画　NetCDF編

MSMデータの取得

本記事では、京都大学の気象庁データを利用させていただきます。

オリジナルのMSMデータは、国際気象通報式FM92 GRIB 二進形式格子点資料気象通報式(第2版)というデータ形式になっています。

本記事では、そのオリジナルデータ(GRIB2形式)から描画させたいと思います。

Z__C_RJTD_yyyyMMddhhmmss_MSM_GPV_Rjp_Lsurf_FHhh-hh_grib2.binは、地表面データです。

Z__C_RJTD_yyyyMMddhhmmss_MSM_GPV_Rjp_L-pall_FHhh-hh_grib2.binは、気圧面データです。

本記事では、気圧面データに着目します。

地表面データはこちら

GrADS　MSM地表面データ　取得＆描画　バイナリファイル編

yyyyは西暦、MMは月、ddは日、hhは時間、㎜は分、ssは秒です。

wgetコマンドで取得しましょう。

$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin

上では、2021年7月17日の00UTC初期値の気圧面データであるMSMデータを取得しています。予測時間は、15時間後までです。

51時間先までほしいので、合計で4つのファイルを取得します。また、初期値が03,06,09,15,18,21時の場合、39時間先までですので、合計のファイルは、３つです。

$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH18-33_grib2.bin
$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH36-39_grib2.bin
$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH42-51_grib2.bin

wgrib2コマンドでバイナリファイルに変換します。

$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin -no_header -bin out0015.bin
$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH18-33_grib2.bin -no_header -bin out1833.bin
$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH36-39_grib2.bin -no_header -bin out3639.bin
$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_L-pall_FH42-51_grib2.bin -no_header -bin out4251.bin

バイナリファイルの編集

この時のバイナリファイルout????.binは、1000hPaのジオポテンシャル高度、東西風、南北風、気温、鉛直流、相対湿度、975hPaのジオポテンシャル高度、東西風、南北風、….100hPaの相対湿度という並び方をしています。

このままだと、GrADSで描画する際、ある高度のジオポテンシャル高度を見たいのに気温が表示されたりとめちゃくちゃになるので、それぞれの変数でまとめるという感じに変更する必要があります。

また、51時間先まで４つのファイルに分かれているので、1つのファイルに統合しようと思います。

Fortranでやります。Pythonもいつかやろうと思います。

バイナリファイルの大きさを確認します。

$ ls -l
-rw-r--r-- 1 hoge users 134628384 Jul 27 20:00 out0015.bin
-rw-r--r-- 1 hoge users 134628384 Jul 27 20:00 out1833.bin
-rw-r--r-- 1 hoge users  44876128 Jul 27 20:00 out3639.bin
-rw-r--r-- 1 hoge users  89752256 Jul 27 20:00 out4251.bin

out0015.binとout1833.binは、ともに134628384byte、out3639.binは、44876128byte、out4251.binは、89752256byteとなっています。

本当にファイルサイズと想定しているファイル形式と一致しているのか計算します。

気象業務センターで確認すると、MSMデータは、北緯22.4~47.6度と東経120~150度で、気圧面データでは、x方向(経度)格子点間隔0.125の格子点数241、y方向(緯度)格子間隔0.1度の格子点数253です。

z方向（1000hPa・975hPa・950hPa・925hPa・900hPa・850hPa・800hPa・700hPa・600hPa・500hPa・400hPa・300hPa・250hPa・200hPa・150hPa・100hPa）の16個 あります。

物理量として、 ジオポテンシャル高度、東西風、南北風、気温、鉛直流、相対湿度(~300hPaまで)があります。

out0015.binとout1833.binは、時間間隔3時間で18時間先までのファイルですので、t方向に６です。

1データが4byte、物理量は、6つありますが、相対湿度が300hPaまでしかありません。

4×5×241×253×16×6＝117,068,160

4×1×241×253×12×6＝17,560,224

117,068,160+17560,224=134,628,384となり一致しました。

同様にout3639.binは、時間間隔3時間で6時先までのファイルですのでt方向に2です。

4×5×241×253×16×2＝39,022,720

4×1×241×253×12×2＝5,853,408

39,022,720+ 5,853,408 =44,876,128となり一致しました。

同様にout4251.binは、時間間隔3時間で12時間先までのファイルですのでt方向に4です。

4×5×241×253×16×4＝78,045,440

4×1×241×253×12×4＝11,706,816

78,045,440+11,706,816=89,752,256となり一致しました。

想定したファイルサイズになったので、プログラムを書きます。

implicit none
real,allocatable::out0015(:),out1833(:),out3639(:),out4251(:)
real,allocatable::gz0015(:,:,:,:),uu0015(:,:,:,:),vv0015(:,:,:,:),temp0015(:,:,:,:),ww0015(:,:,:,:),rh0015(:,:,:,:)
real,allocatable::gz1833(:,:,:,:),uu1833(:,:,:,:),vv1833(:,:,:,:),temp1833(:,:,:,:),ww1833(:,:,:,:),rh1833(:,:,:,:)
real,allocatable::gz3639(:,:,:,:),uu3639(:,:,:,:),vv3639(:,:,:,:),temp3639(:,:,:,:),ww3639(:,:,:,:),rh3639(:,:,:,:)
real,allocatable::gz4251(:,:,:,:),uu4251(:,:,:,:),vv4251(:,:,:,:),temp4251(:,:,:,:),ww4251(:,:,:,:),rh4251(:,:,:,:)
real,allocatable::gz(:,:,:,:),uu(:,:,:,:),vv(:,:,:,:),temp(:,:,:,:),ww(:,:,:,:),rh(:,:,:,:)
real,allocatable::data(:,:,:,:,:)
integer::xx,yy,zz,tt,num,tdata,rnum
integer,parameter::xxdim=241,yydim=253,zzdim=16,humdim=12
integer,parameter::out0015ttdim=6,out1833ttdim=6,out3639ttdim=2,out4251ttdim=4,ttdim=16
integer::out0015recl,out1833recl,out3639recl,out4251recl
out0015recl = 4 * 5 * xxdim * yydim * zzdim * out0015ttdim + 4 * 1 * xxdim * yydim * zzdim * out0015ttdim
out1833recl = 4 * 5 * xxdim * yydim * zzdim * out1833ttdim + 4 * 1 * xxdim * yydim * zzdim * out1833ttdim
out3639recl = 4 * 5 * xxdim * yydim * zzdim * out3639ttdim + 4 * 1 * xxdim * yydim * zzdim * out3639ttdim
out4251recl = 4 * 5 * xxdim * yydim * zzdim * out4251ttdim + 4 * 1 * xxdim * yydim * zzdim * out4251ttdim
allocate(gz0015(xxdim,yydim,zzdim,out0015ttdim),uu0015(xxdim,yydim,zzdim,out0015ttdim),vv0015(xxdim,yydim,zzdim,out0015ttdim),&
temp0015(xxdim,yydim,zzdim,out0015ttdim),ww0015(xxdim,yydim,zzdim,out0015ttdim),rh0015(xxdim,yydim,zzdim,out0015ttdim))
allocate(gz1833(xxdim,yydim,zzdim,out1833ttdim),uu1833(xxdim,yydim,zzdim,out1833ttdim),vv1833(xxdim,yydim,zzdim,out1833ttdim),&
temp1833(xxdim,yydim,zzdim,out1833ttdim),ww1833(xxdim,yydim,zzdim,out1833ttdim),rh1833(xxdim,yydim,zzdim,out1833ttdim))
allocate(gz3639(xxdim,yydim,zzdim,out3639ttdim),uu3639(xxdim,yydim,zzdim,out3639ttdim),vv3639(xxdim,yydim,zzdim,out3639ttdim),&
temp3639(xxdim,yydim,zzdim,out3639ttdim),ww3639(xxdim,yydim,zzdim,out3639ttdim),rh3639(xxdim,yydim,zzdim,out3639ttdim))
allocate(gz4251(xxdim,yydim,zzdim,out4251ttdim),uu4251(xxdim,yydim,zzdim,out4251ttdim),vv4251(xxdim,yydim,zzdim,out4251ttdim),&
temp4251(xxdim,yydim,zzdim,out4251ttdim),ww4251(xxdim,yydim,zzdim,out4251ttdim),rh4251(xxdim,yydim,zzdim,out4251ttdim))
allocate(gz(xxdim,yydim,zzdim,ttdim),uu(xxdim,yydim,zzdim,ttdim),vv(xxdim,yydim,zzdim,ttdim),temp(xxdim,yydim,zzdim,ttdim),&
ww(xxdim,yydim,zzdim,ttdim),rh(xxdim,yydim,zzdim,ttdim))
allocate(data(xxdim,yydim,zzdim,ttdim,6))
allocate(out0015(out0015recl/4),out1833(out1833recl/4),out3639(out3639recl/4),out4251(out4251recl/4))
open(10,file="out0015.bin",access="direct",recl=out0015recl)
open(11,file="out1833.bin",access="direct",recl=out1833recl)
open(12,file="out3639.bin",access="direct",recl=out3639recl)
open(13,file="out4251.bin",access="direct",recl=out4251recl)
read(10,rec=1)out0015(:)
read(11,rec=1)out1833(:)
read(12,rec=1)out3639(:)
read(13,rec=1)out4251(:)
close(10)
close(11)
close(12)
close(13)
!read
num=1
do tt=1,out0015ttdim
do zz=1,humdim
call sub0015(gz0015)
call sub0015(uu0015)
call sub0015(vv0015)
call sub0015(temp0015)
call sub0015(ww0015)
call sub0015(rh0015)
enddo
do zz=humdim+1,zzdim
call sub0015(gz0015)
call sub0015(uu0015)
call sub0015(vv0015)
call sub0015(temp0015)
call sub0015(ww0015)
do yy=1,yydim
do xx=1,xxdim
rh0015(xx,yy,zz,tt)=0.0
enddo
enddo
enddo
enddo
num=1
do tt=1,out1833ttdim
do zz=1,humdim
call sub1833(gz1833)
call sub1833(uu1833)
call sub1833(vv1833)
call sub1833(temp1833)
call sub1833(ww1833)
call sub1833(rh1833)
enddo
do zz=humdim+1,zzdim
call sub1833(gz1833)
call sub1833(uu1833)
call sub1833(vv1833)
call sub1833(temp1833)
call sub1833(ww1833)
do yy=1,yydim
do xx=1,xxdim
rh1833(xx,yy,zz,tt)=0.0
enddo
enddo
enddo
enddo

num=1
do tt=1,out3639ttdim
do zz=1,humdim
call sub3639(gz3639)
call sub3639(uu3639)
call sub3639(vv3639)
call sub3639(temp3639)
call sub3639(ww3639)
call sub3639(rh3639)
enddo
do zz=humdim+1,zzdim
call sub3639(gz3639)
call sub3639(uu3639)
call sub3639(vv3639)
call sub3639(temp3639)
call sub3639(ww3639)
do yy=1,yydim
do xx=1,xxdim
rh3639(xx,yy,zz,tt)=0.0
enddo
enddo
enddo
enddo
num=1
do tt=1,out4251ttdim
do zz=1,humdim
call sub4251(gz4251)
call sub4251(uu4251)
call sub4251(vv4251)
call sub4251(temp4251)
call sub4251(ww4251)
call sub4251(rh4251)
enddo
do zz=humdim+1,zzdim
call sub4251(gz4251)
call sub4251(uu4251)
call sub4251(vv4251)
call sub4251(temp4251)
call sub4251(ww4251)
do yy=1,yydim
do xx=1,xxdim
rh4251(xx,yy,zz,tt)=0.0
enddo
enddo
enddo
enddo

!write

do tdata=1,16
do zz=1,zzdim
do yy=1,yydim
do xx=1,xxdim
if (tdata <= out0015ttdim)then
gz(xx,yy,zz,tdata)=gz0015(xx,yy,zz,tdata)
uu(xx,yy,zz,tdata)=uu0015(xx,yy,zz,tdata)
vv(xx,yy,zz,tdata)=vv0015(xx,yy,zz,tdata)
temp(xx,yy,zz,tdata)=temp0015(xx,yy,zz,tdata)
ww(xx,yy,zz,tdata)=ww0015(xx,yy,zz,tdata)
rh(xx,yy,zz,tdata)=rh0015(xx,yy,zz,tdata)

elseif(tdata <= out0015ttdim+out1833ttdim)then
tt=tdata - out0015ttdim
gz(xx,yy,zz,tdata)=gz1833(xx,yy,zz,tt)
uu(xx,yy,zz,tdata)=uu1833(xx,yy,zz,tt)
vv(xx,yy,zz,tdata)=vv1833(xx,yy,zz,tt)
temp(xx,yy,zz,tdata)=temp1833(xx,yy,zz,tt)
ww(xx,yy,zz,tdata)=ww1833(xx,yy,zz,tt)
rh(xx,yy,zz,tdata)=rh1833(xx,yy,zz,tt)
elseif(tdata <= out0015ttdim+out1833ttdim+out3639ttdim)then
tt=tdata - out0015ttdim - out1833ttdim
gz(xx,yy,zz,tdata)=gz3639(xx,yy,zz,tt)
uu(xx,yy,zz,tdata)=uu3639(xx,yy,zz,tt)
vv(xx,yy,zz,tdata)=vv3639(xx,yy,zz,tt)
temp(xx,yy,zz,tdata)=temp3639(xx,yy,zz,tt)
ww(xx,yy,zz,tdata)=ww3639(xx,yy,zz,tt)
rh(xx,yy,zz,tdata)=rh3639(xx,yy,zz,tt)
else
tt=tdata - out0015ttdim - out1833ttdim - out3639ttdim
gz(xx,yy,zz,tdata)=gz4251(xx,yy,zz,tt)
uu(xx,yy,zz,tdata)=uu4251(xx,yy,zz,tt)
vv(xx,yy,zz,tdata)=vv4251(xx,yy,zz,tt)
temp(xx,yy,zz,tdata)=temp4251(xx,yy,zz,tt)
ww(xx,yy,zz,tdata)=ww4251(xx,yy,zz,tt)
rh(xx,yy,zz,tdata)=rh4251(xx,yy,zz,tt)
endif
enddo
enddo
enddo
enddo
!write(*,*)ww(:,:,:,:)
open(20,file="msm.bin",status="unknown",access="direct",form="unformatted",recl=4*xxdim*yydim*zzdim)
rnum=1
do tt=1,ttdim
write(20,rec=rnum)gz(:,:,:,tt)
rnum=rnum+1
write(20,rec=rnum)uu(:,:,:,tt)
rnum=rnum+1
write(20,rec=rnum)vv(:,:,:,tt)
rnum=rnum+1
write(20,rec=rnum)temp(:,:,:,tt)
rnum=rnum+1
write(20,rec=rnum)ww(:,:,:,tt)
rnum=rnum+1
write(20,rec=rnum)rh(:,:,:,tt)
rnum=rnum+1
write(*,*)ww(:,:,:,tt)
enddo
close(20)

stop

contains

subroutine sub0015(var)
implicit none
real::var(xxdim,yydim,zzdim,out0015ttdim),fread(out0015recl/4)
do yy=1,yydim
do xx=1,xxdim
var(xx,yy,zz,tt)=out0015(num)
num=num+1
enddo
enddo
end subroutine
subroutine sub1833(var)
implicit none
real::var(xxdim,yydim,zzdim,out1833ttdim),fread(out1833recl/4)
do yy=1,yydim
do xx=1,xxdim
var(xx,yy,zz,tt)=out1833(num)
num=num+1
enddo
enddo
end subroutine
subroutine sub3639(var)
implicit none
real::var(xxdim,yydim,zzdim,out3639ttdim),fread(out3639recl/4)
do yy=1,yydim
do xx=1,xxdim
var(xx,yy,zz,tt)=out3639(num)
num=num+1
enddo
enddo
end subroutine
subroutine sub4251(var)
implicit none
real::var(xxdim,yydim,zzdim,out4251ttdim),fread(out4251recl/4)
do yy=1,yydim
do xx=1,xxdim
var(xx,yy,zz,tt)=out4251(num)
num=num+1
enddo
enddo
end subroutine

end

コントロールファイルの作成

GrADSで描画するために、コントロールファイルを作る必要があります。

気象業務センターで確認すると、MSMデータは、北緯22.4~47.6度と東経120~150度で、気圧面データでは、x方向(経度)格子点間隔0.125の格子点数241、y方向(緯度)格子間隔0.1度の格子点数253です。

また、気圧面データは、z方向（1000hPa・975hPa・950hPa・925hPa・900hPa・850hPa・800hPa・700hPa・600hPa・500hPa・400hPa・300hPa・250hPa・200hPa・150hPa・100hPa）の16個あります。

この情報の下、コントロールファイルを作成します。

気圧面データのコントロールファイル

DSET msm.bin
UNDEF 999999999E+20
XDEF 241 LINEAR 120 0.125
YDEF 253 LINEAR 22.4 0.1
ZDEF 16 LEVELS 1000 975 950 925 900 850 800 700 600 500 400 300 250 200 150 100
TDEF 16 LINEAR 00z01jan0000 3hr
VARS 6
gz 16 99 gz
u 16 99 u
v 16 99 v
temp 16 99 temp
w 16 99 w
rh 16 99 rh
ENDVARS

GrADSで描画

本記事では、GrADSで描画します。

openでコントロールファイルを開きます。

open msm_p.ctl

例として、以下のように実行することで気温の図(1000hPa)が描画されます。

ga ->d temp

GrADSの詳細はこちらで

前の記事

compositeコマンド　透過処理

次の記事

GrADS　MSM地表面データ　取得＆描画　バイナリファイル編