気象庁全球数値予報モデルGPV (GSM) の初期値データを京都大学のサイトからダウンロードして、清水慎吾さん作成のgsm2binでgrib2形式からgrads形式に変換するモジュール。
以下がそのモジュール。サンプルのメインプログラムでは、２０１１年１月１日6時のデータを変換し、python interface to gradsで海面気圧をプロットしている。
get_gsm.py
def get_gsm(date=[2011,1,1,0],dir="."):
"""
download jma gsm initial data from http://database.rish.kyoto-u.ac.jp/arch/jmadata/da
ta/gpv/original/
and
convert to grads file using gsm2bin by Dr. Shingo Shimizu
(http://shimizus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%E2%A5%C7%A5%EB%B4%D8%CF%A2%A5%E1%A
5%E2)
usage:
get_gsm(date,dir)
input
date: [year,month,day,hour] hour must be multiple of 6
default [2011,1,1,0]
dir : directory to save data
"""
import os,urllib
#check
assert date[3] % 6 ==0, "Error, hour must be multiple of 6!!!"
#names from date
urldir,original,gradsbinname,gradsdate=name_maker(date)
# get original data
if not os.path.exists(dir+"/"+original):
if not os.path.exists(dir): os.mkdir(dir)
# Download the data
print 'Downloading data, please wait'
opener = urllib.urlopen(urldir+original)
open(dir+"/"+original, 'w').write(opener.read())
print 'Downloaded'
# convert to grads file
os.system("gsm2bin %(dir)s/%(original)s %(dir)s/%(gradsbinname)s" %locals())
# make ctlfile
make_ctl(dir,gradsbinname,gradsdate)
############################################
def name_maker(date):
"""
make filenames from date
"""
import datetime
time=datetime.datetime(*date)
cyear=time.strftime("%Y")
cmonth=time.strftime("%m")
cmonth2=(time.strftime("%b")).lower()
cday=time.strftime("%d")
chour=time.strftime("%H")
# URL of original data
urldir="http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/%(cyear)s/%(
cmonth)s/%(cday)s/" %locals()
# File name of original data
original="Z__C_RJTD_%(cyear)s%(cmonth)s%(cday)s%(chour)s0000_GSM_GPV_Rgl_FD0000_grib2.
bin" %locals()
# grads bin name after conversion
gradsbinname="gsm%(cyear)s%(cmonth)s%(cday)s%(chour)s.bin" %locals()
# date in grads format
gradsdate="%(chour)s:00Z%(cday)s%(cmonth2)s%(cyear)s" %locals()
return urldir,original,gradsbinname,gradsdate
###########################################
def make_ctl(dir,gradsbinname,gradsdate):
"""
make grads ctl file for GSM data
"""
#ctl file text
ctltxt="""
dset ^%(gradsbinname)s
options template
undef -999
xdef 720 LINEAR 0.0 0.5
ydef 361 LINEAR -90.0 0.5
zdef 17 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70 50
30 20 10
tdef 1 LINEAR %(gradsdate)s 6hr
vars 16
slp 0 0 sea level pressure
sp 0 0 surface pressure
su 0 0 surface westerly wind
sv 0 0 surface southerly wind
st 0 0 surface temperature
srh 0 0 surface relative humidity
lca 0 0 LowCloudAmount
mca 0 0 MidCloudAmount
hca 0 0 HighCloudAmount
tca 0 0 TotalCloudAmount
z 17 0 height
u 17 0 westerly wind
v 17 0 southerly wind
t 17 0 temperature
rh 17 0 relative humidity
w 17 0 p-velocity
endvars
"""
# write to file
gradsctlname=dir+"/"+gradsbinname.replace(".bin",".ctl")
f=open(gradsctlname,"w")
f.write(ctltxt %locals())
f.close()
#-----------------------------------------------------------------------
if __name__ == '__main__':
# get gsm data
get_gsm(date=[2011,1,1,6],dir="DATA")
# plot
from grads.galab import GaLab # for python interface for grads
ga = GaLab(Bin='grads',Window=False,Echo=True,Verb=0,Port=False)
ga.open("DATA/gsm2011010106.ctl")
script="""
set grads off
set gxout shaded
d slp/100
cbarn
draw title Sea level pressure (hpa)
printim gsm.png white
"""
ga(script)
