東北大学 外洋域新世代海面水温

東北大学外洋域新世代海面水温（http://www.ocean.caos.tohoku.ac.jp/~merge/sstbinary/actvalbm.cgi)（分解能1/20度）を読むモジュールである。
ファイル名としてはgzipされたファイルを指定でき、解凍する必要はない。
(ftpのURLも直接指定できるがネットワークの負荷を考えること。）
メインプログラムとして走らせた時は領域一部をmatplotlibでプロットする。

readNGSST.py

import numpy as np
import matplotlib.pyplot as plt

def readNGSST(dir,file):
    """
   read Tohoku univ. New Generation Sea Surface Temp.

   usage:
      lon,lat,sst=readNGSST(dir,file)
     
   input
      dir working dir
      file filename 
           gzipped file is fine.
           URL also works
   output
      lon longitude
      lat latitude
      SST sst data (masked array)          
    """
    # read data
    ds=np.DataSource(dir)
    f=ds.open(file, 'r')
    f.seek(200) # skip header
    sst=np.frombuffer(f.read(),dtype="u1")
    f.close()

    # as masked array
    sst=np.ma.masked_array(sst,mask=np.logical_or(sst==255, sst==0))
    sst=sst.reshape(1000,1000)*0.15-3.

    # grid
    lon=116.+np.arange(1000)*0.05
    lat=13.+np.arange(1000)*0.05

    return lon,lat, sst

############################################
if __name__ == '__main__':

    # example plot

    #read data
    filename="msst2009010112.raw.gz"
    lon,lat,sst=readNGSST(dir=".",file=filename)

    #grid salection
    xrange=np.logical_and(lon>=130,lon<=145.)
    yrange=np.logical_and(lat>=30,lat<=45.)
    # plot
    plt.figure()
    plt.pcolor(lon[xrange],lat[yrange],(sst[yrange,:])[:,xrange]) 
    plt.savefig("ngsst.png")
    plt.show() 

参考
書籍
Matplotlib for Python Developers: Build Remarkable Publication Quality Plots the Easy Way

Rでwavelet

前回までのようにRにNINO3のデータを渡し、今回は連続wavelet解析を行ってみる。

Wavelet解析にはwmtsaパッケージを用いる。

以下がプログラムである。
プロットする時、series=TRUEでもとの時系列もプロットする。

from read_NINO import read_NINO

#for R
import rpy2.robjects as robjects

data=read_NINO()
x=data["NINO3 ANOM"]

year=str(x.year[0])
month=str(x.month[0])

#data to R
nino3=robjects.FloatVector(x)
robjects.globalEnv["nino3"] = nino3

#r script
r=robjects.r
rscripts="""
nino3ts=ts(nino3,start=c(%(year)s,%(month)s),frequency=12)
png('NINO3_R_wavelet.png')

#import library wmtsa
library("wmtsa")
## calculate the CWT of the sunspots series using
## a Mexican hat wavelet (gaussian2)
a<-wavCWT(nino3ts)    
## print the result
print(a)
## plot an image of the modulus of the CWT and the
## time series
plot(a,series=TRUE) 

dev.off()
"""
r(rscripts %locals())

結果

Continuous Wavelet Transform of nino3ts
---------------------------------------
Wavelet           : Mexican Hat (Gaussian, second derivative)
Wavelet variance  : 1
Length of series  : 724
Sampling interval : 0.08333333
Number of scales  : 73
Range of scales   : 0.0833333333333333 to 60.3333333333333

x軸は年。Y軸は0が1年（2^0)、２が4年(2^2)に対応する。

参考
A Practical Guide to Wavelet Analysis
http://paos.colorado.edu/research/wavelets/


書籍
Wavelet Methods for Time Series Analysis (Cambridge Series in Statistical and Probabilistic Mathematics)

Rで自己相関

前回に続いてRにNINO3のデータを渡し、今回は自己相関を求めてみる。

Rで自己相関を求める関数はacf である。
自己相関を求めるlagは７２点（７２ヶ月）までとした（lag.max=72)
自己相関の他に
自己共分散　type="covariance"
偏自己相関　type="partial"
も求めることができる。
以下、プログラム

#coding: utf-8
from read_NINO import read_NINO

#for R
import rpy2.robjects as robjects

data=read_NINO()
x=data["NINO3 ANOM"]

year=str(x.year[0])
month=str(x.month[0])

#data to R
nino3=robjects.FloatVector(x)
robjects.globalEnv["nino3"] = nino3

#r script
r=robjects.r
rscripts="""
nino3ts=ts(nino3,start=c(%(year)s,%(month)s),frequency=12)
png('NINO3_R_acf.png')

op<- par(mfrow=c(1,3)) #  描画設定：1行3列画面表示
acf(nino3ts,lag.max=72,xlab="lag (year)")              # 自己相関プロット
acf(nino3ts,type="covariance",lag.max=72,xlab="lag (year)") #自己共分散プロット
acf(nino3ts,type="partial",lag.max=72,xlab="lag (year)")    #偏相関プロット
par(op)                  #作業前の描画設定に戻す

dev.off()
"""
r(rscripts %locals())

図の横の点線は95%の信頼区間をしめす（これより絶対値の小さい相関は95%の信頼区間で無相関ではないとは言えない）。　これを90%などにしたいときはci=0.9などとする。

参考
R の基本パッケージ base, stats 中の時系列オブジェクトの簡易解説
(RjpWiki)

 R言語による時系列分析
(hamadakoichi blog)

R による時系列分析入門
（書籍）

過去のRに関する記事

Rで時系列

以前につくったpythonで取り込むNINO監視領域のデータを、統計ソフトRに持って行き、時系列データとする。Rに取り込むとRを使って統計解析ができるから便利だ。

Pythonの世界からRの世界に持っていくにはrpy2を用いる。
NINO3データを渡して、Rではtsを使って時系列オブジェクトを作り、
plot.tsでプロットする。

以下、プログラム。

from read_NINO import read_NINO


#for R
import rpy2.robjects as robjects


data=read_NINO()
x=data["NINO3 ANOM"]

year=str(x.year[0])
month=str(x.month[0])

#data to R
nino3=robjects.FloatVector(x)
robjects.globalEnv["nino3"] = nino3

#r script
r=robjects.r
rscripts="""
png('NINO3_R.png')
nino3ts=ts(nino3,start=c(%(year)s,%(month)s),frequency=12)
plot.ts(nino3ts)
dev.off()
"""
r(rscripts %locals())

参考
R の基本パッケージ base, stats 中の時系列オブジェクトの簡易解説
(RjpWiki)

 R言語による時系列分析
(hamadakoichi blog)

Rによる時系列分析入門
（書籍）

過去のRに関する記事

k-平均法

k-平均法は非階層型クラスタリング手法の一つである。

参考 
クラスタリングの定番アルゴリズム「K-means法」をビジュアライズしてみた（てっく煮ブログ）

以前つくった串本浦神の潮位差と黒潮緯度の分布（下図）を２つのグループに分けてみる。

道具はscipyのk-means2関数を用いる。分けられたグループで色分けする。大きな●はそれぞれのグループの重心である(kmeans2で出てくる重心はNormarizeされたデータのものなのでそのままはつかわないこと）。重心が収束しているかの確認もおこなっている。

#read module
from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

import matplotlib.pyplot as plt
import numpy as np

from scipy.cluster.vq import  kmeans2, whiten,vq # functions for k-means

#read data
kuroshio_lat=read_kuroshio_lat()
kushimoto_uragami=read_kushimoto_uragami()   

#kmeans
nomissing=np.logical_and(~kuroshio_lat.mask,~kushimoto_uragami.mask) #missing value should 
be excluded.
lat=kuroshio_lat.data[nomissing]
level=kushimoto_uragami.data[nomissing]
data=np.column_stack((level,lat))

whiteneddata=whiten(data) #whitened beforehand
centroids,index=kmeans2(whiteneddata,2,iter=20) #k-means

#check centroids are converged
label = vq(whiteneddata, centroids)[0]
for i in range(2):
    print 'group=',i
    print 'estimated=', centroids[i,:]
    print 'actual   =',whiteneddata[label==i,:].mean(axis=0)

#### plot
plt.figure()
#1stgroup 
plt.plot(level[index==0],lat[index==0],'r.')
plt.plot(level[index==1].mean(),lat[index==1].mean(),'go',markersize=12)
#2nd group
plt.plot(level[index==1],lat[index==1],'g.')
plt.plot(level[index==0].mean(),lat[index==0].mean(),'ro',markersize=12)

plt.xlabel("Kushimoto-Uragami")
plt.ylabel("Kuroshio Latitude")
plt.savefig("kurosho_kmeans.png")
plt.show()

出力

group= 0 estimated= [ 1.58989573 42.05977249]      actual = [ 1.58989573 42.05977249] group= 1 estimated= [ 0.14166519 40.3829155 ]      actual = [ 0.14166519 40.3829155 ]

NINO3のスペクトル密度を信頼区間つきで

前回つくった信頼区間つきでスペクトル密度をもとめる関数を、NINO3監視指数に適用する。NINO3の読み込みは以前に作ったread_NINOモジュールを用いる。

下の例ではtaperにはHammingフィルタ、周波数空間でのスムージングに1-2-1フィルタ、信頼度は90%を用いた。
黒線がスペクトル密度、赤点線が信頼区間を示す。縦の黒点線は2年と7年に対応する補助線である。

from read_NINO import read_NINO

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

from  psd_withError import psd_withError #import psd_withError module  

data=read_NINO()
x=data["NINO3 ANOM"]

# case1
Pxx,f,upper,lower=psd_withError(x, NFFT=len(x),Fs=12,window=np.hamming(len(x)),smooth＝[1,2,1],Confidencelevel=0.90)  
plt.figure()
plt.loglog(f,Pxx,'k-o')
plt.loglog(f,upper,'r--')
plt.loglog(f,lower,'r--')
plt.loglog([1./2.,1./2.],[1.E-5,1.E2],'k--') #vertical line for 2 years
plt.loglog([1./7.,1./7.],[1.E-5,1.E2],'k--') #vertical line for 7 years
plt.savefig("NINO_psd_withError.png")
plt.show() 

psd関数をスペルトル密度誤差も出せるよう拡張

前回の見延さんのスペクトル密度の推定誤差理論を使って、以前用いたmatlabplotlib.mlab.psd関数を拡張したモジュールを作った(psd_withError)。

stand aloneで用いる時(if __name__ == '__main__':以下)は、見延さんのスクリプトと同じことを確認するプログラム。結果が同じにならないのは、

• 見延さんのスクリプトではfftの後データ数で割ってから二乗をとっている
• matplotlibのpsd関数でonesidedを取るとき2倍している
• matplotlibのpsd関数では、データテーパーをかけることによる振幅の減少を補償している

ためである。もし同じ結果を得たければ、プログラム中のコメントアウトしている部分をはずすことで可能である。

 
psd_withError.py

import numpy as np
from numpy.random import randn
import scipy as sp
import scipy.signal as signal
import scipy.stats as stats

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


#------------------------------------------
def psd_withError(x, NFFT=256, Fs=2, detrend=mlab.detrend_none, window=mlab.window_hanning,noverlap=0, pad_to=None, sides='def
ault', scale_by_freq=None,smooth=None,Confidencelevel=0.9):
    """
Extention of matplotlib.mlab.psd
The power spectral density with upper and lower limits within confidence level
You should not use Weltch's method here, instead you can use smoother in frequency domain with *smooth* 

Same input as matplotlib.mlab.psd except the following new inputs
*smooth*:
    smoothing window in frequency domain  
    for example, 1-2-1 filter is [1,2,1]
    default is None
*Confidencelevel*
    Confidence level to estimate upper and lower estimates.
    Default is 0.9

Returns the tuple (*Pxx*, *freqs*,*upper*,*lower*).
    upper and lower are limits of psd within confidence level.   
"""

    Pxxtemp,freqs = mlab.psd( x, NFFT, Fs, detrend, window, noverlap, pad_to, sides,scale_by_freq)

# commented out, if you want to get Minobe-san's result  Pxxtemp=Pxxtemp/float(NFFT)/2.*((np.abs(window)**2).mean())

#smoothing
    smooth=np.asarray(smooth)
    if smooth is not None:
        avfnc=smooth/float(sum(smooth))
        Pxx=np.convolve(Pxxtemp[:,0],avfnc,mode="same")
    else:
        Pxx=Pxxtemp[:,0]
        avfnc=np.asarray([1.])

#estimate uppper and lower estimate woth equivalent degree of freedom
    if pad_to is None:
        pad_to = NFFT

    if cbook.iterable(window):
        assert(len(window) == NFFT)
        windowVals = window
    else:
        windowVals = window(np.ones((NFFT,), x.dtype))

    edof=(1.+(1./sum(avfnc**2)-1.)*float(NFFT)/float(pad_to)*np.mean(windowVals))# equivalent degree of freedom
    a1=(1.-Confidencelevel)/2.
    a2=1.-a1
    lower=Pxx*stats.chi2.ppf(a1,2*edof)/stats.chi2.ppf(0.50,2*edof)
    upper=Pxx*stats.chi2.ppf(a2,2*edof)/stats.chi2.ppf(0.50,2*edof)
    return Pxx,freqs,upper,lower

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

if __name__ == '__main__':
    plt.figure()
    for npnl in range(4):
        tlng=100; nsmpl=1000
        if npnl==0 :
            nfft=100
            taper=signal.boxcar(tlng)
            avfnc=[1, 1, 1, 1, 1];
            cttl='no padding, boxcar, 5-running'
        elif npnl==1:
            nfft=100
            taper=signal.hamming(tlng)
            avfnc=[1, 1, 1, 1, 1]
            cttl='no padding, hamming, 5-running'
        elif npnl==2:
            nfft=200
            taper=signal.boxcar(tlng)
            avfnc=[1, 1, 1, 1, 1];
            cttl='double padding, boxcar, 5-running'
        elif npnl==3:
            nfft=200
            taper=signal.hamming(tlng)
            avfnc=np.convolve([1, 2, 1],[1, 2, 1])
            cttl='double padding, hamming, double 1-2-1'
    
        tsrs=randn(tlng,nsmpl)
        ds_psd=np.zeros([nfft/2+1,nsmpl])
        upper_psd=np.zeros([nfft/2+1,nsmpl])
        lower_psd=np.zeros([nfft/2+1,nsmpl])
        for n in range(nsmpl):
            a,b,u,l=psd_withError(tsrs[:,n],NFFT=tlng,pad_to=nfft,Fs=1,window=taper,smooth=avfnc,Confidencelevel=0.9)
            ds_psd[:,n]=a[:]
            upper_psd[:,n]=u[:]
            lower_psd[:,n]=l[:]
        frq=b[:]

# 90% confidence level by Monte-Carlo
        srt_psd=np.sort(ds_psd,axis=1)
        c=np.zeros([frq.size,2])                
        c[:,0]=np.log10(srt_psd[:,nsmpl*0.05])
        c[:,1]=np.log10(srt_psd[:,nsmpl*0.95])

# estimate from extended degree of freedom 
        ce=np.zeros([frq.size,2]) 
        ce[:,0]=np.log10(np.sort(upper_psd,axis=1)[:,nsmpl*0.5])
        ce[:,1]=np.log10(np.sort(lower_psd,axis=1)[:,nsmpl*0.5])
#plot
        plt.subplot(2,2,npnl+1)
        plt.plot(frq,c,'b',frq,ce,'r')
        plt.title(cttl)
        plt.xlabel('frq')
        plt.ylabel('psd')
        if (npnl==0): plt.legend(('Monte-carlo','','Theory'),'lower center',labelspacing=0.05)
    plt.subplots_adjust(wspace=0.6,hspace=0.4)
    plt.savefig("psd_withError.png")
    plt.show()