2010年3月30日火曜日

ENSOの周期 psd関数のオプション2 Zero padding



前回の続き。
pad_toによってzero paddingの量を変えることができる。

pad_to: integer The number of points to which the data segment is padded when performing the FFT. This can be different from NFFT, which specifies the number of data points used. While not increasing the actual resolution of the psd (the minimum distance between resolvable peaks), this can give more points in the plot, allowing for more detail. This corresponds to the n parameter in the call to fft(). The default is None, which sets pad_to equal to NFFT

下がその例である。
差を良く見るために2年から7年の周期(周波数1/2から1/7)の間だけプロットしている。
黒線がzero paddingを行わない場合。
赤線がデータの長さまで、zero paddingを行う場合。つまり、これはzero paddingを行わないのと同じだ。
青線がデータ長の2倍まで zero paddingを行う場合。より多くの点が加わっていることがわかる。

from read_NINO import read_NINO

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
import numpy as np

data=read_NINO()
x=data["NINO3 ANOM"]
Pxx,f=mlab.psd(x, NFFT=len(x),Fs=12)  
# zero padding (length of data) 
Pxx2,f2=mlab.psd(x, NFFT=len(x),pad_to=len(x),Fs=12)
# zero padding (length of data) *2
Pxx3,f3=mlab.psd(x, NFFT=len(x),pad_to=2*len(x),Fs=12)
plt.figure()
plt.loglog(f,Pxx,'k-o')
plt.loglog(f2,Pxx2,'r-+')
plt.loglog(f3,Pxx3,'b-x')
plt.xlim([1./7.,1./2.])
plt.ylim([10E-3,10e1])
plt.savefig("NINO_psd_zeropadding.png")
plt.show()

2010年3月29日月曜日

ENSOの周期 psd関数のオプション1 Taper



前回使ったmatplotlib.mlabのpsd関数にはいくつかオプションがある。それを試してみよう。
デフォルトではもとのデータにHanningウィンドウ(データ・ウィンドウ、データ・テーパ)をかけている。
これはオプションで変えることができる。

ウィンドウをはずす場合
前回のコードに続けて。
window=mlab.window_noneがウィンドウ無しの場合である。
#window_none()
Pxx_no,f_no=mlab.psd(x, NFFT=len(x),Fs=12,window=mlab.window_none) 
plt.figure()
plt.loglog(f,Pxx,'k')
plt.loglog(f_no,Pxx_no,'r')
plt.loglog([1./7.,1./7.],[1.E-9,1.E2],'k--')
plt.loglog([1./2.,1./2.],[1.E-9,1.E2],'k--')
plt.savefig("NINO_psd_no.png")
plt.show()
データの端の影響でスペクトルの漏れがあり、フィルタ無しの場合(赤線)にはピークの立つ場所が増える。

次は、Hammingウィンドウをかける場合。
window=np.hamming(len(x))でデータと同じ長さのHammingウィンドウをかけている。
#numpy hamming
Pxx_ham,f_ham=mlab.psd(x, NFFT=len(x),Fs=12,window=np.hamming(len(x))) 
plt.figure()
plt.loglog(f,Pxx,'k')
plt.loglog(f_ham,Pxx_ham,'r')
plt.loglog([1./7.,1./7.],[1.E-9,1.E2],'k--')
plt.loglog([1./2.,1./2.],[1.E-9,1.E2],'k--')
plt.savefig("NINO_psd_ham.png")
plt.show()
このケースではHanningフィルター(黒線)とHammingフィルター(赤線)とではあまり差がないが、一般にはHammingフィルターの方が良いようだ(見延2001)。

参考文献
見延 2001:大気海洋物理学特論4 大気海洋統計データ解析 第3章


2010年3月28日日曜日

Matplotlib for Python Developers



このブログでよく使っているmatplotlibの本が出ているので、紹介。
こちらに書評が出ている。

2010年3月12日金曜日

ENSOの周期



ENSO(エルニーニョ)の周期を見てみる。

データは前回作ったモジュールでNINO3領域ののアノマリを得る。

スペクトル密度を求める関数がmatplotlib.mlab(Matlab風の関数群)にpsdとしてある。
(図まで描いてくれるmatplotlib.pyplotのpsdというのもある。)

そのpsdを用いて作ったスペクトル密度の図とプログラムが以下の通り。
データ数(NFFT)はすべてを使う。
時間の単位は年とする。毎月のデータがあるから単位時間あたりのサンプル数(Fs)は12。
エルニーニョは2年から7年の周期があると言われているから補助線として2年(周波数1/2)と7年(周波数1/7)に当たるところに補助線(点線)を引いている。実際、それらの間にスペクトルピークが見られている。
from read_NINO import read_NINO

import matplotlib.mlab as mlab
import matplotlib.pyplot as plt

data=read_NINO()
x=data["NINO3 ANOM"]
Pxx,f=mlab.psd(x, NFFT=len(x),Fs=12)  
plt.figure()
plt.loglog(f,Pxx)
plt.loglog([1./7.,1./7.],[1.E-9,1.E2],'k--')
plt.loglog([1./2.,1./2.],[1.E-9,1.E2],'k--')
plt.savefig("NINO_psd.png")
plt.show()

2010年3月6日土曜日

エルニーニョ: NINO監視領域の温度偏差



寒暖の差大きく、大雪も 12~2月、気象庁まとめ (日経新聞)
気象庁は1日、この冬(昨年12月~今年2月)の天候まとめを発表した。3カ月間の平均気温は全国的に高めとなったものの、暖冬をもたらす「エルニーニョ現象」と北極域の気圧が高くなる「北極振動」の影響で、例年以上に寒暖の差が大きくなった。  気象庁によると、今冬はエルニーニョ現象の影響で西高東低の冬型の気圧配置が長続きせず、平均気温は西日本で平年より1.0度、東日本で0.9度、北日本と沖縄・奄美は0.6度高くなった。  季節外れの暖かさとなったのは1月下旬と2月下旬で、2月25日には全国約160カ所で、2月の最高気温の記録を更新した。(01日 19:53)
エルニーニョを話題とする。エルニーニョを監視するために、東太平洋のエルニーニョ監視領域の海面温度がよく使われる。エルニーニョ監視領域
エルニーニョ監視領域の温度とその平年値からの偏差(anomaly)のデータはNOAAのclimate predictionセンターから手に入るここの情報によれば、平年値は1971年から2000年までの平均で定義される。

以下は、これらのエルニーニョ監視領域の絶対値と偏差をネットから取得するモジュールと、4つの監視領域の偏差をデータ全期間もしくは過去48ヶ月分プロットするプログラムである。今回は正で赤、負で青で塗り分けるために、matplotlibのfill_betweenを使った。

NINO_all.png
 
NINO_48months.png
read_NINO.py
import numpy as np
import scikits.timeseries as ts

def read_NINO():
    """
    output output NINO index data
    data souce is http://www.cpc.ncep.noaa.gov/data/indices/sstoi.indices
    Explanation is http://www.cpc.ncep.noaa.gov/data/indices/
    usage:  from read_NINO import read_NINO 
            data=read_NINO()
    outputs are     
    data["NINO1+2"]     : NINO1+2  
    data["NINO1+2 ANOM"]: NINO1+2 ANOMARY  
    data["NINO3"]       : NINO3
    data["NINO3 ANOM"]  : NINO3 ANOMALY
    data["NINO4"]       : NINO4
    data["NINO4 ANOM"]  : NINO4 ANOMALY
    data["NINO3.4"]     : NINO3.4
    data["NINO3.4 ANOM"]: NINO3.4 ANOMALY  
    """  
    ninodata={}
    names=('YEAR','MONTH','NINO1+2','NINO1+2 ANOM','NINO3','NINO3 ANOM','NINO4','NINO4 ANOM','NINO3.4','NINO3.4 ANOM')
    ds = np.DataSource(None)
    f=ds.open('http://www.cpc.ncep.noaa.gov/data/indices/sstoi.indices')
    data=np.loadtxt(f,skiprows=1,\
         dtype={'formats':('i4','i4','f4','f4','f4','f4','f4','f4','f4','f4'),\
                'names': names})
    firstday=ts.Date(freq='M',year=data['YEAR'][0],month=data['MONTH'][0])
    for n in names[2:]:
          ninodata[n]=ts.time_series(data[n],start_date=firstday,freq='M')
    return ninodata

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

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import matplotlib.dates as mdates

    data=read_NINO()

#plot
    for nt in range(2):
        temp=data['NINO3'].dates # dummy to get dates of time series

        if nt==0:   # all data
            t=ts.date_array(start_date=temp[0],end_date=temp[-1])
            figname="NINO_all.png"
        else:       # last 48 months
            t=ts.date_array(start_date=temp[-1]-47,end_date=temp[-1])
            figname="NINO_48months.png"

        fig=plt.figure() 
        names=('NINO1+2 ANOM','NINO3 ANOM','NINO4 ANOM','NINO3.4 ANOM')
        for n,name in enumerate(names):
            x=(data[name])[t].dates.tolist()
            y=(data[name])[t].series
            ax=fig.add_subplot(len(names),1,n+1)
            ax.plot(x,y,'k-')
            ax.fill_between(x,y1=y,y2=0,facecolor='red',where=y>=0)  #positive red
            ax.fill_between(x,y1=y,y2=0,facecolor='blue',where=y<=0)  #positive red
            if (nt==0): # axis for all data 
                years    = mdates.YearLocator(5) #major tick every 5 years 
                ax.xaxis.set_major_locator(years)
            else:       # axis for latst 48 months
                years    = mdates.YearLocator(1) #major tick every 1 year
                ax.xaxis.set_major_locator(years)
                months    = mdates.MonthLocator() # minor tick every month
                ax.xaxis.set_minor_locator(months) 
            plt.title(name)                  #tile
        plt.subplots_adjust(hspace=0.5)             #adjust subplot scape 
        plt.savefig(figname)
        plt.show()

2010年3月4日木曜日

串本と浦神の潮位差 (日平均データ)



以前、串本と浦神の潮位差の月平均データを読んだが、今回は日データを読む。
データは気象庁のサイトで手に入れることができる(速報値)。
http://www.data.kishou.go.jp/kaiyou/db/kobe/kuroshio/chouisa/chouisa.php

以下は、データを読むモジュールと時系列プロットするプログラム。

2010/4/6 追記 データのフォーマットと値に修正があり、それに伴いプログラムと図を変更してある。

 

import numpy as np
import scikits.timeseries as ts

def read_kushimoto_uragami_daily():
    day=[]
    value=[]
    ds = np.DataSource(None)
    f=ds.open('http://www.data.kishou.go.jp/kaiyou/db/kobe/kuroshio/chouisa/chou
isa.dat')
    l=0
    for line in f:
         l+=1
         if (l$gt;1): # skip one line
             fields =line.strip().split()
             if len(fields)$gt;=1: day.append(fields[0])
             if len(fields)==1:
                 value.append(-999.)
             elif len(fields)==2:
                 try:
                     a=float(fields[1]) # check if it is number
                     value.append(a)
                 except:
                     value.append(-999.) # if data is not number
    f.close()
    valuea=np.asarray(value,dtype='float32')
    kushimoto_uragami=ts.time_series(valuea,start_date=day[0],freq='d',mask=(val
uea<-990.))    
    return kushimoto_uragami

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import scikits.timeseries.lib.plotlib as tpl

    kushimoto_uragami=read_kushimoto_uragami_daily()

    fig = tpl.tsfigure()
    fsp = fig.add_tsplot(111)
    fsp.tsplot(kushimoto_uragami, '-')
    plt.xlabel("YEAR")
    plt.ylabel("Kushimoto-Uragami (cm)")
    plt.savefig("kushimoto_uragami_daily.png")
    plt.show()