2010年9月1日水曜日

MODIS aqua Global Level 3 Mapped Mid-IR SST (4km)



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 件のコメント: