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
を求める。
と比較。
