2010年12月17日金曜日

気象庁全国合成レーダー降雨



気象庁の全国合成レーダー降雨データがこちらで配布されている。
データは10分ごとにある。
データはtarでアーカイブされており、その中に含まれるデータはgrib2形式である。
清水慎吾さんがこちらで、データをgrads形式に変換するプログラムを公開している(全国合成気象庁レーダ雨量,GRIB2形式)。

以下のプログラムはファイルをネットから読み込んで、清水さんのツールでgrads形式に変換するpythonスクリプトである。

RadarToGradsはファイル(tarファイル)をネットから呼んで、tarファイルから降雨データのバイナリファイル(grib2)を取り出し、gradsファイルを作る。
grib2からgradsへのファイル変換は清水さんのjmaradar2binで行っている。gradsコントロールファイル("view_radar.ctl")が自動的にできるようにmain.cで"gradsout=1"としている。

PlotByGradsは作ったgradsデータを使って、gradsでプロットしている。gradsを操作するのに、pythonインターフェイスを使っている。

メインでは、これらを使って、2010年10月20日午前8時00分世界標準時(日本時間午後5時)の降水を東経125-135度北緯24-35度でgradsデータに変換して、それをプロットしている。これは奄美大島で豪雨が発生したころにあたる。

import_jmaradar.py
def RadarToGrads(dir=".",gradsbinname="temp.bin",year=2010,mon=1,day=1,hour=0,mi
n=0,slon=118.,elon=150.,slat=20.,elat=48.):
    """
   Import JMA rain radar data (every 10min) from http://database.rish.kyoto-u.ac
.jp/arch/jmadata/synthetic-original.html
   Then convert data to grads file using tool by Dr. Shingo Shimizu (http://shim
izus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%EC%A1%BC%A5%C0%B4%D8%CF%A2%A5%E1%A5%E2)
 

   usage:
       RadarToGrads(dir=".",gradsbinname="temp.bin",year=2010,mon=1,day=1,hour=0
,min=0,slon=118.,elon=150.,slat=20.,elat=48.):
     
   input
      dir: directory to import tar file
      gradsbinname: file name for grads binary
      year, mon, day, hour, min  : time for data (every 10 min) in UTC
      slon,elon  Western,East bounday
      slat,elat  South,North bounday


   output files
      grads binary: gradsbinname
      grads control: view_radar.ctl        
         
    """


    import os
    import tarfile
    import urllib

    # filenames
    syear="%04d" % year
    smon="%02d" % mon
    sday="%02d" % day
    shour="%02d" % hour
    smin="%02d" % min
    date=syear+smon+sday+shour+smin
    tarname="Z__C_RJTD_%(date)s00_RDR_JMAGPV__grib2.tar" %locals()
    urldir="http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/jma-radar/synth
etic/original/%(syear)s/%(smon)s/%(sday)s/" %locals()
    filename="Z__C_RJTD_%(date)s00_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin" %l
ocals()

    # get tar file
    if not os.path.exists(dir+"/"+tarname):
        if not os.path.exists(dir): os.mkdir(dir) 
        # Download the data
        print 'Downloading data, please wait'
        opener = urllib.urlopen(urldir+tarname )
        open(dir+"/"+tarname, 'w').write(opener.read())
        print 'Downloaded'

    # extract radar precipitation file
    tarfile.open(dir+"/"+tarname, 'r').extract(filename)

    # convert to grads file
    os.system("jmaradar2bin %(filename)s %(gradsbinname)s %(slon)s %(elon)s %(sl
at)s %(elat)s" %locals())
    os.system("rm %(filename)s" %locals())

####################################################################

def PlotByGrads(pngname):
    """
      plot converted JMA radar rain data
    """
    from grads.galab import GaLab  # for python interface for grads
    ga = GaLab(Bin='grads',Window=False,Echo=False,Verb=0,Port=False)
    ga.open("view_radar.ctl")
    script="""
    set grads off
    set mpdset hires
    set gxout shaded
    set clevs 0 1 5 10 20 30 50 80
    d rr
    cbarn
    printim %(pngname)s white
"""
    ga(script %locals())


#####################################################################

if __name__ == '__main__':

# radar file to grads 
    RadarToGrads(dir="DATA",gradsbinname="jmaradar.bin",year=2010,mon=10,day=20,
hour=8,min=00,slon=125.,elon=135.,slat=24.,elat=35.)

# plot by grads
    PlotByGrads("jmaradar.png")


0 件のコメント: