EORC MODIS SST

JAXA EORC提供のMODIS SSTをプロットするスクリプトを作成した。
データのバイナリはEORCに提供していただいた。
下の例ではは2010年8月2日の領域６のデータから紀伊半島付近をプロットしている。 A2GL11008020431OD1_OSTAQ_01000_01000_sst_06
ファイルのフォーマットの情報はここ。

import numpy as np

def read_MODIS_EORC(file,lon1,lon2,lat1,lat2):
    """
   read MODIS_SST from EORC

   usage:
      lon,lat,sst=read_MODIS_EORC(file,lon1,lon2,lat1,lat2)
     
   input
     
      file filename 
      lon1,lon2  Western,East bounday
      lat1,lat2  South,West bounday
   output
      lon longitude
      lat latitude
      SST sst data (masked array) in C          
    """
    # open file
    f=open(file, 'r')
    # read header
    pixel=f.read(6); print "pixel=",pixel
    line=f.read(6); print "line=",line
    nwlon=f.read(8); print "NW lon.=",nwlon
    nwlat=f.read(8); print "NW lat.=",nwlat
    res=f.read(8);   print "resolution=",res
    slope=f.read(9);   print "slope=",slope
    offset=f.read(9);  print "ofset=",offset
    dummy=f.read(1)
    parameter=f.read(8); print "parameter name=",parameter
    dummy=f.read(1)
    filename=f.read(45); print "file name=",filename
    f.close()
    # read data
    f=open(file,'r')
    f.seek(int(pixel)*2) #skip header
    data=np.fromfile(f,dtype=">u2")
    data=data.reshape((int(line),int(pixel)))
    msk=np.logical_or(data==0,data>=65534) # 0 cloud, 65534 out of obs., 65535 land 
    data=np.ma.masked_array(data[::-1],mask=msk[::-1,:])
    data=data*float(slope)+float(offset)-273.15 # to degree C
    f.close()
    # grid
    lon=float(nwlon)+np.arange(int(pixel))*float(res)
    lat=float(nwlat)-np.arange(int(line))*float(res)
    lat=lat[::-1]
    # selection
    xrange=np.logical_and(lon>=lon1,lon<=lon2)
    yrange=np.logical_and(lat>=lat1,lat<=lat2)
    sst=(data[yrange,:])[:,xrange] 

    return lon[xrange],lat[yrange], sst

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

    #read data
    filename="A2GL11008020431OD1_OSTAQ_01000_01000_sst_06"
    lon1=135.
    lon2=138.
    lat1=32.
    lat2=35.
    lon,lat,sst=read_MODIS_EORC(filename,lon1,lon2,lat1,lat2)

    #plot
    tmax=int(sst.max())
    tmin=int(sst.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,sst,vmin=tmin,vmax=tmax)
# 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("MODIS_EORC_SST.png")
    plt.show()

謝辞
ここで用いたSSTバイナリは宇宙航空研究開発機構(JAXA)地球観測研究センター(EORC)から提供されたものである。ここに感謝する。

関連ポスト
MODIS aqua Global Level 3 Mapped Mid-IR SST (4km)
http://oceansciencehack.blogspot.com/2010/09/modis-aqua-global-level-3-mapped-mid-ir.html

参考
MODIS EORC SST FAQ
http://kuroshio.eorc.jaxa.jp/ADEOS/mod_nrt_new/html/05_faq.html
JAXA　MODISページ
http://www.eorc.jaxa.jp/hatoyama/satellite/sendata/modis_j.html
MODIS wikipedia
http://ja.wikipedia.org/wiki/MODIS