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_Lsurf_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_Lsurf_FH16-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_Lsurf_FH34-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_Lsurf_FH40-51_grib2.bin

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

$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin -no_header -bin out0015.bin
$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_Lsurf_FH16-33_grib2.bin -no_header -bin out1633.bin
$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_Lsurf_FH34-39_grib2.bin -no_header -bin out3439.bin
$ wgrib2 Z__C_RJTD_20210717000000_MSM_GPV_Rjp_Lsurf_FH40-51_grib2.bin -no_header -bin out4051.bin

バイナリファイルの編集

この時のバイナリファイルout????.binは、海面更生気圧、地上気圧、東西風、南北風、気温、相対湿度、1時間降水量、上層雲量、中層雲量、下層雲量、全雲量、日射量、１時先の海面更生気圧、…..51時間先の日射量という並び方をしています。

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

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

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

$ ls -l
-rw-r--r-- 1 hoge users 69241393 Jul 27 20:00 out0015.bin
-rw-r--r-- 1 hoge users 78716561 Jul 27 20:00 out1633.bin
-rw-r--r-- 1 hoge users 26238929 Jul 27 20:00 out3439.bin
-rw-r--r-- 1 hoge users 52477745 Jul 27 20:00 out4051.bin

out0015.binは、184607800byte、out1633.binは、209869920byte、out3439.binは、69956640byte、out4051.binは、139913280byteとなっています。

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

気象業務センターで確認すると、MSMデータは、北緯22.4~47.6度と東経120~150度で、地表面データでは、x方向(経度)格子点間隔0.0625の格子点数481、y方向(緯度)格子間隔0.05度の格子点数505です。

物理量として、 海面更生気圧、地上気圧、東西風、南北風、気温、相対湿度、1時間降水量、上層雲量、中層雲量、下層雲量、全雲量、日射量 があります。

0015.binは、時間間隔1時間で16時間先までのファイルですので、t方向に15です。

1データが4byte、物理量は12個あります。しかし、初期場には、1時間降水量及び日射量がないため、結果10個の物理量がt方向に16あり、2個の物理量がt方向に15あります。

4×10×481×505×16=155459200

4×2×481×505×15=29148600

155459200+29148600=184607800

out1633.binでは、t方向に18です。

4×12×481×505×18=209869920

out3439.binでは、t方向に6です。

4×12×481×505×6=69956640

out4051.binでは、t方向に12です。

4×12×481×505×12=139913280

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

implicit none
real,allocatable::out0015(:),out1633(:),out3439(:),out4051(:)
real,allocatable::psea0015(:,:,:),sp0015(:,:,:),uu0015(:,:,:),vv0015(:,:,:),temp0015(:,:,:),rh0015(:,:,:),r1h0015(:,:,:)
real,allocatable::ncld_upper0015(:,:,:),ncld_mid0015(:,:,:),ncld_low0015(:,:,:),clda0015(:,:,:),dswrf0015(:,:,:)
real,allocatable::psea1633(:,:,:),sp1633(:,:,:),uu1633(:,:,:),vv1633(:,:,:),temp1633(:,:,:),rh1633(:,:,:),r1h1633(:,:,:)
real,allocatable::ncld_upper1633(:,:,:),ncld_mid1633(:,:,:),ncld_low1633(:,:,:),clda1633(:,:,:),dswrf1633(:,:,:)
real,allocatable::psea3439(:,:,:),sp3439(:,:,:),uu3439(:,:,:),vv3439(:,:,:),temp3439(:,:,:),rh3439(:,:,:),r1h3439(:,:,:)
real,allocatable::ncld_upper3439(:,:,:),ncld_mid3439(:,:,:),ncld_low3439(:,:,:),clda3439(:,:,:),dswrf3439(:,:,:)
real,allocatable::psea4051(:,:,:),sp4051(:,:,:),uu4051(:,:,:),vv4051(:,:,:),temp4051(:,:,:),rh4051(:,:,:),r1h4051(:,:,:)
real,allocatable::ncld_upper4051(:,:,:),ncld_mid4051(:,:,:),ncld_low4051(:,:,:),clda4051(:,:,:),dswrf4051(:,:,:)
real,allocatable::psea(:,:,:),sp(:,:,:),uu(:,:,:),vv(:,:,:),temp(:,:,:),rh(:,:,:),r1h(:,:,:)
real,allocatable::ncld_upper(:,:,:),ncld_mid(:,:,:),ncld_low(:,:,:),clda(:,:,:),dswrf(:,:,:)
real,allocatable::data(:,:,:)
integer::xx,yy,zz,tt,num,tdata,rnum
integer,parameter::xxdim=481,yydim=505
integer,parameter::out0015ttdim=16,out1633ttdim=18,out3439ttdim=6,out4051ttdim=12,ttdim=52
integer::out0015recl,out1633recl,out3439recl,out4051recl
out0015recl = 4 * 10 * xxdim * yydim * out0015ttdim + 4 * 2 * xxdim * yydim * (out0015ttdim -1)
out1633recl = 4 * 12 * xxdim * yydim * out1633ttdim 
out3439recl = 4 * 12 * xxdim * yydim * out3439ttdim 
out4051recl = 4 * 12 * xxdim * yydim * out4051ttdim
allocate(psea0015(xxdim,yydim,out0015ttdim),sp0015(xxdim,yydim,out0015ttdim),uu0015(xxdim,yydim,out0015ttdim),&
vv0015(xxdim,yydim,out0015ttdim),temp0015(xxdim,yydim,out0015ttdim),rh0015(xxdim,yydim,out0015ttdim),&
r1h0015(xxdim,yydim,out0015ttdim),ncld_upper0015(xxdim,yydim,out0015ttdim),&
ncld_mid0015(xxdim,yydim,out0015ttdim),ncld_low0015(xxdim,yydim,out0015ttdim),&
clda0015(xxdim,yydim,out0015ttdim),dswrf0015(xxdim,yydim,out0015ttdim))

allocate(psea1633(xxdim,yydim,out1633ttdim),sp1633(xxdim,yydim,out1633ttdim),uu1633(xxdim,yydim,out1633ttdim),&
vv1633(xxdim,yydim,out1633ttdim),temp1633(xxdim,yydim,out1633ttdim),rh1633(xxdim,yydim,out1633ttdim),&
r1h1633(xxdim,yydim,out1633ttdim),ncld_upper1633(xxdim,yydim,out1633ttdim),&
ncld_mid1633(xxdim,yydim,out1633ttdim),ncld_low1633(xxdim,yydim,out1633ttdim),&
clda1633(xxdim,yydim,out1633ttdim),dswrf1633(xxdim,yydim,out1633ttdim))

allocate(psea3439(xxdim,yydim,out3439ttdim),sp3439(xxdim,yydim,out3439ttdim),uu3439(xxdim,yydim,out3439ttdim),&
vv3439(xxdim,yydim,out3439ttdim),temp3439(xxdim,yydim,out3439ttdim),rh3439(xxdim,yydim,out3439ttdim),&
r1h3439(xxdim,yydim,out3439ttdim),ncld_upper3439(xxdim,yydim,out3439ttdim),&
ncld_mid3439(xxdim,yydim,out3439ttdim),ncld_low3439(xxdim,yydim,out3439ttdim),&
clda3439(xxdim,yydim,out3439ttdim),dswrf3439(xxdim,yydim,out3439ttdim))

allocate(psea4051(xxdim,yydim,out4051ttdim),sp4051(xxdim,yydim,out4051ttdim),uu4051(xxdim,yydim,out4051ttdim),&
vv4051(xxdim,yydim,out4051ttdim),temp4051(xxdim,yydim,out4051ttdim),rh4051(xxdim,yydim,out4051ttdim),&
r1h4051(xxdim,yydim,out4051ttdim),ncld_upper4051(xxdim,yydim,out4051ttdim),&
ncld_mid4051(xxdim,yydim,out4051ttdim),ncld_low4051(xxdim,yydim,out4051ttdim),&
clda4051(xxdim,yydim,out4051ttdim),dswrf4051(xxdim,yydim,out4051ttdim))

allocate(psea(xxdim,yydim,ttdim),sp(xxdim,yydim,ttdim),uu(xxdim,yydim,ttdim),vv(xxdim,yydim,ttdim),&
temp(xxdim,yydim,ttdim),rh(xxdim,yydim,ttdim),r1h(xxdim,yydim,ttdim),ncld_upper(xxdim,yydim,ttdim),&
ncld_mid(xxdim,yydim,ttdim),ncld_low(xxdim,yydim,ttdim),clda(xxdim,yydim,ttdim),dswrf(xxdim,yydim,ttdim))
allocate(data(xxdim,yydim,ttdim))
allocate(out0015(out0015recl/4),out1633(out1633recl/4),out3439(out3439recl/4),out4051(out4051recl/4))
open(10,file="out0015.bin",access="direct",recl=out0015recl)
open(11,file="out1633.bin",access="direct",recl=out1633recl)
open(12,file="out3439.bin",access="direct",recl=out3439recl)
open(13,file="out4051.bin",access="direct",recl=out4051recl)
read(10,rec=1)out0015(:)
read(11,rec=1)out1633(:)
read(12,rec=1)out3439(:)
read(13,rec=1)out4051(:)
close(10)
close(11)
close(12)
close(13)
!read
num=1
do tt=1,out0015ttdim
   if(tt==1)then
      call sub0015(psea0015)
      call sub0015(sp0015)
      call sub0015(uu0015)
      call sub0015(vv0015)
      call sub0015(temp0015)
      call sub0015(rh0015)
      do yy=1,yydim
      do xx=1,xxdim
         r1h0015(xx,yy,tt)=0.0
      enddo
      enddo
      call sub0015(ncld_upper0015)
      call sub0015(ncld_mid0015)
      call sub0015(ncld_low0015)
      call sub0015(clda0015)
      do yy=1,yydim
      do xx=1,xxdim
         dswrf0015(xx,yy,tt)=0.0
      enddo
      enddo
   else
      call sub0015(psea0015)
      call sub0015(sp0015)
      call sub0015(uu0015)
      call sub0015(vv0015)
      call sub0015(temp0015)
      call sub0015(rh0015)
      call sub0015(r1h0015)
      call sub0015(ncld_upper0015)
      call sub0015(ncld_mid0015)
      call sub0015(ncld_low0015)
      call sub0015(clda0015)
      call sub0015(dswrf0015)
   endif
enddo
num=1
do tt=1,out1633ttdim
   call sub1633(psea1633)
   call sub1633(sp1633)
   call sub1633(uu1633)
   call sub1633(vv1633)
   call sub1633(temp1633)
   call sub1633(rh1633)
   call sub1633(r1h1633)
   call sub1633(ncld_upper1633)
   call sub1633(ncld_mid1633)
   call sub1633(ncld_low1633)
   call sub1633(clda1633)
   call sub1633(dswrf1633)
enddo
num=1
do tt=1,out3439ttdim
   call sub3439(psea3439)
   call sub3439(sp3439)
   call sub3439(uu3439)
   call sub3439(vv3439)
   call sub3439(temp3439)
   call sub3439(rh3439)
   call sub3439(r1h3439)
   call sub3439(ncld_upper3439)
   call sub3439(ncld_mid3439)
   call sub3439(ncld_low3439)
   call sub3439(clda3439)
   call sub3439(dswrf3439)
enddo
num=1
do tt=1,out4051ttdim
   call sub4051(psea4051)
   call sub4051(sp4051)
   call sub4051(uu4051)
   call sub4051(vv4051)
   call sub4051(temp4051)
   call sub4051(rh4051)
   call sub4051(r1h4051)
   call sub4051(ncld_upper4051)
   call sub4051(ncld_mid4051)
   call sub4051(ncld_low4051)
   call sub4051(clda4051)
   call sub4051(dswrf4051)
enddo
!write

do tdata=1,ttdim
do yy=1,yydim
do xx=1,xxdim
   if (tdata <= out0015ttdim)then
      psea(xx,yy,tdata)=psea0015(xx,yy,tdata)
      sp(xx,yy,tdata)=sp0015(xx,yy,tdata)
      uu(xx,yy,tdata)=uu0015(xx,yy,tdata)
      vv(xx,yy,tdata)=vv0015(xx,yy,tdata)
      temp(xx,yy,tdata)=temp0015(xx,yy,tdata)
      rh(xx,yy,tdata)=rh0015(xx,yy,tdata)
      r1h(xx,yy,tdata)=r1h0015(xx,yy,tdata)
      ncld_upper(xx,yy,tdata)=ncld_upper0015(xx,yy,tdata)
      ncld_mid(xx,yy,tdata)=ncld_mid0015(xx,yy,tdata)
      ncld_low(xx,yy,tdata)=ncld_low0015(xx,yy,tdata)
      clda(xx,yy,tdata)=clda0015(xx,yy,tdata)
      dswrf(xx,yy,tdata)=dswrf0015(xx,yy,tdata)
   elseif(tdata <= out0015ttdim+out1633ttdim)then
      tt=tdata - out0015ttdim
      psea(xx,yy,tdata)=psea1633(xx,yy,tt)
      sp(xx,yy,tdata)=sp1633(xx,yy,tt)
      uu(xx,yy,tdata)=uu1633(xx,yy,tt)
      vv(xx,yy,tdata)=vv1633(xx,yy,tt)
      temp(xx,yy,tdata)=temp1633(xx,yy,tt)
      rh(xx,yy,tdata)=rh1633(xx,yy,tt)
      r1h(xx,yy,tdata)=r1h1633(xx,yy,tt)
      ncld_upper(xx,yy,tdata)=ncld_upper1633(xx,yy,tt)
      ncld_mid(xx,yy,tdata)=ncld_mid1633(xx,yy,tt)
      ncld_low(xx,yy,tdata)=ncld_low1633(xx,yy,tt)
      clda(xx,yy,tdata)=clda1633(xx,yy,tt)
      dswrf(xx,yy,tdata)=dswrf1633(xx,yy,tt)
   elseif(tdata <= out0015ttdim+out1633ttdim+out3439ttdim)then
      tt=tdata - out0015ttdim - out1633ttdim
      psea(xx,yy,tdata)=psea3439(xx,yy,tt)
      sp(xx,yy,tdata)=sp3439(xx,yy,tt)
      uu(xx,yy,tdata)=uu3439(xx,yy,tt)
      vv(xx,yy,tdata)=vv3439(xx,yy,tt)
      temp(xx,yy,tdata)=temp3439(xx,yy,tt)
      rh(xx,yy,tdata)=rh3439(xx,yy,tt)
      r1h(xx,yy,tdata)=r1h3439(xx,yy,tt)
      ncld_upper(xx,yy,tdata)=ncld_upper3439(xx,yy,tt)
      ncld_mid(xx,yy,tdata)=ncld_mid3439(xx,yy,tt)
      ncld_low(xx,yy,tdata)=ncld_low3439(xx,yy,tt)
      clda(xx,yy,tdata)=clda3439(xx,yy,tt)
      dswrf(xx,yy,tdata)=dswrf3439(xx,yy,tt)
   else
      tt=tdata - out0015ttdim - out1633ttdim - out3439ttdim
      psea(xx,yy,tdata)=psea4051(xx,yy,tt)
      sp(xx,yy,tdata)=sp4051(xx,yy,tt)
      uu(xx,yy,tdata)=uu4051(xx,yy,tt)
      vv(xx,yy,tdata)=vv4051(xx,yy,tt)
      temp(xx,yy,tdata)=temp4051(xx,yy,tt)
      rh(xx,yy,tdata)=rh4051(xx,yy,tt)
      r1h(xx,yy,tdata)=r1h4051(xx,yy,tt)
      ncld_upper(xx,yy,tdata)=ncld_upper4051(xx,yy,tt)
      ncld_mid(xx,yy,tdata)=ncld_mid4051(xx,yy,tt)
      ncld_low(xx,yy,tdata)=ncld_low4051(xx,yy,tt)
      clda(xx,yy,tdata)=clda4051(xx,yy,tt)
      dswrf(xx,yy,tdata)=dswrf4051(xx,yy,tt)
   endif
enddo
enddo
enddo
!write(*,*)ww(:,:,:,:)
open(20,file="msm.bin",status="unknown",access="direct",form="unformatted",recl=4*xxdim*yydim)
rnum=1
do tt=1,ttdim
   write(20,rec=rnum)psea(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)sp(:,:,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)rh(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)r1h(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)ncld_upper(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)ncld_mid(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)ncld_low(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)clda(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)dswrf(:,:,tt)
   rnum=rnum+1
enddo
close(20)

stop

contains

subroutine sub0015(var)
  implicit none
  real::var(xxdim,yydim,out0015ttdim),fread(out0015recl/4)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,tt)=out0015(num)
        num=num+1
  enddo
  enddo
end subroutine 
subroutine sub1633(var)
  implicit none
  real::var(xxdim,yydim,out1633ttdim),fread(out1633recl/4)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,tt)=out1633(num)
        num=num+1
  enddo
  enddo
end subroutine 
subroutine sub3439(var)
  implicit none
  real::var(xxdim,yydim,out3439ttdim),fread(out3439recl/4)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,tt)=out3439(num)
        num=num+1
  enddo
  enddo
end subroutine 
subroutine sub4051(var)
  implicit none
  real::var(xxdim,yydim,out4051ttdim),fread(out4051recl/4)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,tt)=out4051(num)
        num=num+1
  enddo
  enddo
end subroutine 
  
end

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

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

気象業務センターで確認すると、MSMデータは、北緯22.4~47.6度と東経120~150度で、地表面データでは、x方向(経度)格子点間隔0.0625の格子点数481、y方向(緯度)格子間隔0.05度の格子点数505です。

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

地表面データのコントロールファイル

DSET msm.bin
UNDEF 999999999E+20
XDEF 481 LINEAR 120 0.0625
YDEF 505 LINEAR 22.4 0.05
ZDEF 1 LEVELS 1
TDEF 51 LINEAR 00z01jan0000 1hr
VARS 12
psea 1 99 sp
sp 1 99 u
u 1 99 u
v 1 99 v
temp 1 99 temp
rh 1 99 rh
r1h 1 99 prec
ncld_upper 1 99 ncld_upper
ncld_mid 1 99 ncld_mid
ncld_low 1 99 ncld_low
clda 1 99 clda
dswrf 1 99 dswrf
ENDVARS

GrADSで描画

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

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

open msm_s.ctl

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

ga ->d temp

GrADSの詳細はこちらで

前の記事

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

次の記事

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