今年8月の500hpa高度場

異常高温の夏、今後増加？　
エルニーニョや偏西風蛇行で
１８９８年の統計開始後最高だった今夏の高温の原因について、気象庁の「異常気象分析検討会」（会長、木本昌秀東大大気海洋研究所教授）が３日、見解をまとめた。春まで継続したエルニーニョと、夏に発生したラニーニャ両現象の相乗効果や、日本上空付近を吹く偏西風の蛇行などの複合要因を指摘した。(2010/09/03　47News)　web魚拓

偏西風の位置の指標としての500hpa高度場(geopotential height)をGoogle Earthでプロットしてみた。
GradsにはGoogle Earth用のデータを作るオプションがある(version 2.0a5以降）
http://www.iges.org/grads/gadoc/gradcomdsetkml.html
それを用いて、NCEP NCAR/Reanalysisの今年8月のデータから画像を作るGradsスクリプトが以下の通り。

*************KML
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/hgt.mon.mean.nc"
"set lev 500"
"set time aug2010"
"set x 1 144"
"set y 1 73"
"set kml hgt2010aug"
"set gxout kml"
"d hgt"
************FIGURE
"c"
"set parea 0 11 0 8.5"
"set grid off"
"set grads off"
"set mproj scaled"
"set mpdraw off"
"set gxout shaded"
"color 5400 5900 20 -kind blue->skyblue->yellow->red->brown"
"d hgt"
"printim  hgt2010aug.png x2048 y1024 -t 0"

前半で、
hgt2010aug.kml
hgt2010aug.tif
という二つのファイルが出来る。

本来なら、これだけで十分なはずだが、これだと図がいまいちなので、
http://gradsusr.org/pipermail/gradsusr/2010-January/010117.html
を参考に後半で図だけを好きなように作り直す。
その後、hgt2010aug.kmlをエディタで開いて、hgt2010aug.tifで指定しているところをhgt2010aug.pngに書き換える。

後は、hgt2010aug.kmlをGoogle Earthで開けばよい。

ちなみに、”set x 1 144”としているが、本来このデータはサイクリック用のグリッドが145まで用意されているが、そのままだとうまく端つながらないので、144までにしている。

また、データはOPeNDAPを使いネットから直接読み、シェイドの色をcolorスクリプト
http://kodama.fubuki.info/wiki/wiki.cgi/GrADS/script/color.gs?lang=jp
で指定した。

できたデータをブログに貼るために、
Embed your Google Earth Project in your Website
http://earth.google.com/outreach/tutorial_kmlembed.html
を利用した。
データは一旦Google Earthで開いたファイルをkmz方式で保存し、こちらにおいた。


(注.firefoxでは見れないことがあるようだ。その場合は他のブラウザでためしてみてください。）

TRMM 3B42 降水データ

人工衛星TRMMの降水データ3B42をプロットする。
この降水データは水平分解能0.25度で3時間ごとにデータがある。
このデータをNASA MiradorからOPeNDAPを通して入手する。
以下の例は7月27日0時世界標準時（日本標準時9時）の日本付近のデータをプロットする例である。この日は東シナ海に熱帯低気圧が存在し、それによる雨が見える。
OPeNDAPからはPyNIOを使ってデータを入手する。
関数name_makerはファイルのあるURLはこういう規則だろうと推測して作ったものである。間違いかもしれないので確認すること。

import Nio
import numpy as np

def read_TRMM_3B42(date,lon1,lon2,lat1,lat2):
    """
   read TRMM 3B42 3 hourly rainfall  from NASA Mirador

   usage:
      lon,lat,sst=read__TRMM_3B42(date,lon1,lon2,lat1,lat2)
     
   input
     
      date (year,month,day,hour) 
      lon1,lon2  Western,East bounday
      lat1,lat2  South,West bounday
   output
      lon longitude
      lat latitude
      rainfall data (mm/hour)         
    """
#check
    assert date[3] % 3 ==0, "Error, hour must be multiple of 3!!!"
# open OPENDAP
    fname=name_maker(date[0],date[1],date[2],date[3])
    f = Nio.open_file(fname)
## variable name check
    #variables = f.variables.keys()
    #print f
    #print "variables",variables

# read data
    temp=f.variables["DATA_GRANULE.PlanetaryGrid.precipitation"]

#grid
    # grid
    lon=-179.875+np.arange(temp.shape[1])*0.25
    lat=-49.875 +np.arange(temp.shape[2])*0.25
    # selection
    xrange=np.where(np.logical_and(lon>=lon1,lon<=lon2))
    yrange=np.where(np.logical_and(lat>=lat1,lat<=lat2))
    
#precipitation
    x=xrange[0]
    y=yrange[0]
    precip=(temp[0,x[0]:x[-1]+1,y[0]:y[-1]+1].T).reshape((len(y),len(x)))

# closing 
    f.close()
    return lon[xrange],lat[yrange],precip

############################################

def name_maker(year,month,day,hour):
    import scikits.timeseries as ts
    #
    t=ts.Date('H', year=year,month=month,day=day,hour=hour)
    if t.hour==0:
        tc=t-1
    else:
        tc=t
    cyear = "%04d"  % tc.year
    cyear2= "%02d"  % (t.year-2000)
    cmonth2= "%02d" % t.month
    cday2=   "%02d" % t.day
    shour = str(t.hour)
    yd3= "%03d" % tc.day_of_year
    #
    fname="http://disc2.nascom.nasa.gov/opendap/TRMM_L3//TRMM_3B42//%(cyear)s/%(yd3)s/3B42.%(cyear2)s%(cmonth2)s%(cday2)s.%(shour)s.6A.HDF.Z" %locals()
    return fname

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

# range
    lon1=120.
    lon2=145.
    lat1=20.
    lat2=40.

#date
    year=2009
    month=7
    day=27
    hour=0 
    lon,lat,precip=read_TRMM_3B42(date=(year,month,day,hour),lon1=lon1,lon2=lon2,lat1=lat1,lat2=lat2)

   #plot
    vmax=int(precip.max())
    vmin=int(precip.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,precip,vmin=vmin,vmax=vmax,cmap=plt.cm.Blues)
    ti=str(year)+"/"+str(month)+"/"+str(day)+":"+str(hour)+"Z"
    plt.title(ti)
# 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("TRMM_3B42.png")
    plt.show()

参考
NASA Giovanniのサイト
http://disc2.nascom.nasa.gov/Giovanni/tovas/TRMM_V6.3B42.2.shtml

TRMM 3B42 readme
ftp://meso-a.gsfc.nasa.gov/pub/trmmdocs/3B42_3B43_doc.pdf

Sageで方程式を解く

Sageで方程式を解いてみる。

x,a,b,c=var('x a b c')
solve(x^2+1==0,x)

f=a*x^2+b*x+c==0
solve(f,x)

Q=solve(f,x)
for q in Q:
    q.substitute({a:1,b:-5,c:6})

Latexの書式を出すときには

latex(_)

などとする。

方程式の解を指定範囲で全て探すプログラム

以下のプログラムは、方程式を関数f(x)=0で与えたときに、指定した範囲内で解を全て求めることを意図したものである。
基本は、範囲をdivの数で分割し、それぞれの領域でScipyのbrentqを適用し、解を集めている。
しかし、関数f(x)が解の前後で符号を変えないケースがあり、このような場合の解を求めるためにNewton法も用いている。

f1d_solve_range.py

import numpy as np

def f1d_solve_range(f,range,div=1000,args=(),newtontirigger=1.E-4):
    """
Solve f(x)=0 to otain all solutions between range
range is devided by div to find solution in each subrange

When |f(x)| < 1.E-4, Newton Method is also used
because brentq can not find solution if f(x) does not change sign at a solution.
"""
    from scipy.optimize import brentq,newton

    x=np.linspace(range[0],range[1],div)
    dx=x[1]-x[0] 

    #solution 1 by bysection
    solution1=[]
    for x0  in x[:-1]:
        try:
            a=brentq(f,x0,x0+dx,args=args)
            solution1.append(a)
        except:
            continue
    solution1=np.unique(solution1)

    # solution2 by Netwon method just in case f(a) f(b) has same sign
    solution2=[] 
    for x0 in x:
        if np.abs(f(x0,*args))<newtontirigger:
            b=newton(f,x0,args=args)
            solution2.append(b)
    solution2=np.unique(solution2)
    solution2.clip(range[0],range[1])  # remove out of range

    # final solution
    solution=solution1
    for s2 in solution2:
        if np.alltrue(np.abs(solution-s2)>dx):
            print "New solution added by Newtonian!:",s2
            solution=np.append(solution,s2)
    
    return np.sort(solution)  

例１　

#example 1
    def f(x):
        return x**3 
    
    print "ex1=",f1d_solve_range(f,(-10.,10))

結果:

ex1= [ 0.]

例2　

#example 2
    def f(x):
        return x**3-x 
    print "ex2=",f1d_solve_range(f,(-10.,10)) 

結果：
ex2= [-1.  0.  1.]


例3

#example 3
    def f(x):
        return x**3-x**2 
    print "ex3=",f1d_solve_range(f,(-10.,10)) 
    x=np.linspace(-2,2,100)

結果：
ex3= New solution added by Newtonian!: 1.92717530048e-08
[  1.92717530e-08   1.00000000e+00]
解の前後で符号を変えないケースでNewton法で解が求まっている。

のグラフを参照

例4

#example 4
    def f(x):
        return x**2+1 
    print "ex4=",f1d_solve_range(f,(-10.,10)) 

結果：

ex4= []

解なし。

例5

#example 5
    def f(x,a,b,c):
        return a*x**2 + b*x +c 
    print "ex5=",f1d_solve_range(f,(-10.,10),args=(1,2,1)) 

結果:
ex5= New solution added by Newtonian!: -0.999999982865
[-0.99999998]

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

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()

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