CEReS NOAA/AVHRR SST

CEReSで配布しているNOAA AVHRR SSTをプロットしてみた。

FTPサイトからデータをダウンロードして、解凍する。
ここでは n1810080316.sst.giを利用する。
(NOAA 18 ,2010年,8月3日,16世界標準時）
以下がスクリプト。
このデータには雲や陸地のマスクはないようだ。
前回はPyNGLでプロットしたが、今回はmatplotlib basemap toolkitを利用した。
matlabのpplot関数に対応する関数を使いたかったからだ（PyNGLにはあるのかな？）。

import numpy as np
import matplotlib.pyplot as plt

def read_CEReS_SST(dir,file,lon1,lon2,lat1,lat2):
    """
   read CERes SST
   usage:
      lon,lat,sst=read_CEReS_SST(dir,file,lon1,lon2,lat1,lat2)     
   input
      dir working dir
      file filename 
           gzipped file is fine.
           URL also works
      lon1,lon2  Western,East bounday
      lat1,lat2  South,West bounday
   output
      lon longitude
      lat latitude
      SST sst data (masked array) in C          
    """
    # read data
    ds=np.DataSource(dir)
    f=ds.open(file, 'r')
    f.seek(80) # skip header
    ssttemp=np.frombuffer(f.read(),dtype=">i2",count=6378*5562)
    f.close()
    # as masked array
    sst=ssttemp.reshape((5562,6378))*0.1
    sst=np.ma.masked_array(sst[::-1,:],mask=(sst[::-1,:]==0  ))

    # grid
    lon=100.+np.arange(6378)*0.01097869
    lat=9.97971+np.arange(5562)*0.00899322
    # selection
    xrange=np.logical_and(lon>=lon1,lon<=lon2)
    yrange=np.logical_and(lat>=lat1,lat<=lat2)
    sst=(sst[yrange,:])[:,xrange]-273.15 # Kelvin to C

    return lon[xrange],lat[yrange], sst
############################################
if __name__ == '__main__':
    from mpl_toolkits.basemap import Basemap
    # example plot

    #read data
    filename="n1810080316.sst.gi"
    lon1=130.
    lon2=140.
    lat1=30.
    lat2=36.
    lon,lat,sst=read_CEReS_SST(".",filename,lon1,lon2,lat1,lat2)

    #plot
    tmax=30.
    tmin=20.
    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("CEReS_SST.png")
    plt.show()