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

SHARE

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

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

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

GSMデータの取得

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

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

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

Z__C_RJTD_yyyyMMddhhmmss_GSM_GPV_Rjp_Lsurf_FDddhh-ddhh_grib2.binは、地表面データです。

Z__C_RJTD_yyyyMMddhhmmss_GSM_GPV_Rjp_L-pall_FDddhh-ddhh_grib2.binは、気圧面データです。

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

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_GSM_GPV_Rjp_Lsurf_FD0000-0312_grib2.bin

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

264時間先までほしいので、合計で3つのファイルを取得します。また、初期値が06,18時の場合、132時間先までですので、合計のファイルは、2つです。

$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0315-0512_grib2.bin
$ wget  http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2021/07/17/Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0515-1100_grib2.bin

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

$ wgrib2 Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0000-0312_grib2.bin -no_header -bin out00000312.bin
$ wgrib2 Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0318-0512_grib2.bin -no_header -bin out03150512.bin
$ wgrib2 Z__C_RJTD_20210717000000_GSM_GPV_Rjp_L-pall_FD0518-1100_grib2.bin -no_header -bin out05151100.bin

バイナリファイルの編集

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

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

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

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

$ ls -l
-rw-r--r-- 1 tanaka users 74399512 Aug 20 16:11 out00000312.bin
-rw-r--r-- 1 tanaka users 14032128 Aug 20 16:11 out03150512.bin
-rw-r--r-- 1 tanaka users 38588352 Aug 20 16:13 out05151100.bin

out00000312.binのファイルサイズは74399512byte、out03150512.binでは、14032128byte、out05151100.binでは、38588352byteです。

気象業務センターで確認すると、GSMデータは、北緯20~50度と東経120~150度で、気圧面データでは、x方向(経度)格子点間隔0.25の格子点数121、y方向(緯度)格子間隔0.2度の格子点数151です。

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

out00000312.binは、1時間ごとにデータがあり、t方向に85あります。また、初期場には、積算降水量と日射量がありません。このことを考慮すると、

4×2×151×121×84=62,121,400

4×10×151×121×85=12,278,112

62,121,400+12,278,112=74,399,512

out03150512.binは、3時間ごとにデータがあり、t方向に16あります。

4×12×151×121×16= 14,032,128

out05151100.binは、3時間ごとにデータがあり、t方向に44あります。

4×12×151×121×44= 38,588,352

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

implicit none
real,allocatable::out0312(:),out0512(:),out1100(:)
real,allocatable::psea0312(:,:,:),uu0312(:,:,:),vv0312(:,:,:),temp0312(:,:,:),rh0312(:,:,:),rsh0312(:,:,:)
real,allocatable::ncld_upper0312(:,:,:),ncld_mid0312(:,:,:),ncld_low0312(:,:,:),clda0312(:,:,:),sp0312(:,:,:),dswrf0312(:,:,:)
real,allocatable::psea0512(:,:,:),uu0512(:,:,:),vv0512(:,:,:),temp0512(:,:,:),rh0512(:,:,:),rsh0512(:,:,:)
real,allocatable::ncld_upper0512(:,:,:),ncld_mid0512(:,:,:),ncld_low0512(:,:,:),clda0512(:,:,:),sp0512(:,:,:),dswrf0512(:,:,:)

real,allocatable::psea1100(:,:,:),uu1100(:,:,:),vv1100(:,:,:),temp1100(:,:,:),ww1100(:,:,:),rh1100(:,:,:),rsh1100(:,:,:)
real,allocatable::ncld_upper1100(:,:,:),ncld_mid1100(:,:,:),ncld_low1100(:,:,:),clda1100(:,:,:),sp1100(:,:,:),dswrf1100(:,:,:)

real,allocatable::psea(:,:,:),uu(:,:,:),vv(:,:,:),temp(:,:,:),rh(:,:,:),rsh(:,:,:)
real,allocatable::ncld_upper(:,:,:),ncld_mid(:,:,:),ncld_low(:,:,:),clda(:,:,:),sp(:,:,:),dswrf(:,:,:)

integer::xx,yy,tt,num,tdata,rnum
integer,parameter::xxdim=121,yydim=151
integer,parameter::out0312ttdim=85,out0512ttdim=16,out1100ttdim=44,ttdim=145
integer::out0312recl,out0512recl,out1100recl

out0312recl = 4 * 10 * xxdim * yydim * out0312ttdim + 4 * 2 * xxdim * yydim * ( out0312ttdim - 1 )
out0512recl = 4 * 12 * xxdim * yydim * out0512ttdim 
out1100recl = 4 * 12 * xxdim * yydim * out1100ttdim 

allocate(psea0312(xxdim,yydim,out0312ttdim),uu0312(xxdim,yydim,out0312ttdim),vv0312(xxdim,yydim,out0312ttdim),&
temp0312(xxdim,yydim,out0312ttdim),rh0312(xxdim,yydim,out0312ttdim),rsh0312(xxdim,yydim,out0312ttdim),&
ncld_upper0312(xxdim,yydim,out0312ttdim),ncld_mid0312(xxdim,yydim,out0312ttdim),ncld_low0312(xxdim,yydim,out0312ttdim),&
clda0312(xxdim,yydim,out0312ttdim),sp0312(xxdim,yydim,out0312ttdim),dswrf0312(xxdim,yydim,out0312ttdim))
allocate(psea0512(xxdim,yydim,out0512ttdim),uu0512(xxdim,yydim,out0512ttdim),vv0512(xxdim,yydim,out0512ttdim),&
temp0512(xxdim,yydim,out0512ttdim),rh0512(xxdim,yydim,out0512ttdim),rsh0512(xxdim,yydim,out0512ttdim),&
ncld_upper0512(xxdim,yydim,out0512ttdim),ncld_mid0512(xxdim,yydim,out0512ttdim),ncld_low0512(xxdim,yydim,out0512ttdim),&
clda0512(xxdim,yydim,out0512ttdim),sp0512(xxdim,yydim,out0512ttdim),dswrf0512(xxdim,yydim,out0512ttdim))
allocate(psea1100(xxdim,yydim,out1100ttdim),uu1100(xxdim,yydim,out1100ttdim),vv1100(xxdim,yydim,out1100ttdim),&
temp1100(xxdim,yydim,out1100ttdim),rh1100(xxdim,yydim,out1100ttdim),rsh1100(xxdim,yydim,out1100ttdim),&
ncld_upper1100(xxdim,yydim,out1100ttdim),ncld_mid1100(xxdim,yydim,out1100ttdim),ncld_low1100(xxdim,yydim,out1100ttdim),&
clda1100(xxdim,yydim,out1100ttdim),sp1100(xxdim,yydim,out1100ttdim),dswrf1100(xxdim,yydim,out1100ttdim))

allocate(psea(xxdim,yydim,ttdim),uu(xxdim,yydim,ttdim),vv(xxdim,yydim,ttdim),temp(xxdim,yydim,ttdim))
allocate(rh(xxdim,yydim,ttdim),rsh(xxdim,yydim,ttdim),ncld_upper(xxdim,yydim,ttdim),ncld_mid(xxdim,yydim,ttdim))
allocate(ncld_low(xxdim,yydim,ttdim),clda(xxdim,yydim,ttdim),sp(xxdim,yydim,ttdim),dswrf(xxdim,yydim,ttdim))
allocate(out0312(out0312recl/4),out0512(out0512recl/4),out1100(out1100recl/4))
open(10,file="out00000312.bin",access="direct",recl=out0312recl)
open(11,file="out03150512.bin",access="direct",recl=out0512recl)
open(12,file="out05151100.bin",access="direct",recl=out1100recl)
read(10,rec=1)out0312(:)
read(11,rec=1)out0512(:)
read(12,rec=1)out1100(:)
close(10)
close(11)
close(12)
!read
num=1
do tt=1,out0312ttdim
   if(tt==1)then
      call sub0312(psea0312)
      call sub0312(sp0312)
      call sub0312(uu0312)
      call sub0312(vv0312)
      call sub0312(temp0312)
      call sub0312(rh0312)
      do yy=1, yydim
      do xx=1, xxdim
         rsh0312(xx,yy,tt)=0.0
      enddo
      enddo
      call sub0312(ncld_upper0312)
      call sub0312(ncld_mid0312)
      call sub0312(ncld_low0312)
      call sub0312(clda0312)
      do yy=1, yydim
      do xx=1, xxdim
         dswrf0312(xx,yy,tt)=0.0
      enddo
      enddo
   else
      call sub0312(psea0312)
      call sub0312(sp0312)
      call sub0312(uu0312)
      call sub0312(vv0312)
      call sub0312(temp0312)
      call sub0312(rh0312)
      call sub0312(rsh0312)
      call sub0312(ncld_upper0312)
      call sub0312(ncld_mid0312)
      call sub0312(ncld_low0312)
      call sub0312(clda0312)
      call sub0312(dswrf0312)
   endif
enddo
num=1
do tt=1,out0512ttdim
   call sub0512(psea0512)
   call sub0512(sp0512)
   call sub0512(uu0512)
   call sub0512(vv0512)
   call sub0512(temp0512)
   call sub0512(rh0512)
   call sub0512(rsh0512)
   call sub0512(ncld_upper0512)
   call sub0512(ncld_mid0512)
   call sub0512(ncld_low0512)
   call sub0512(clda0512)
   call sub0512(dswrf0512)
enddo
num=1
do tt=1,out1100ttdim
   call sub1100(psea1100)
   call sub1100(sp1100)
   call sub1100(uu1100)
   call sub1100(vv1100)
   call sub1100(temp1100)
   call sub1100(rh1100)
   call sub1100(rsh1100)
   call sub1100(ncld_upper1100)
   call sub1100(ncld_mid1100)
   call sub1100(ncld_low1100)
   call sub1100(clda1100)
   call sub1100(dswrf1100)
enddo

!write

do tdata=1,ttdim
do yy=1,yydim
do xx=1,xxdim
   if (tdata <= out0312ttdim)then
      psea(xx,yy,tdata)=psea0312(xx,yy,tdata)
      uu(xx,yy,tdata)=uu0312(xx,yy,tdata)
      vv(xx,yy,tdata)=vv0312(xx,yy,tdata)
      temp(xx,yy,tdata)=temp0312(xx,yy,tdata)
      rh(xx,yy,tdata)=rh0312(xx,yy,tdata)
      rsh(xx,yy,tdata)=rsh0312(xx,yy,tdata)
      ncld_upper(xx,yy,tdata)=ncld_upper0312(xx,yy,tdata)
      ncld_mid(xx,yy,tdata)=ncld_mid0312(xx,yy,tdata)
      ncld_low(xx,yy,tdata)=ncld_low0312(xx,yy,tdata)
      clda(xx,yy,tdata)=clda0312(xx,yy,tdata)
      sp(xx,yy,tdata)=sp0312(xx,yy,tdata)
      dswrf(xx,yy,tdata)=dswrf0312(xx,yy,tdata)
   elseif(tdata <= out0312ttdim+out0512ttdim)then
      tt=tdata - out0312ttdim
      psea(xx,yy,tdata)=psea0512(xx,yy,tt)
      uu(xx,yy,tdata)=uu0512(xx,yy,tt)
      vv(xx,yy,tdata)=vv0512(xx,yy,tt)
      temp(xx,yy,tdata)=temp0512(xx,yy,tt)
      rh(xx,yy,tdata)=rh0512(xx,yy,tt)
      rsh(xx,yy,tdata)=rsh0512(xx,yy,tt)
      ncld_upper(xx,yy,tdata)=ncld_upper0512(xx,yy,tt)
      ncld_mid(xx,yy,tdata)=ncld_mid0512(xx,yy,tt)
      ncld_low(xx,yy,tdata)=ncld_low0512(xx,yy,tt)
      clda(xx,yy,tdata)=clda0512(xx,yy,tt)
      sp(xx,yy,tdata)=sp0512(xx,yy,tt)
      dswrf(xx,yy,tdata)=dswrf0512(xx,yy,tt)
   else
      tt=tdata - out0312ttdim - out0512ttdim
      psea(xx,yy,tdata)=psea1100(xx,yy,tt)
      uu(xx,yy,tdata)=uu1100(xx,yy,tt)
      vv(xx,yy,tdata)=vv1100(xx,yy,tt)
      temp(xx,yy,tdata)=temp1100(xx,yy,tt)
      rh(xx,yy,tdata)=rh1100(xx,yy,tt)
      rsh(xx,yy,tdata)=rsh1100(xx,yy,tt)
      ncld_upper(xx,yy,tdata)=ncld_upper1100(xx,yy,tt)
      ncld_mid(xx,yy,tdata)=ncld_mid1100(xx,yy,tt)
      ncld_low(xx,yy,tdata)=ncld_low1100(xx,yy,tt)
      clda(xx,yy,tdata)=clda1100(xx,yy,tt)
      sp(xx,yy,tdata)=sp1100(xx,yy,tt)
      dswrf(xx,yy,tdata)=dswrf1100(xx,yy,tt)
   endif
enddo
enddo
enddo
open(20,file="gsm.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)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)rsh(:,:,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)sp(:,:,tt)
   rnum=rnum+1
   write(20,rec=rnum)dswrf(:,:,tt)
   rnum=rnum+1
enddo
close(20)

stop

contains

subroutine sub0312(var)
  implicit none
  real::var(xxdim,yydim,out0312ttdim)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,tt)=out0312(num)
        num=num+1
  enddo
  enddo
end subroutine 
subroutine sub0512(var)
  implicit none
  real::var(xxdim,yydim,out0512ttdim)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,tt)=out0512(num)
        num=num+1
  enddo
  enddo
end subroutine 
subroutine sub1100(var)
  implicit none
  real::var(xxdim,yydim,out1100ttdim)
  do yy=1,yydim
  do xx=1,xxdim
        var(xx,yy,tt)=out1100(num)
        num=num+1
  enddo
  enddo
end subroutine 
  
end

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

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

気象業務センターで確認すると、日本域のGSMデータは、日本とその近海(北緯20度～50度、東経120度〜150度)で、x方向(経度)格子間隔0.25度の格子点数121個、y方向(緯度)格子間隔0.2度の格子点数151個です。

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

GSM(日本域)

DSET gsm.bin
UNDEF 999999999E+20
XDEF 121 LINEAR 120 0.25
YDEF 151 LINEAR 20 0.2
ZDEF 1 LEVELS 1000
TDEF 145 LINEAR 00z01jan0000 3hr
VARS 12
psea 1 99 psea
u 1 99 u
v 1 99 v
temp 1 99 temp
rh 1 99 rh
rsh 1 99 rsh
ncld_upper 1 99 ncld_upper
ncld_mid 1 99 ncld_mid
ncld_low 1 99 ncld_low
clda 1 99 clda
sp 1 99 sp
dswrf 1 99 dswrf
ENDVARS

GrADSで描画

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

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

open gsm_s.ctl

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

ga ->d temp

GrADSの詳細はこちらで

前の記事

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

次の記事

PHPとLaravelの環境構築