MSMの描画

気象庁が作成しているMSM(Meso scale model)による数値予報をGrADSで描画します。図はいつか忘れましたが、ある時刻を初期値とした今後24時間降水量と海面気圧と地上風をかいたものです。

データは、ここでは簡単のため、京大生存研にてnetcdfに変換されたデータ
http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/latest/
を用います。元データはgrib2形式のため何らかの形で変換する必要があります。これは現在こちらのページではサポートしていません。そのためGSMは京大生存研にてnetcdfに変換するサービスがないのでかけないことになリます。

ともかく、データをダウンロードします。
MSM(日付)P.ncがpressure level（風、気温など）
MSM(日付)S.ncがsurface level（地表気温、降水量など）
が入っています。COARSE規約でかかれているのでsdfopenで読めます。

netcdfですので、データさえ手に入れば描画するのは比較的容易だと思われますので、とりあえずここまででかきかけとします。

MSMで使用される地形データ。TOPO.MSMから描画できます。

GMTで描画する場合

GMTで描画する場合は、grd2xyzを使います。ただし3次元データにはこのままだとできないので、surfaceデータしかいまのところは対応していません。
たとえば、以下のような天気図がかけます。

スクリプトは 
#!/bin/sh
file=MSM2013110921S
region=120/150/23/47
proj=q1:14000000
size=s0.07
#region=125/148/30/46
#proj=q1:11000000
#size=s0.08
#region=136/141/34/38
#proj=q1:2500000
#size=s0.4
out=test.eps
cpt=cloudamount.cpt

gmtset ANNOT_FONT_SIZE_PRIMARY 20
gmtset BASEMAP_TYPE plain
t=13

echo 0 255 255 255 1 50 50 153 >tmp.cpt
makecpt -Ctmp.cpt -T0/100/10 >$cpt

grd2xyz $file.nc?psea[$t] > psea2.dat
awk '{print $1,$2,$3*0.01}' psea2.dat >psea.dat
grd2xyz $file.nc?u[$t] > u.dat
grd2xyz $file.nc?v[$t] > v.dat
grd2xyz $file.nc?ncld_low[$t] > clda.dat
xyz2grd u.dat -Gu.grd -I0.5/0.6 -R$region
xyz2grd v.dat -Gv.grd -I0.5/0.6 -R$region
xyz2grd psea.dat -Gpsea.grd -I0.0625/0.05 -R$region


psxy -J$proj -R$region clda.dat -S$size -C$cpt -K -BSWneg5a5 -X2 -Y1.5>$out
psscale -D24.5/9.3/18.5/0.4 -K -L -E -N -O -C$cpt>>$out
pscoast -J$proj -R$region -W3,153/51/0 -K -O -Di>>$out
grdcontour psea.grd -J -R -C2 -W2,black -K -O>>$out
grdcontour psea.grd -J -R -C10 -A10 -W4,black -K -O>>$out
grdvector -J -R u.grd v.grd -K -O -Gblack -Q0.015/0.04/0.04 -S50 >>$out

convert test.eps test.png
convert -rotate 90 test.png test.png
ついにできた。fortranでデータ読み込みを制御することでpressure dataからも同時に生成。
この図は海面気圧、850気温、500高度。

スクリプトは 
#!/bin/sh
day=2013112200
files=MSM"$day"S.nc
filep=MSM"$day"P.nc
region=120/150/23/47
proj=q1:14000000
size=s0.07
out=test.eps
cpt=cloudamount.cpt
fortran=read850500.f90
timep=1 # x3h
times=`expr $timep \* 3 - 3`
echo $times
echo program netcdf_example >$fortran
echo  implicit none >>$fortran
echo  include "'"/usr/include/netcdf.inc"'" >>$fortran
echo  integer, parameter :: xsz=241, ysz=253, zsz=16 ,tsz=12  >>$fortran
echo  integer :: ncid, status, rhid >>$fortran
echo " real(4) :: t850(xsz,ysz,zsz,tsz)"  >>$fortran
echo " real(4) :: z500(xsz,ysz,zsz,tsz)" >>$fortran
echo " real(4) :: lon(xsz),lat(ysz)" >>$fortran
echo character :: filename*60 >>$fortran
echo "character(20) :: var1,var2,var3,var4" >>$fortran
echo integer :: i,j,k,l,m,n >>$fortran
echo var1="'"z"'" >>$fortran
echo var2="'"temp"'" >>$fortran
echo var3="'"lon"'" >>$fortran
echo var4="'"lat"'" >>$fortran
echo  filename="'"$filep"'" >>$fortran
echo " status = nf_open(filename, nf_nowrite, ncid)">>$fortran
echo " status = nf_inq_varid(ncid, var1, rhid)" >>$fortran
echo " status = nf_get_var_real(ncid, rhid, z500)">>$fortran
echo " status = nf_inq_varid(ncid, var2, rhid)">>$fortran
echo " status = nf_get_var_real(ncid, rhid, t850)">>$fortran
echo "status = nf_inq_varid(ncid, var3, rhid)" >>$fortran
echo "status = nf_get_var_real(ncid, rhid, lon)" >>$fortran
echo "status = nf_inq_varid(ncid, var4, rhid)" >>$fortran
echo "status = nf_get_var_real(ncid, rhid, lat)" >>$fortran
echo " status = nf_close(ncid)" >>$fortran
echo "open(91,file="'"z500.dat"'")" >>$fortran
echo "open(92,file="'"t850.dat"'")" >>$fortran
echo do j=1,ysz >>$fortran
echo do i=1,xsz >>$fortran
echo "write(91,*) lon(i),lat(j),z500(i,j,10,"$timep")" >>$fortran
echo "write(92,*) lon(i),lat(j),t850(i,j,6,"$timep")*0.002613+255.40-273.15" >>$fortran
echo end do >>$fortran
echo end do >>$fortran
echo  stop>>$fortran
echo end program>>$fortran
gfortran $fortran -L/usr/local/netcdf-4.3.0/lib -lnetcdf -L/usr/lib -lnetcdff
./a.out

gmtset ANNOT_FONT_SIZE_PRIMARY 20
gmtset BASEMAP_TYPE plain

makecpt -Cwysiwyg -T4900/5900/50 >z500.cpt

grd2xyz $files?psea[$t] > psea2.dat
awk '{print $1,$2,$3*0.01}' psea2.dat >psea.dat
xyz2grd psea.dat -Gpsea.grd -I0.125/0.1 -R$region
xyz2grd z500.dat -Gz500.grd -I0.125/0.1 -R$region
xyz2grd t850.dat -Gt850.grd -I0.125/0.1 -R$region
psxy -J$proj -R$region clda.dat -S$size -C$cpt -K -BSWneg5a5 -X2 -Y1.5>$out
psscale -D24.5/9.3/18.5/0.4 -K -L -E -N -O -Cz500.cpt>>$out
grdimage z500.grd -J -R -Cz500.cpt -K -O >>$out
grdcontour t850.grd -J -R -C2 -A2+s16 -W2ta -K -O>>$out
grdcontour psea.grd -J -R -C2 -W4,black -K -O>>$out
grdcontour psea.grd -J -R -C10 -A10+s16 -W6,black -K -O>>$out
pscoast -J$proj -R$region -W3,153/51/0 -K -O -Di -BSWneg5a5>>$out

convert test.eps test.png
convert -rotate 90 test.png test.png
スクリプト内でいちいちechoするので頭悪いですが。
本当はtempはshort=>integer(2)で定義されているのでinteger(2)で読むべきなのですがreal(4)でうまく読めてしまったのでこうしています。ncdump -cしたときに、 
	short temp(time, p, lat, lon) ;
		temp:scale_factor = 0.002613491f ;
		temp:add_offset = 255.4005f ;
		temp:long_name = "temperature" ;
		temp:units = "K" ;
		temp:standard_name = "air_temperature" ;
と出ているので、正しい気温に戻すためには0.002613491をかけて255.4をたさなければいけないことに注意してください。

降水量 
#!/bin/sh
day=$1
file=MSM"$1"S.nc
region=120/150/22.4/47.6
proj=q135/1/1:12500000
size=s0.084
#region=128/140/30/39
#proj=q1:5600000
#size=s0.18
cpt=rain.cpt
gmtset ANNOT_FONT_SIZE_PRIMARY 24
gmtset ANNOT_FONT_PRIMARY 0
t=1
out=test.eps
echo 1	50 50	255	5	50 50	255 >$cpt
echo 5	0	128	255	10	0	128	255>>$cpt
echo 10	64	192	255	20	64	192	255>>$cpt
echo 20	64	255	255	30	64	255	255>>$cpt
echo 30	128	255	64	50	128	255	64>>$cpt
echo 50	255	255	64	100	255	255	64>>$cpt
echo 100	255	160	64	200	255	160	64>>$cpt
echo 200	255	32	64	400	255	32	64>>$cpt
echo B 255 255 255 >>$cpt
echo F 180 0 0 >>$cpt

echo 45 > snowtemp.dat
gfortran rain.f90
rm -f temp.dat
rm -f r1h.dat
rm -f snowr1h.dat
while [ $t -lt 34 ]
do
echo $t
grd2xyz $file?temp[$t] >> temp.dat
grd2xyz $file?r1h[$t] >> r1h.dat
t=`expr $t + 1`
done
./a.out
gmtset GRID_PEN_PRIMARY  0.25p,,0.25_2:0
echo 0 0 | psxy -Jx1 -R0/1/0/1 -Ss100 -G255/255/255 -K -X2 -Y-2 -N >test.eps 

gmtset ANNOT_FONT_PRIMARY 31
gmtset ANNOT_FONT_SIZE_PRIMARY 28
psscale -D28/7.5/12/0.8 -K -L -E1 -N -O -C$cpt>>test.eps
gmtset ANNOT_FONT_PRIMARY 0
pscoast -J$proj -R$region -G200/200/200 -K -O -Di>>test.eps
gmtset ANNOT_FONT_SIZE_PRIMARY 20
psxy -J$proj -R$region snowr1h.dat -S$size -C$cpt -K -BSWneg2a2/g2a2 -O>>test.eps
pscoast -J$proj -R$region -W3,40/40/40 -K -O -Dl>>test.eps
grd2xyz $file?psea[33] > psea.dat
awk '{print $1,$2,$3*0.01}' psea.dat >tmp
xyz2grd tmp -I0.125/0.1 -R$region -Gtmp.grd
grdcontour tmp.grd -J$proj -R$region -C1.0 -A4+s12+f1 -K -O>>test.eps
echo 21 -2 15 1.2 | psxy -Jx1 -R0/15/0/15 -Sr -G220/220/220 -K -O -N>>$out
echo 21 -2 22 0 31 "MC" 33hr Prcp.[mm] Init=$day"UTC" | pstext -J -R -K -O -N>>$out
cp test.eps test2.eps
awk '{if($1=="%%BoundingBox:" || $1=="%%HiResBoundingBox:"){print $1, $2,$3,720,950}else{print $0}}' < test2.eps > test.eps

convert test.eps test.png
convert -rotate 90 test.png test.png
back