TRMM 3B42 降水データ

人工衛星TRMMの降水データ3B42をプロットする。
この降水データは水平分解能0.25度で3時間ごとにデータがある。
このデータをNASA MiradorからOPeNDAPを通して入手する。
以下の例は7月27日0時世界標準時（日本標準時9時）の日本付近のデータをプロットする例である。この日は東シナ海に熱帯低気圧が存在し、それによる雨が見える。
OPeNDAPからはPyNIOを使ってデータを入手する。
関数name_makerはファイルのあるURLはこういう規則だろうと推測して作ったものである。間違いかもしれないので確認すること。

import Nio
import numpy as np

def read_TRMM_3B42(date,lon1,lon2,lat1,lat2):
    """
   read TRMM 3B42 3 hourly rainfall  from NASA Mirador

   usage:
      lon,lat,sst=read__TRMM_3B42(date,lon1,lon2,lat1,lat2)
     
   input
     
      date (year,month,day,hour) 
      lon1,lon2  Western,East bounday
      lat1,lat2  South,West bounday
   output
      lon longitude
      lat latitude
      rainfall data (mm/hour)         
    """
#check
    assert date[3] % 3 ==0, "Error, hour must be multiple of 3!!!"
# open OPENDAP
    fname=name_maker(date[0],date[1],date[2],date[3])
    f = Nio.open_file(fname)
## variable name check
    #variables = f.variables.keys()
    #print f
    #print "variables",variables

# read data
    temp=f.variables["DATA_GRANULE.PlanetaryGrid.precipitation"]

#grid
    # grid
    lon=-179.875+np.arange(temp.shape[1])*0.25
    lat=-49.875 +np.arange(temp.shape[2])*0.25
    # selection
    xrange=np.where(np.logical_and(lon>=lon1,lon<=lon2))
    yrange=np.where(np.logical_and(lat>=lat1,lat<=lat2))
    
#precipitation
    x=xrange[0]
    y=yrange[0]
    precip=(temp[0,x[0]:x[-1]+1,y[0]:y[-1]+1].T).reshape((len(y),len(x)))

# closing 
    f.close()
    return lon[xrange],lat[yrange],precip

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

def name_maker(year,month,day,hour):
    import scikits.timeseries as ts
    #
    t=ts.Date('H', year=year,month=month,day=day,hour=hour)
    if t.hour==0:
        tc=t-1
    else:
        tc=t
    cyear = "%04d"  % tc.year
    cyear2= "%02d"  % (t.year-2000)
    cmonth2= "%02d" % t.month
    cday2=   "%02d" % t.day
    shour = str(t.hour)
    yd3= "%03d" % tc.day_of_year
    #
    fname="http://disc2.nascom.nasa.gov/opendap/TRMM_L3//TRMM_3B42//%(cyear)s/%(yd3)s/3B42.%(cyear2)s%(cmonth2)s%(cday2)s.%(shour)s.6A.HDF.Z" %locals()
    return fname

############################################
if __name__ == '__main__':
    import matplotlib.pyplot as plt
    from mpl_toolkits.basemap import Basemap
    # example plot

# range
    lon1=120.
    lon2=145.
    lat1=20.
    lat2=40.

#date
    year=2009
    month=7
    day=27
    hour=0 
    lon,lat,precip=read_TRMM_3B42(date=(year,month,day,hour),lon1=lon1,lon2=lon2,lat1=lat1,lat2=lat2)

   #plot
    vmax=int(precip.max())
    vmin=int(precip.min())
    fig=plt.figure()
    ax = fig.add_axes([0.1,0.1,0.8,0.8])
    map = Basemap(projection='merc',llcrnrlat=lat1,urcrnrlat=lat2,\
                      llcrnrlon=lon1,urcrnrlon=lon2,resolution='i')
    lons,lats=np.meshgrid(lon[:],lat[:])
    x,y=map(lons,lats)
    map.drawcoastlines()
    map.drawmapboundary()
    map.pcolor(x,y,precip,vmin=vmin,vmax=vmax,cmap=plt.cm.Blues)
    ti=str(year)+"/"+str(month)+"/"+str(day)+":"+str(hour)+"Z"
    plt.title(ti)
# colorbar
    pos = ax.get_position()
    l, b, w, h = pos.bounds
    cax = plt.axes([l+w+0.02, b, 0.03, h]) # setup colorbar axes.
    plt.colorbar(cax=cax) # draw colorbar
    plt.savefig("TRMM_3B42.png")
    plt.show()

参考
NASA Giovanniのサイト
http://disc2.nascom.nasa.gov/Giovanni/tovas/TRMM_V6.3B42.2.shtml

TRMM 3B42 readme
ftp://meso-a.gsfc.nasa.gov/pub/trmmdocs/3B42_3B43_doc.pdf