気象庁全球数値予報モデルGPV (GSM)

気象庁全球数値予報モデル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)