MODIS aqua Global Level 3 Mapped Mid-IR 海水面温度をPyNIOで読み込み、PyNGLでプロットした(daily, 4km)。
データはNASA PO.DAACのサイトのFTPサイトより、8月3日のデータを入手し(A2010215.L3m_DAY_SST4_4.bz2)、bunzip2で解凍しておく。
データはhdf4フォーマットである。
以下がスクリプト。
SSTは変数名"l3m_data"でunsigned 16bit integer
クオリティフラッグが変数名"l3m__qual"でunsigned 8bit integer
これらはHDFVIEW で確かめた。
クオリティフラッグが0以外は欠損値とした。
データは、南北が逆になっているので反転してある。
import Nio
import Ngl
import numpy as np
import scikits.timeseries as ts
#day
year=2010
month=8
day=3
t=ts.Date("D",year=year,month=month,day=day)
days='%03d' % t.day_of_year
td=str(year)+"/"+str(month)+"/"+str(day)
# open file
nfile=Nio.open_file("A%(year)s%(days)s.L3m_DAY_SST4_4" %locals(), mode='r',format="h4")
nfile.set_option("MaskedArrayMode","MaskedNever")
#grid
lon=nfile.Westernmost_Longitude+nfile.Longitude_Step[0]*np.arange(nfile.Number_of_Columns)
lat=nfile.Southernmost_Latitude+nfile.Latitude_Step[0]*np.arange(nfile.Number_of_Lines)
#grid salection
xrange=np.logical_and(lon>=130,lon<=140.)
yrange=np.logical_and(lat>=30,lat<=36.)
temp=nfile.variables["l3m_data"]
flag=nfile.variables["l3m_qual"]
temp2=((temp[::-1,:])[yrange,:])[:,xrange].astype("uint16")
flag2=((flag[::-1,:])[yrange,:])[:,xrange].astype("uint8")
sst=temp2*temp.Slope[0]+temp.Intercept[0]
sst=np.ma.masked_array(sst,mask=flag2>0)
nfile.close()
# plot by PyNGL
# Open a workstation.
#
wks_type = "png"
wks = Ngl.open_wks(wks_type,"MODIS")
resources = Ngl.Resources()
resources.tiXAxisString = "~F25~longitude"
resources.tiYAxisString = "~F25~latitude"
resources.cnFillOn = True # Turn on contour fill.
resources.cnLineLabelsOn = False # Turn off line labels.
resources.cnInfoLabelOn = False # Turn off info label.
resources.mpGridAndLimbOn = False
resources.sfXCStartV = float(min(lon[xrange])) # Define where contour plot
resources.sfXCEndV = float(max(lon[xrange])) # should lie on the map plot.
resources.sfYCStartV = float(min(lat[yrange]))
resources.sfYCEndV = float(max(lat[yrange]))
resources.mpLimitMode = "LatLon" # Limit the map view.
resources.mpMinLonF = float(min(lon[xrange]))
resources.mpMaxLonF = float(max(lon[xrange]))
resources.mpMinLatF = float(min(lat[yrange]))
resources.mpMaxLatF = float(max(lat[yrange]))
resources.mpDataBaseVersion = "MediumRes"
resources.tiMainString = "MODIS aqua Global Level 3 Mapped Mid-IR SST day=%(td)s" %locals()# Set a title.
resources.tiMainFontHeightF=0.015
resources.cnLevelSelectionMode = "ExplicitLevels" # Define own levels.
resources.cnLevels = np.arange(20.,30.,0.5)
#
# draw contours over map.
#
map = Ngl.contour_map(wks,sst,resources)
Ngl.end()
MODISにはAquaとTerra、それぞれに海面水温データとしてはMid-IR SSTとthermal IR SSTがある。どのように使いわけたらいいのだろうか?
参照
JAXA MODISページ
http://www.eorc.jaxa.jp/hatoyama/satellite/sendata/modis_j.html
EORC
http://kuroshio.eorc.jaxa.jp/ADEOS/mod_nrt_new/index.html
MODIS wikipedia
http://ja.wikipedia.org/wiki/MODIS

0 件のコメント:
コメントを投稿