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

東北大学外洋域新世代海面水温（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に関する記事