Ocean Science Hack

GrADSで作成するpostscriptの不要な横線

2012-02-06T16:48:00.001+09:00

GrADSでpostscript形式の図を作ると不要な横線が入ってしまうことがある。

従来の対策は、例えば
http://www.sci.hokudai.ac.jp/~sasakiyo/tips_grads.html

GrADS 2.0.0から従来のgxout shaded に加えて、gxout shade2が加わっている。
http://www.iges.org/grads/gadoc/gradcomdsetgxout.html
これで図を作ったところ、従来横線が入っていた場合に、横線が入らなくなった。常にそうであるかは要確認。

フィルタの応答関数

2011-10-17T08:51:00.002+09:00

フィルタの周波数応答関数をプロットしてみる。

以下の例は、京都大学の講義「MATLABで時系列解析」の線形システムの伝達関数を求めるを翻案したものである。

想定としては、1ヶ月ずつのデータがあり、シグナルとしては、48ヶ月周期と6ヶ月周期とホワイトノイズが重なったシグナルに対して13ヶ月の移動平均をかけている。

理論的な周波数応答関数を求める関数はscipyのfreqz関数である(matlabのfreqz関数に対応する)。応答関数は入力 x のパワースペクトル密度関数 Pxx と出力のクロススペクトル密度関数 Pyx の比としても推定できる。スペクトル密度関数を求める関数psdに関しては何度か紹介した。

#!/usr/bin/env python

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as plab
import matplotlib.mlab as mlab
from scipy.signal import freqz,lfilter

#--- Signal

Fs = 1   # frequency 
T = 1000.    # final time
t = np.arange(0.,T+1./Fs,1./Fs)  # time
x = np.sin(2.*np.pi*t/6.) + np.sin(2.*np.pi*t/48.)+np.random.normal(0,0.1,len(t)) # input signal

#---- power spectral by Welch's method 

plt.figure()
nfft = 256
window = np.hanning(256) 
noverlap = 128 
[Pxx,f] = mlab.psd(x,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # Power spectral
plt.subplot(2,2,1) 
plt.plot(f,10*np.log10(Pxx)) # plot in decibel 
plt.legend(('Pxx',))
plt.plot([1./12.,1./12],[-70.,20.],'k--')
plt.plot([1./6.,1./6],[-70.,20.],'k--')
plt.plot([1./48.,1./48],[-70.,20.],'k--')
plt.ylim((-70,20))

#--- Moving average

h = np.ones(13)/13.

#--- Calculate theoritical response

fs,Res = freqz(h)
plt.subplot(2,2,2)
plt.plot(fs/2./np.pi,np.abs(Res))
plt.legend(('theoretical response',))
plt.plot([1./12.,1./12],[0.,1.],'k--')
plt.plot([1./6.,1./6],[0.,1.],'k--')
plt.plot([1./48.,1./48],[0.,1.],'k--')
plt.ylim((0,1.))


#--- Apply Filter

y = lfilter(h,1,x) # filter 

[Pxy,f] = mlab.csd(x,y,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # cross scectral density

[Pyy,f] = mlab.psd(y,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none)  # power spectral density 

plt.subplot(2,2,3)
plt.plot(f,10.*np.log10(abs(Pyy)))
plt.legend(('Pyy',))
plt.plot([1./12.,1./12],[-70.,20.],'k--')
plt.plot([1./6.,1./6],[-70.,20.],'k--')
plt.plot([1./48.,1./48],[-70.,20.],'k--')
plt.ylim((-70,20))


#--- Calculate estimated response

Resc = Pxy/Pxx 
plt.subplot(2,2,4)
plt.plot(f,abs(Resc)); 
plt.legend(('estimated response',)) 
plt.plot([1./12.,1./12],[0.,1.],'k--')
plt.plot([1./6.,1./6],[0.,1.],'k--')
plt.plot([1./48.,1./48],[0.,1.],'k--')
plt.ylim([0,1])

#--- Save fig
plt.savefig("test_response.png")
plt.show()

以下がフィルタを12ヶ月移動平均にした場合。

以下はシグナルからホワイトノイズを除いた場合。この時はスペクトルからの応答関数の推定はうまくいかない（移動平均はこの場合は13ヶ月）。

参考）　Scipy Cookbook FIRFilter

気象庁全球数値予報モデルGPV (GSM)

2011-05-06T18:27:00.000+09:00

気象庁全球数値予報モデルGPV (GSM) の初期値データを京都大学のサイトからダウンロードして、清水慎吾さん作成のgsm2binでgrib2形式からgrads形式に変換するモジュール。

以下がそのモジュール。サンプルのメインプログラムでは、２０１１年１月１日6時のデータを変換し、python interface to gradsで海面気圧をプロットしている。

get_gsm.py

def get_gsm(date=[2011,1,1,0],dir="."):
    """
   download jma gsm initial data from   http://database.rish.kyoto-u.ac.jp/arch/jmadata/da
ta/gpv/original/
   and
   convert to grads file using gsm2bin by Dr. Shingo Shimizu
   (http://shimizus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%E2%A5%C7%A5%EB%B4%D8%CF%A2%A5%E1%A
5%E2)

   usage:
      get_gsm(date,dir)
     
   input
     
      date: [year,month,day,hour]  hour must be multiple of 6
            default [2011,1,1,0]
      dir : directory to save data    
    """
    import os,urllib

#check
    assert date[3] % 6 ==0, "Error, hour must be multiple of 6!!!"

#names from date
    urldir,original,gradsbinname,gradsdate=name_maker(date)

# get original data
    if not os.path.exists(dir+"/"+original):
        if not os.path.exists(dir): os.mkdir(dir) 
        # Download the data
        print 'Downloading data, please wait'
        opener = urllib.urlopen(urldir+original)
        open(dir+"/"+original, 'w').write(opener.read())
        print 'Downloaded'

# convert to grads file
    os.system("gsm2bin  %(dir)s/%(original)s %(dir)s/%(gradsbinname)s" %locals())

# make ctlfile
    make_ctl(dir,gradsbinname,gradsdate)

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

def name_maker(date):
    """
    make filenames from date
"""
    
    import datetime
    time=datetime.datetime(*date)
    cyear=time.strftime("%Y")
    cmonth=time.strftime("%m")
    cmonth2=(time.strftime("%b")).lower()
    cday=time.strftime("%d")
    chour=time.strftime("%H")

    # URL of original data
    urldir="http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/%(cyear)s/%(
cmonth)s/%(cday)s/" %locals()

    # File name of original data
    original="Z__C_RJTD_%(cyear)s%(cmonth)s%(cday)s%(chour)s0000_GSM_GPV_Rgl_FD0000_grib2.
bin" %locals()

    # grads bin name after conversion
    gradsbinname="gsm%(cyear)s%(cmonth)s%(cday)s%(chour)s.bin" %locals()

    # date in grads format
    gradsdate="%(chour)s:00Z%(cday)s%(cmonth2)s%(cyear)s"  %locals()
    return urldir,original,gradsbinname,gradsdate

###########################################
def make_ctl(dir,gradsbinname,gradsdate):
    """
    make grads ctl file for GSM data 
"""
    #ctl file text

    ctltxt="""
dset ^%(gradsbinname)s
options template
undef -999
xdef 720 LINEAR 0.0 0.5
ydef 361 LINEAR -90.0 0.5
zdef 17  LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70 50
30 20 10
tdef 1 LINEAR %(gradsdate)s 6hr

vars 16
slp 0  0 sea level pressure
sp  0  0 surface pressure
su  0  0 surface westerly  wind
sv  0  0 surface southerly wind
st  0  0 surface temperature
srh 0  0 surface relative humidity
lca 0  0 LowCloudAmount
mca 0  0 MidCloudAmount
hca 0  0 HighCloudAmount
tca 0  0 TotalCloudAmount
z   17 0 height  
u   17 0 westerly wind
v   17 0 southerly wind
t   17 0 temperature
rh  17 0 relative humidity
w   17 0 p-velocity
endvars
"""

    # write to file
    gradsctlname=dir+"/"+gradsbinname.replace(".bin",".ctl")
    f=open(gradsctlname,"w")
    f.write(ctltxt %locals())
    f.close()

#-----------------------------------------------------------------------

if __name__ == '__main__':

    # get gsm data
    get_gsm(date=[2011,1,1,6],dir="DATA")

    # plot
    from grads.galab import GaLab  # for python interface for grads
    ga = GaLab(Bin='grads',Window=False,Echo=True,Verb=0,Port=False)
    ga.open("DATA/gsm2011010106.ctl")
    script="""
    set grads off
    set gxout shaded
    d slp/100
    cbarn
    draw title Sea level pressure (hpa)
    printim gsm.png white
    """
    ga(script)

気象庁MGDSST

2011-03-05T09:52:00.000+09:00

気象庁の海面水温データ（0.25x0.25度、daily)であるMGDSSTをウェブからダンロードしてから、データを読み込むまでのモジュールをつくる。MGDSSTはNEAR-GOOSのサイトから入手が可能である。
idとpasswordが必要なので、あらかじめ登録しておく。

以下がスクリプト
usrname="xxxxx"
password="*****"
のところを自分のidとpasswordで置き換える。

import numpy as np 
import os

def downloadfile(infile="remotefile",outfile="localfile"):
    """
    download file from NEAR-GOOS data site with id and password
    input
       infile: remote file on the web
       outfile: local file
    
    set password and id below
    """
    import urllib2
    host="goos.kishou.go.jp:80"
    realm="NEAR-GOOS"
    usrname="xxxxx"
    password="*****"
    auth_handler = urllib2.HTTPBasicAuthHandler() 
    auth_handler.add_password(realm, host, usrname, password) 
    opener = urllib2.build_opener(auth_handler) 
    urllib2.install_opener(opener) 
    f=urllib2.urlopen(infile)
    open(outfile, 'w').write(f.read())

#----------------------------------------------------------
def readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)):
    """
    load MGDSST data from NEAR-GOOSE site
    read data with fortran (f2py)
    select region specified by lonrange and latrange 

    Usage:
    lon,lat,sst= readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)

    Input 
       dir: directory for data download
       year,month,day:  date of the data 
       lonrange,latrange = longitude and latitude range for data
    output
       lon, lat: longitude, latitude
       sst:  sst data
"""

# date
    import datetime
    t=datetime.datetime(year,month,day)
    syear=t.strftime("%Y")
    smon=t.strftime("%b")
    sday=t.strftime("%d")
# filename
    dataname="mgdsst."+smon+sday
    filename=dir+"/"+dataname+"."+syear
    filenamegz=filename+".gz"
    print filename
    urlname="http://goos.kishou.go.jp/rrtdb/usr/pub/JMA/mgdsst/"+syear+"/"+dataname+".gz"

# file exist?
    if not os.path.exists(filename):
        if not os.path.exists(filenamegz):
           downloadfile(infile=urlname,outfile=filenamegz)
        os.system("gunzip %(filenamegz)s" %locals())         

# temp.bin is used for temporary file
    if os.path.lexists("temp.bin"):
        raise IOError("temporary file temp.bin somehow exist! Remove it!") 
    else:
        os.system("ln -s %(filename)s temp.bin" %locals())

# data read through fortran (f2py)
    from read_sst import read_sst
    im=1440
    jm=720
    sst=np.zeros([jm,im])

    sst,ierr=read_sst()
    assert ierr==0, "data error in fortran"

    lon=0.125+0.25*np.arange(im)
    lat=-89.875+0.25*np.arange(jm)

    os.system("rm temp.bin")

# range selection
    xrange=np.logical_and(lon>=lonrange[0],lon<=lonrange[1])
    yrange=np.logical_and(lat>=latrange[0],lat<=latrange[1])
    sst=(sst[yrange,:])[:,xrange] 

    return lon[xrange],lat[yrange], sst

#-----------------------------------------------------------------------
if __name__ == '__main__':
    import cdms2,MV2,vcs #cdat

    # read data
    lon,lat,sst=readMGDSST(dir="DATA",year=2011,month=1,day=1,lonrange=(0,360),latrange=(-90,90))

    # hearafter plot with cdat
    # First let's create a variable
    var = MV2.masked_array(sst,mask=(sst<-990.),id="SST")

    # Now create the axis
    lons=cdms2.createAxis(lon)
    lats=cdms2.createAxis(lat)

    # And decorate it
    lons.id = 'longitude'
    lats.id = 'latitude'
    lons.units = 'degrees_east'
    lats.units = 'degrees_north'

    # Just to make sure we will describe this axis as a longitude one
    lons.designateLongitude()
    lats.designateLatitude()

    # Now let's apply these to the variable
    var.setAxisList((lats,lons))

    #var.info() to see info

    #plot 
    v = vcs.init()
    v.plot(var)
    v.png("MGDSST.png")
    del v

上のサンプルメインプログラムでは、2011年1月1日のデータをディレクトリ"DATA"にダウンロードし、経度0から360度、緯度-90度から90度の範囲で読み込んでいる。
その上で、cdatを用いて可視化した。
結果は以下の図。

ダウンロードしたバイナリファイルを読むのにfortranのプログラムをモジュールとしてつかう。清水慎吾さんのデータを読むサンプルプログラムを参考に、以下のようなfortranサブルーチンを作った（配列はpython風に緯度と経度を入れ替えている。欠損値は-999.）。
このプログラムをf2pyを使って以下のようにコンパイルすればpythonのモジュールとして使用が可能になる（上のプログラムのread_sst)。
f2py --f90exec=gfortran -c -m read_sst read_sst.f90
(ここではgfortranを使っている。）

read_sst.ｆ90

subroutine READ_SST(sst,ierr)
  implicit none
  integer   :: year,mon,day
  integer   :: i,j,jj
  
  integer,parameter :: im=1440
  integer,parameter :: jm=720
! 
  real,intent(out):: sst(jm,im)
  integer, intent(out):: ierr
!
  integer   :: isst(im,jm),ist
  real      :: undef = -999.
  
  open(55,file="temp.bin",status='old',form='formatted',iostat=ierr,err=999)
  read(55,*) year,mon,day
  write(6,*) year,mon,day

  do j=1,jm
     read(55,10,iostat=ierr,err=999) (isst(i,j),i=1,im)
  enddo
  close(55)
10 format(1440i3)

  do j=1,jm
     do i=1,im
        jj = jm -j +1
        ist = isst(i,jj)
        if(ist .eq. 888 .or. ist.eq.999) then
           sst(j,i) = undef
        else
           sst(j,i) = real(ist)*0.1 
        endif        
     enddo
  enddo
  return 
999 print *, "Read Error!"

end subroutine READ_SST

Gradsでregrid

2011-02-25T10:37:00.001+09:00

GrADSにおいて、二つのデータセットのグリッドをそろえて、差を見てみることをおこなった。

OpenGrADSのre関数を用いた。
下のソース例、
define sstre=re(sst.1,180,linear,0,2,89,linear,-88,2,ba)
では、変数sst.1を
経度方向に、180グリッド、0度から2度刻みで、線形
緯度方向に、89グリッド、南緯88度から2度刻みで、線形
box averageでグリッドに落としている。

例ではNOAA OISST version 2(1)とNOAA Extended Reconstructed SST version 3b(2)という二つの１月の海面水温データセットから２０１１年１月の値をそれぞれプロットし、3番目の図で、その差(2-1)を引いている。データはNOAAのサイトからOPeNDAPで取得している。

下がgradsのスクリプト。
"grads -p"でportraitモードで使う。
他に、mul.gs, color.gs, draws.gs という Chihiro Kodama氏のgradsスクリプトを用いている。

"reinit"
*NOAA OISST V2
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.oisst.v2/sst.m
nmean.nc"
"set time jan2011"
"mul 1 3 1 3 -yoffset -0.5"
"set grads off"
"color 0 30 2 -kind grainbow -gxout grfill"
"d sst"
"set strsiz 0.15"
"draws 1) NOAA OISST V2"

*NOAA Extended Reconstructed SST V3b
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.ersst/sst.mnme
an.nc"
"set time jan2011"
"mul 1 3 1 2  -yoffset -0.5 "
"set grads off"
"color 0 30 2 -kind grainbow -gxout grfill"
"d sst.2"
"set strsiz 0.15"
"draws 2) NOAA Extended Reconstructed SST V3b"
"cbarn 0.7 1 7.6 7.5"

*diff
"define sstre=re(sst.1,180,linear,0,2,89,linear,-88,2,ba)"
"mul 1 3 1 1  -yoffset -0.5"
"set grads off"
"color -levs -2 -1.5 -1 -0.5 0.5 1.0 1.5 2  -kind blue->white->red -gxout grfill"
"d sst.2-sstre"
"set strsiz 0.15"
"draws 2)-1)"
"cbarn 0.7 1 1. 2."
"printim sst_regrid_compare.png white"

GaLab in Python Interface to GrADS

2011-02-18T19:02:00.001+09:00

Python Interface to GrADSのGaLabを用いてプロットしてみた。
プロットしたデータはEarth System Research Laboratoryで配布されているNCEP/NCAR再解析の2011年１月の海水面気圧をOPeNDAPで取得している。

import matplotlib.pyplot as plt
from grads.galab import GaLab  # for python interface for grads

ga = GaLab(Bin='grads',Window=False,Echo=True,Verb=0,Port=False)

ga.open("http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/slp.mon.mean.nc")

ga.blue_marble("on")
ga("set time jan2011")
ga.basemap('ortho',opts=(130.,40.))

ga.contour("slp")
del ga
plt.title("Sea Level Pressure (Jan2011)")
plt.savefig("galab_slp.png")

Matplotlib Basemap Toolkitよりも手軽にプロットできる半面、緯度・経度線を取り除いたりとかの自由はききにくいようだ。

気象庁全国合成レーダー降雨

2010-12-17T15:45:00.000+09:00

気象庁の全国合成レーダー降雨データがこちらで配布されている。
データは１０分ごとにある。
データはtarでアーカイブされており、その中に含まれるデータはgrib2形式である。
清水慎吾さんがこちらで、データをgrads形式に変換するプログラムを公開している（全国合成気象庁レーダ雨量,GRIB2形式)。

以下のプログラムはファイルをネットから読み込んで、清水さんのツールでgrads形式に変換するpythonスクリプトである。

RadarToGradsはファイル(tarファイル）をネットから呼んで、tarファイルから降雨データのバイナリファイル(grib2)を取り出し、gradsファイルを作る。
grib2からgradsへのファイル変換は清水さんのjmaradar2binで行っている。gradsコントロールファイル("view_radar.ctl"）が自動的にできるようにmain.cで"gradsout=1"としている。

PlotByGradsは作ったgradsデータを使って、gradsでプロットしている。gradsを操作するのに、pythonインターフェイスを使っている。

メインでは、これらを使って、2010年10月20日午前8時００分世界標準時（日本時間午後5時）の降水を東経125-１３５度北緯24-35度でgradsデータに変換して、それをプロットしている。これは奄美大島で豪雨が発生したころにあたる。

import_jmaradar.py

def RadarToGrads(dir=".",gradsbinname="temp.bin",year=2010,mon=1,day=1,hour=0,mi
n=0,slon=118.,elon=150.,slat=20.,elat=48.):
    """
   Import JMA rain radar data (every 10min) from http://database.rish.kyoto-u.ac
.jp/arch/jmadata/synthetic-original.html
   Then convert data to grads file using tool by Dr. Shingo Shimizu (http://shim
izus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%EC%A1%BC%A5%C0%B4%D8%CF%A2%A5%E1%A5%E2)
 

   usage:
       RadarToGrads(dir=".",gradsbinname="temp.bin",year=2010,mon=1,day=1,hour=0
,min=0,slon=118.,elon=150.,slat=20.,elat=48.):
     
   input
      dir: directory to import tar file
      gradsbinname: file name for grads binary
      year, mon, day, hour, min  : time for data (every 10 min) in UTC
      slon,elon  Western,East bounday
      slat,elat  South,North bounday


   output files
      grads binary: gradsbinname
      grads control: view_radar.ctl        
         
    """


    import os
    import tarfile
    import urllib

    # filenames
    syear="%04d" % year
    smon="%02d" % mon
    sday="%02d" % day
    shour="%02d" % hour
    smin="%02d" % min
    date=syear+smon+sday+shour+smin
    tarname="Z__C_RJTD_%(date)s00_RDR_JMAGPV__grib2.tar" %locals()
    urldir="http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/jma-radar/synth
etic/original/%(syear)s/%(smon)s/%(sday)s/" %locals()
    filename="Z__C_RJTD_%(date)s00_RDR_JMAGPV_Ggis1km_Prr10lv_ANAL_grib2.bin" %l
ocals()

    # get tar file
    if not os.path.exists(dir+"/"+tarname):
        if not os.path.exists(dir): os.mkdir(dir) 
        # Download the data
        print 'Downloading data, please wait'
        opener = urllib.urlopen(urldir+tarname )
        open(dir+"/"+tarname, 'w').write(opener.read())
        print 'Downloaded'

    # extract radar precipitation file
    tarfile.open(dir+"/"+tarname, 'r').extract(filename)

    # convert to grads file
    os.system("jmaradar2bin %(filename)s %(gradsbinname)s %(slon)s %(elon)s %(sl
at)s %(elat)s" %locals())
    os.system("rm %(filename)s" %locals())

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

def PlotByGrads(pngname):
    """
      plot converted JMA radar rain data
    """
    from grads.galab import GaLab  # for python interface for grads
    ga = GaLab(Bin='grads',Window=False,Echo=False,Verb=0,Port=False)
    ga.open("view_radar.ctl")
    script="""
    set grads off
    set mpdset hires
    set gxout shaded
    set clevs 0 1 5 10 20 30 50 80
    d rr
    cbarn
    printim %(pngname)s white
"""
    ga(script %locals())


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

if __name__ == '__main__':

# radar file to grads 
    RadarToGrads(dir="DATA",gradsbinname="jmaradar.bin",year=2010,mon=10,day=20,
hour=8,min=00,slon=125.,elon=135.,slat=24.,elat=35.)

# plot by grads
    PlotByGrads("jmaradar.png")

Google Graphで棒グラフ

2010-11-07T20:36:00.000+09:00

GE-GraphというGoogle Graph上に棒グラフをプロットできるフリーソフトがある。http://www.sgrillo.net/googleearth/gegraph.htm
これを使ってみた。

結果から。下の図は正味海洋-大気熱フラックス（海から大気）の１月長期平均をGoogle Earth上で棒グラフにしたものである。正値のみ使っている。

このkmzファイルはこちら。

GE-Graphに入力するデータとして
場所名,緯度,経度,値
で構成されるcsvファイルを用意する必要がある。
上の例のcsvファイルはコレ（その１,その２）
GE-Graphが5000地点以上を受け付けないのでファイルを２つに分けている。場所名は特に意味は無いので全て0を入れている。

上のデータをGE-Graphに読み込ませる。設定は以下の図のようにした。

それぞれのデータで"Run"を行ってkmzファイルを作った。Google Earth上で二つのファイルをまとめて、saveしたのが上のkmzファイルである。

最後にcsvファイルを作るpythonプログラムを添付する。
NCEP/NCARのデータをESRLのサイトからPython Interface to GrADS を使って読み込む。
latent/sensible heat flux, longwave/shortwave radiationを全て足して正味の熱フラックスを計算する。
陸・海マスクを読み、海のデータだけを使う。また正（海から大気）のデータだけ残している。データは4000点づつセーブする。

from grads.ganum import GaNum  # for python interface for grads
import numpy as np

#====== DATA read
ga = GaNum(Bin='grads')
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/lhtfl.mon.ltm.nc") #net latent heat flux
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/shtfl.mon.ltm.nc") #net sensible heart flux
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/nswrs.mon.ltm.nc") #net shortwave radiation
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/nlwrs.mon.ltm.nc") #net longwave radiation
ga("set t 1")  # January
ga("netheat=lhtfl.1+shtfl.2+ nswrs.3+nlwrs.4") # all net heat
netheat=ga.exp("netheat")
ga("close 4")
ga("close 3")
ga("close 2")
ga("close 1")
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/land.sfc.gauss.nc")  #land sea mask
land=ga.exp("land")
del ga

lon=netheat.grid.lon
lat=netheat.grid.lat

#===== DATA manupilation
data=netheat.data
data[land>0]=-999.  # remove values on land

x,y=np.meshgrid(lon,lat)

vl=0.               # data criteria 

xout=x[data>vl].flatten()      # only  values >vl  are used
xout[xout>180]=xout[xout>180]-360. #   -180.=<lon<=180.  
yout=y[data>vl].flatten()      # only  values >vl  are used
zout=data[data>vl].flatten()   # only  values >vl  are used
points=np.zeros(data[data>vl].size) # position number (anythin is ok)

#==== DATA output
nn=zout.size/4000+1
for n in range(nn):
    st=4000*n
    en=4000*(n+1)
    np.savetxt("data%(n)s.txt" %locals(),np.transpose((points[st:en],yout[st:en],xout[st:en],zout[st:en])),fmt='%4i,%7.3f,%7.3f,%7.3f')

Google Earthでアニメーション

2010-10-29T00:12:00.005+09:00

以前、Google Earthで図を描いたが、今回はアニメーションを作ってみる。図が3枚だけのアニメーションだ。

手順１
Gradsで図を3枚つくる。GradsでGoogle Earth向きの図をつくるやり方は以前の記事を参照のこと。下の例はNCEP/NCARの最解析データを使って、2008年から2010年の8月の表面大気温度のアノマリをプロットしている。colorbarはxcbarを使って図の上に一緒に乗せた。

"reinit"
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/air.mon.mean.nc"
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/air.mon.ltm.nc"
year=2008
while (year<=2010)
"set time aug"year
"set x 1 144"
"set y 1 73"
"set parea 0 11 0 8.5"
"set grid off"
"set grads off"
"set mproj scaled"
"set mpdraw off"
"set gxout shaded"
"color -5 5 0.5 -kind blue->gray->red"
"d air.1-air.2(t=8)"
"xcbar  4.5 6.5 4 4.2 -edge triangle -fstep 10 "
"printim  airtemp_"year"aug.png x2048 y1024 -t 0"
"c"
year=year+1
endwhile

こんな感じの図ができる。

手順2
基になるkmlファイルをGradsで作成する。やり方は以前の記事と同じだ。

*************KML
"reinit"
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/air.mon.mean.nc"
"set time aug2010"
"set x 1 144"
"set y 1 73"
"set kml airtemp"
"set gxout  kml"
"d air"

以下のようなkmlファイルが出来る。

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <GroundOverlay>
    <name>Monthly Mean Air Temperature</name>
    <Icon>
      <href>airtemp.tif</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>
</kml>

手順３
kmlファイルを以下の様に書き換える。
ポイントは

・GroundOverlayから/GraoundOverlayの部分を図の数だけ繰り返す。
・それぞれ図の名前を指定する。
・図がどの時間を占めるべきかTimeSpan /TimeSpanの部分で指定する。なお、今の場合、図は8月だけの図であるが、アニメーションで図が続く様に１年中同じ図が表示されるよう、１月１日から12月31日をTimeSpanで設定している。
・繰り返しの全体をFolder と/Folderではさむ。　

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Folder>
  <GroundOverlay>
    <name>Monthly Longterm Mean of Air Temperature</name>
    <TimeSpan>
      <begin>2008-01-01</begin>
      <end>  2008-12-31</end>
    </TimeSpan>
    <Icon>
      <href>airtemp_2008aug.png</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>

  <GroundOverlay>
    <name>Monthly Longterm Mean of Air Temperature</name>
    <TimeSpan>
      <begin>2009-01-01</begin>
      <end>  2009-12-31</end>
    </TimeSpan>
    <Icon>
      <href>airtemp_2009aug.png</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>

  <GroundOverlay>
    <name>Monthly Longterm Mean of Air Temperature</name>
    <TimeSpan>
      <begin>2010-01-01</begin>
      <end>  2010-12-31</end>
    </TimeSpan>
    <Icon>
      <href>airtemp_2010aug.png</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>
</Folder>
</kml>

出来たファイルをkmzにまとめたものをこちらに置いた。kmzファイルはkmlと図をまとめたもので、Google Earthで一度kmlファイルを開いて、「名前をつけて保存」で作成できる。

Google Earthで開くとこんな感じ。

Ubuntu 10.10にアップグレードメモ

2010-10-25T08:54:00.000+09:00

Ubuntuを10.10にアップグレードした。
順調にアップグレードしてくれると思ったらあまかった。
再起動後GUIの画面が始まらず、CUIでログインできるだけ。
リカバリーモードさえうまくいかない。
その夜はそれ以上進まず、諦めて寝ることにした。
次の日ネットで調べたら以下の情報を見つけた。

After upgrade to 10.10 no GUI just command line
http://ubuntuforums.org/showthread.php?t=1575262

これに書いてある中で、 fglrxを消すというのが効果があった。
これで、リカバリーモードのSafe Graphics ModeでGUIが出てくるようになった。
あとは再起動やアップグレードをしているうちに、ちゃんと正常に動くようになった。
良かった、良かった。

富士山の初冠雪日(2) 続・トレンド

2010-10-07T23:48:00.000+09:00

前回の続きとして、私自身あまりよく理解しているとは言えないが、ココを参考にして、ロバスト推定でも線形回帰をとってみた。

参考：
ロバスト推定（統計数理研究所）
http://www.ism.ac.jp/~fujisawa/research/robust.html
ロバスト推定（朱鷺の社）
http://ibisforest.org/index.php?%E3%83%AD%E3%83%90%E3%82%B9%E3%83%88%E6%8E%A8%E5%AE%9A

Rのパッケージrobustbaseを使うために

install.packages("robustbase")

を行っておく。
以下がプログラム。

import numpy as np
import scikits.timeseries as ts
#module
from read_FirstSnowFuji import *

#for R
import rpy2.robjects as robjects

ayear,amonth,aday= read_FirstSnowFuji()
days=np.array([],dtype="i4")
for  n in range(ayear.size):
    td=ts.Date("D",year=ayear[n],month=amonth[n],day=aday[n])
    days=np.append(days,td.day_of_year)


#for R
ryear=robjects.FloatVector(ayear)
rdays=robjects.FloatVector(days)
robjects.globalEnv["ryear"] = ryear
robjects.globalEnv["rdays"] = rdays

r=robjects.r

scripts="""
library("robustbase")   #import robustbase
png('FujiSnow_trend_robust.png')
plot(ryear,rdays,xlab="Year",ylab="Day of year") #scatter plot 
fm<-lmrob(rdays ~ ryear)  # robust linear model
abline(fm)  #plot regression line
abline(lm(rdays ~ ryear),lty=2) # usual linear model
dev.off()
"""
r(scripts)

# summary
summary=r('sm<-summary(fm)')
print summary

以下が図。点線が前回の回帰直線、実線がrobustbaseで得られた直線である。

出力(print summary)は

Call:
lmrob(formula = rdays ~ ryear)

Weighted Residuals:
    Min      1Q Median      3Q     Max
-57.989 -7.872 -0.125   5.986 25.291

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) 76.10515   51.85081   1.468 0.144896
ryear        0.10154    0.02673   3.798 0.000234 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 10.03
Convergence in 12 IRWLS iterations

Robustness weights:
2 observations c(21,115) are outliers with |weight| <= 0.00057 ( < 0.00085);
13 weights are ~= 1. The remaining 102 ones are summarized as
   Min. 1st Qu. Median    Mean 3rd Qu.    Max.
0.00365 0.88600 0.94480 0.88410 0.98680 0.99880
Algorithmic parameters:
tuning.chi         bb tuning.psi refine.tol    rel.tol
1.5476400 0.5000000 4.6850610 0.0000001 0.0000001
nResample     max.it     groups    n.group   best.r.s   k.fast.s      k.max
       500         50          5        400          2          1        200
trace.lev compute.rd
         0          0
seed : int(0)

これによると回帰係数は0.10154（100年で約10日冠雪日が遅くなる）で、しかもp値が非常に小さく統計的に有意である。
ただし、c(21,115)が外れ値、すなわち１９１４年（ayear[21-1]=1914)と２００８年（ayear[115-1]=2008）が外れ値になるのは良いとして、weightが１に近い点は13だけで、あとの外れ値を除く１０２点はweightが１以下に落とされている。

これを見るために、回帰直線の係数、weightなどの数値をpythonに戻し、プロットしてみる（同じことをRでもできるはずだが、pythonのmatplotlibのほうがなれているので）。
以下プログラム。
robust推定で得られた回帰直線を緑線で、
weightが1に近い点を青○で、
weightが0.5から0.9988の間にある点を黒点、
weightが0.00085から0.5の間にある点を赤点、
外れ値を赤xでプロットした。

#weights
weights=r('weights<-sm$weights')
weight=np.array(weights)

#plot
import matplotlib.pyplot as plt
plt.plot(ayear,ayear*a+b,'g-')
range=(weight>0.9988)
plt.plot(ayear[range],days[range],"bo")
range=np.logical_and(weight<=0.9988,weight>=0.5)
plt.plot(ayear[range],days[range],"k.")
range=np.logical_and(weight<0.5,weight>=0.00085)
plt.plot(ayear[range],days[range],"r.")
range=(weight<0.00085)
plt.plot(ayear[range],days[range],"rx")
plt.xlim(ayear[0],ayear[-1]) 
plt.xlabel("year")
plt.ylabel("day of year")
plt.savefig("FujiSnow_trend_robust_2.png")
plt.show()

ちなみに今年２０１０年の重みはweight[-1]=0.86999219949571671
このように多数の点の重みを落とすのは統計的に良くても、物理的な理由は自明ではない。
今回、この方法を用いたのは適当ではなかったかもしれない。

前回と今回の解析で、初冠雪日にトレンドがあるかは統計的にははっきりしなかったが、初冠雪日が遅くなっていっているかもしれないという作業仮説自体は興味深いものがある。これから、将来、どのような傾向があらわれてくるだろうか？

参考
異常気象を追う初冠雪
http://www.bioweather.net/column/essay2/aw15.htm

縮小する富士山の永久凍土
http://www.fujisan-net.jp/data/article/1583.html

富士山の初冠雪日(2) トレンド

2010-10-06T00:25:00.001+09:00

前回作ったモジュール(read_FirstSnowFuji)で初冠雪日を読み、年初からの日にち（1月１日を１とする）に変換したうえで、線形回帰を当てはめてみる。
以前やったように、rpy2を用いて、Rにデータを渡してやる。

import numpy as np
import matplotlib.pyplot as plt
import scikits.timeseries as ts
#module
from read_FirstSnowFuji import *

#for R
import rpy2.robjects as robjects

ayear,amonth,aday= read_FirstSnowFuji()
days=np.array([],dtype="i4")
for  n in range(ayear.size):
    td=ts.Date("D",year=ayear[n],month=amonth[n],day=aday[n])
    days=np.append(days,td.day_of_year)


#for R
ryear=robjects.FloatVector(ayear)
rdays=robjects.FloatVector(days)
robjects.globalEnv["ryear"] = ryear
robjects.globalEnv["rdays"] = rdays
r=robjects.r

scripts="""
png('FujiSnow_trend.png')
plot(ryear,rdays,xlab="Year",ylab="Day of year") # scatter plot
fm<-lm(rdays ~ ryear) # linear model
abline(fm) # plot regression line                 
dev.off()             
"""
r(scripts)
summary=r('sm<-summary(fm)')
print summary

出力(print summary)は

Call:
lm(formula = rdays ~ ryear)

Residuals:
     Min       1Q   Median       3Q      Max
-54.2679 -5.8555   0.9491   8.1001 27.3063

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 138.25153   71.89542   1.923   0.0570 .
ryear         0.06873    0.03683   1.866   0.0645 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.45 on 115 degrees of freedom
Multiple R-squared: 0.0294,    Adjusted R-squared: 0.02096
F-statistic: 3.484 on 1 and 115 DF, p-value: 0.06453

係数が0.06873/yearだから、100年に一週間ほどのトレンド（初冠雪日が遅くなる）となる。ただし、p値が約0.0645であるから危険率5%でトレンドがない可能性を否定できない。データのばらつきが大きいことから、さもありなんといったところである。

富士山の初冠雪日(1) データ取得

2010-10-01T00:04:00.000+09:00

富士山　一気に冬の顔
初冠雪、平年より６日早く
山頂から８合目付近にかけてうっすらと雪をかぶった富士山＝富士吉田市上吉田甲府地方気象台は２５日、富士山が初冠雪したと発表した。平年より６日、昨年より１２日早い。独自に観測している富士吉田市も同日、昨年より１５日早く「富士山初雪化粧」を宣言した。山梨県内は２２日に甲府など３地点で最高気温が３５度以上の猛暑日を記録したばかり。富士山は早くも冬に向けての“衣替え”が始まった。
2010年09月26日山梨日日新聞　 web魚拓

甲府気象台のwebページ表から、富士山初冠雪日を取得するプログラムを作る。
以下がプログラム。
関数 read_FirstSnowFuji()が初冠雪日を取得するモジュールで、年、月、日をnumpyの配列として返す。HTMLから表を読むのにtable parserを用いている。
メインプログラムでは、日付をscikits.timeseriesで年初からの日にち（1月１日を１とする）に変換し、グラフにしている。平均値を黒い実線で、平均値から標準偏差（約14日）の2倍離れた値に黒点線を引いた。ばらつきがかなり大きい印象だ。今年は少しだけ平年より早い。

read_FirstSnowFuji.py

import numpy as np
# table parser from http://simbot.wordpress.com/2006/05/17/html-table-parser-using-python/
from table_parser import *

def read_FirstSnowFuji():
    import urllib
    """
   read first snow cover day of the FUJI mountain from KOFU local meteorological observatory
   usage:
      year,month,day=read_FirstSnowFuji()
   input
      none
   output
      year,month,day:  first snow cover day (numpy array)
   
   """

    f = urllib.urlopen('http://www.jma-net.go.jp/kofu/menu/siryo_hatsukansetsu.dat.html')
    p = TableParser()
    p.feed(f.read())
    f.close()
    data=p.doc # Get to the data

    #get 
    year=[[],[],[]] # 3 column
    month=[[],[],[]] # 3 column
    day=[[],[],[]] # 3 column
    for d in data[0][1:]:
        for n,c in enumerate([0,4,8]):
            if not d[c]=="":
                year[n].append(int(d[c]))
                mi,di=d[c+1].split("/")
                month[n].append(int(mi))
                day[n].append(int(di))
    # to numpy array
    ayear=np.array([],dtype="i4")
    amonth=np.array([],dtype="i4")
    aday=np.array([],dtype="i4")
    for n in range(3):
        ayear=np.append(ayear,year[n])
        amonth=np.append(amonth,month[n])
        aday=np.append(aday,day[n])

    return ayear,amonth,aday

############################################
if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import scikits.timeseries as ts

    ayear,amonth,aday= read_FirstSnowFuji()

    days=np.array([],dtype="i4")
    for  n in range(ayear.size):
        td=ts.Date("D",year=ayear[n],month=amonth[n],day=aday[n])
        days=np.append(days,td.day_of_year)

    #plot
    plt.plot(ayear,days,"b-")
    mean=days.mean()
    stdp=days.mean()+2*days.std()
    stdm=days.mean()-2*days.std()
    plt.plot([ayear[0],ayear[-1]],[mean,mean],"k-")
    plt.plot([ayear[0],ayear[-1]],[stdp,stdp],"k--")
    plt.plot([ayear[0],ayear[-1]],[stdm,stdm],"k--")
    plt.xlim(ayear[0],ayear[-1])
    plt.xlabel("year")
    plt.ylabel("day of year")
    plt.savefig("FujiFirstSnow.png")
    plt.show()

今年8月の500hpa高度場

2010-09-27T23:11:00.000+09:00

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

2010-09-25T20:26:00.000+09:00

人工衛星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で方程式を解く

2010-09-17T22:53:00.000+09:00

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

などとする。

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

2010-09-04T01:57:00.001+09:00

以下のプログラムは、方程式を関数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

2010-09-04T00:17:00.003+09:00

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

2010-09-03T05:44:00.001+09:00

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)

2010-09-01T21:52:00.000+09:00

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

Scipyで常微分方程式を解く (2) 最適化

2010-08-31T21:43:00.000+09:00

前回の応用として、

を解き、

を満たすようなを求める。
真の解は

より

以下がスクリプト。

[1] a が与えられた時に、(x,y)=(0,1)から常微分方程式を積分し、x=x1でのyの値を返す関数y_at_x(a,x1)を定義。
[2] y_at_x(a=1,x1=3)を真の解と比較。
[3] F(a,x1,y1)=y_at_x(a,x1)-y1 という関数を定義。　
[4]　F(a,x1=3,y1=4)=0　を解いてaを求める。aはscipyのbrentq関数を使って0から1の範囲で探す。最後に真の解を使って答え合わせ。

import numpy as np
from scipy import integrate
from scipy import optimize 

#solve dy/dx=a y
# x=0, y=1  
# x=3, y=4
# what is a?

#define derivative
def dy_dx(y,x,a):
    return a*y

#[1]
def y_at_x(a,x1):
    y1, infodict = integrate.odeint(dy_dx, 1., [0,x1],args=(a,),full_output=True)
    #print infodict["message"]
    return y1


# [2] test y_at_x(a,x1)
y_test= y_at_x(a=1.,x1=3.)
print "test y_at_ax(1,3)",y_test[1],np.exp(1.*3.)

# [3] define y_at_x(a,x1)-y1
def F(a,x1,y1):
    return y_at_x(a,x1)[1][0]-y1


# [4] Solve  y_at_x(a,3)-4=0
asolution=optimize.brentq(F,0.,1.,args=(3.,4.))
print "solution of a=",asolution
print "true solution", np.log(4)/3.
print "test solution1",y_at_x(a=asolution,x1=3.)[1][0], np.exp(asolution*3.)

出力は

test y_at_ax(1,3) [ 20.08553873] 20.0855369232
solution of a= 0.462098103279
true solution 0.462098120373
test solution1 4.0 3.99999979487

Scipyで常微分方程式を解く

2010-08-30T14:04:00.001+09:00

Scipyのodeint関数を使い、常微分方程式を解いてみる。

簡単な例として、

を、初期値

で解く。
真の解は

である。

以下がpythonスクリプト。xの0から3まで積分している。数値計算で得られた値を＋で、真の解より得られる値を実線で図示している。

import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt 

#solve dy/dx=y

#define derivative
def dy_dx(y,x):
    return y

#integration
x1=np.linspace(0,3, 10)
x=np.linspace(0,3, 100)
y1, infodict = integrate.odeint(dy_dx, 1., x1, full_output=True)
print infodict

# plot
plt.plot(x1,y1,'+',markersize=12)  #solusion 
plt.plot(x,np.exp(x))            #true solution 
plt.savefig("ode_ex.png")
plt.show()

参照:
LoktaVolterraTutorial
http://www.scipy.org/LoktaVolterraTutorial
A Coupled Spring-Mass System
http://www.scipy.org/Cookbook/CoupledSpringMassSystem

bloggerで数式を書くテスト

2010-08-30T00:07:00.001+09:00

こちらの記事を参照した。
LaTeX for BloggerをSafariと新しい編集インターフェイスのために改造する。
http://satomacoto.blogspot.com/2009/11/latex-for-bloggersafariblogger.html

以下のような数式が書ける。

PyNGLとPyNIOをインストール

2010-08-29T16:48:00.000+09:00

PyNGL(1.31)とPyNIO(1.4.0)をUbuntu10.04にインストールした。

PyNGLは NCAR Command Language (NCL)をもとにしたpythonベースのグラフィックモジュールで、PyNIOはいろいろな形式のファイルを読み込むモジュールである。

以下がインストール手順。

(1)　プレコンパイル済みのファイルを入手
PyNGL-1.3.1.linux-debian-i686-gcc432-py255-numpy130.tar.gz　と　PyNIO-1.4.0.linux-debian-i686-gcc432-py255-numpy130.tar.gz　をEarth System Gridより入手。
Pythonのバージョンは2.6.5を使っているが、numpyは1.3.0を使っているので、プレコンパイル済みのバイナリの中からこれらを選択した。

(2)　以下を実行する。

sudo tar -C /usr/local -xzf PyNIO-1.4.0.linux-debian-i686-gcc432-py255-numpy130.tar.gz
sudo tar -C /usr/local -xzf PyNGL-1.3.1.linux-debian-i686-gcc432-py255-numpy130.tar.gz

(3)　パスを通す。

export PYTHONPATH=$PYTHONPATH:/usr/local/lib/python2.5/site-packages:/usr/local/lib/python2.5/si
te-packages/PyNIO:/usr/local/lib/python2.5/site-packages/PyNGL
export PYNGL_NCARG=/usr/local/lib/python2.5/site-packages/PyNGL/ncarg

(4) pynglexの書き換え
/usr/local/bin/pynglexの一行目を

#!/usr/bin/env python

のように編集する。

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

2010-05-27T22:15:00.003+09:00

東北大学外洋域新世代海面水温（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

2010-05-26T23:38:00.000+09:00

前回までのように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で自己相関

2010-05-15T15:17:00.001+09:00

前回に続いて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で時系列

2010-05-14T18:23:00.000+09:00

以前につくった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-平均法

2010-04-30T21:32:00.000+09:00

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のスペクトル密度を信頼区間つきで

2010-04-27T23:28:00.001+09:00

前回つくった信頼区間つきでスペクトル密度をもとめる関数を、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関数をスペルトル密度誤差も出せるよう拡張

2010-04-26T23:20:00.002+09:00

前回の見延さんのスペクトル密度の推定誤差理論を使って、以前用いた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()

スペクトルの推定誤差

2010-04-20T23:55:00.000+09:00

見延さんのスペクトル誤差に関するmatlab scriptをpythonに移植してみた。

見延　2001:大気海洋物理学特論4 大気海洋統計データ解析　第3章
http://ocw.hokudai.ac.jp/Course/GraduateSchool/Science/MeteorologyAndOceanology/2001/page/materials/MeterologyAndOceanology-2001-Note-03.pdf
Page15-16

numpy/scipy/matplotlibを使えば、matlabのscriptを移植するのはそれほど難しくない。
χ二乗分布の累積分布の逆関数はscipyのstats.chi2.ppfをそのまま用いた。
（注：今の場合関係ないが、この関数は自由度が1以下の場合はうまくいかないバグがあるようだ。すくなくともversion 0.7.0で確認した。）

#coding: utf-8
import numpy as np
from numpy.random import randn
import scipy as sp
import scipy.signal as signal

import matplotlib.pyplot as plt 
#---------------------------------------
def chi2inv_nint(p,df):
    import scipy.stats as stats
    c=stats.chi2.ppf(p,df)
    return c

#------------------------------------------
def psd_ex1(tsrs,nfft,ctaper,avfnc):
    """
function [ds_psd,frq]= psd_ex1(tsrs,nfft,ctaper,avfnc)
Power Spectrum Density を直接法によって求めるpython Script
平滑化は周波数領域で行なう．
複数時系列（sz2$gt;1）についても計算する．
入力
tsrs: 入力時系列 (sz1, sz2)各列についてスペクトルを求める
nfft: fft を計算する時間サンプル数
ctaper: taper 関数('hamming','hanning'など)
avfnc: 平滑化ウインドウ
出力
ds_psd: スペクトル密度．１×nfft/2+1．
frq: 周波数．１×nfft/2+1．
todo
PSD の定義に振幅がきちんと合うようにする．
"""
    if tsrs.ndim==1:
        sz1=tsrs.size #tsrs をsz2*1 のデータに変換．
        sz2=1
        tsrs=tsrs.reshape((sz1,1))
    else:
        sz1,sz2=tsrs.shape
    print sz1,nfft
    assert sz1<=nfft, ' sz1$gt;nfft' #nfft はtlng よりも大きくなくてはならない
#満たさないとエラーメッセージを出して終了
    avfnc=avfnc/float(np.sum(avfnc)) #後で説明する平滑化ウインドウの規格化    
    w = getattr(signal, ctaper)(sz1) #実際のtaper 計算
    for n in range(sz2):
        tsrs[:,n]=tsrs[:,n]*w #時系列にテーパーをかける
    tsrs=np.append(tsrs,np.zeros((nfft-sz1+1,sz2)),axis=0) #ゼロパッティング
    famp=sp.fft(tsrs,axis=0) #フーリエ変換
    famp=famp[0:(nfft/2+1),:]/sz1 #実数時系列のスペクトルなら，この範囲のみが必要．
    d_psd=abs(famp)**2; #直接スペクトル値 (d: direct)
    ds_psd=np.zeros_like(d_psd)
    for n in range(sz2):
        ds_psd[:,n]=np.convolve(d_psd[:,n],avfnc,mode='same'); #畳込みを使って平滑化 (ds: direct-smoothed) #畳み込みで増えた部分を除く
    frq=np.arange(nfft/2+1)/float(nfft); #周波数を求めておく
    return ds_psd,frq

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

if __name__ == '__main__':
    plt.figure()
    for npnl in range(4):
        tlng=100; nsmpl=1000
        if npnl==0 :
            nfft=100
            taper=signal.boxcar(tlng)
            ctaper='boxcar'
            avfnc=np.asarray([1, 1, 1, 1, 1]);
            cttl='no padding, boxcar, 5-running'
        elif npnl==1:
            nfft=100
            taper=signal.hamming(tlng)
            ctaper='hamming'
            avfnc=np.asarray([1, 1, 1, 1, 1])
            cttl='no padding, hamming, 5-running'
        elif npnl==2:
            nfft=200
            taper=signal.boxcar(tlng)
            ctaper='boxcar'
            avfnc=np.asarray([1, 1, 1, 1, 1]);
            cttl='double padding, boxcar, 5-running'
        elif npnl==3:
            nfft=200
            taper=signal.hamming(tlng)
            ctaper='hamming'
            avfnc=np.convolve([1, 2, 1],[1, 2, 1])
            cttl='double padding, hamming, double 1-2-1'
    
        tsrs=randn(tlng,nsmpl);
        [ds_psd,frq]=psd_ex1(tsrs,nfft,ctaper,avfnc)
# Monte-Carlo 法で得る90%信頼限界
        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])
        avfnc=avfnc/float(sum(avfnc)) # frequency domain の平滑化関数の規格化
        c0=srt_psd[:,nsmpl*0.5]

        edof=(1+(1/sum(avfnc**2)-1)*tlng/nfft*np.mean(taper))# 等価自由度の計算
        ce=np.zeros([frq.size,2]) 
        ce[:,0]=np.log10(c0*chi2inv_nint(0.050,2*edof)/chi2inv_nint(0.50,2*edof))
        ce[:,1]=np.log10(c0*chi2inv_nint(0.950,2*edof)/chi2inv_nint(0.50,2*edof))
#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_ex1.png")
    plt.show()

結果の図

ENSO の周期　psd関数のオプション5 scale_by_freq

2010-04-10T00:18:00.001+09:00

前回の続き。
これは単なるテクニカルなオプションだが、メモのために記しておく。

*scale_by_freq*: boolean
Specifies whether the resulting density values should be scaled by the scaling frequency, which gives density in units of Hz^-1. This allows for integration over the returned frequency values. The default is True for MatLab compatibility.

今まで使ってきた下のような例では、もとは月単位のデータを年単位のデータにしたときに(Fs=12)、自動的に12で割っている（scale_by_freq=Trueがデフォルト)。
すなわち、周波数が(1/2,1/3,1/4, ......)（月単位）の時のスペクトル密度（Pxx2,scale_by_freq=False)の場合と、周波数が（6,4,3, ....)(年単位）の時のスペクトル密度Pxx(scale_by_freq=True)の積分が一致するように

Pxx=Pxx2/12

の関係がある。
下のグラフのようにPxx(黒線）とPxx2/12(赤線）とは一致する。

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)  
Pxx2,f2=mlab.psd(x, NFFT=len(x),Fs=12,scale_by_freq=False) 

plt.figure()
plt.loglog(f,Pxx,'k')
plt.loglog(f2,Pxx2/12.,'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_scale.png")
plt.show()

ENSO の周期　psd関数のオプション4 平滑化

2010-04-08T23:16:00.000+09:00

前回の続き。
これはpsd関数のオプションというわけではないが、Welch法ではなく、得られたスペクトル密度を畳込みによって平滑化してみる。
平滑化にはscipyのcookbookの例にあるsmooth関数を用いる。
平滑化にはスペクトルの周波数の解像度を下げるかわりに統計的自由度を上げる働きがある（見延、2001; Hartmann 2008)。
平滑化には1-2-1フィルタを用いる。numpyでは5点bartlett関数を使うと良い（両端が0）。

import numpy as np
print np.bartlett(5)

[ 0. 0.5 1. 0.5 0. ]

以下がプログラム例。
黒線がもとのスペクトル密度。
赤線が平滑化したスペクトル密度。
（れいによって2年目と7年目の周期に当たるところに縦の点線）

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"]

#original
Pxx,f=mlab.psd(x, NFFT=len(x),Fs=12)  

# Smooth by 1-2-1 filter 
from cookb_signalsmooth import smooth # from Scipy cookbook
Pxx2=smooth(Pxx[:,0],5,'bartlett') 


#plot
plt.figure()
plt.loglog(f,Pxx,'k-o')
plt.loglog(f,Pxx2,'r-+')
plt.loglog([1./2.,1./2.],[1.E-6,1.E1],'k--') #vertical line for 2 years
plt.loglog([1./7.,1./7.],[1.E-6,1.E1],'k--') #vertical line for 7 years
plt.savefig("NINO_psd_smooth.png")
plt.show()

参考文献
見延　2001:大気海洋物理学特論4 大気海洋統計データ解析　第3章
Hartmann 2008:ATM552 講義ノート第6a章

ENSOの周期　psd関数のオプション3 Welch法

2010-04-08T00:07:00.000+09:00

前回の続き。
ここで用いているpsd関数は、データ長をブロックにわけ、その平均をおこなうことでスペクトル密度をスムーズにするオプションがある(Welch法）。

The power spectral density by Welch’s average periodogram method. The vector x is divided into NFFT length blocks.
noverlap: integer The number of points of overlap between blocks. The default value is 0 (no overlap).

以下がその例。
黒線がオリジナル。
赤線が、データ長を二分割したケース（オーバラップ0）。
青線が、データ超を二分割した上で、50%のオーバラップを許した例。
スペルトル密度の線がなめらかになるのがわかる。
(補助線として2年と7年に対応する周期に縦点線をひいている。）

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"]

#original
Pxx,f=mlab.psd(x, NFFT=len(x),Fs=12)  

# Date are devided into 2 blocks
#Zero pad to the length of the orignal data sequence.
Pxx2,f2=mlab.psd(x, NFFT=len(x)//2,pad_to=len(x),Fs=12)

#Date are devided into 2 blocks
#Allow 50 % overlap 
#Zero pad to the length of the orignal data sequence.
Pxx3,f3=mlab.psd(x, NFFT=len(x)//2,noverlap=int(0.5*len(x)/2.),pad_to=len(x),Fs=
12)

#plot
plt.figure()
plt.loglog(f,Pxx,'k-o')
plt.loglog(f2,Pxx2,'r-+')
plt.loglog(f3,Pxx3,'b-x')
plt.loglog([1./2.,1./2.],[1.E-6,1.E1],'k--') #vertical line for 2 years
plt.loglog([1./7.,1./7.],[1.E-6,1.E1],'k--') #vertical line for 7 years
plt.savefig("NINO_psd_welch.png")
plt.show()

参考
pylab.psd関数のデモ
pylab_examples example code: psd_demo2.py
MATLABのWelch法の関数pwelch
http://www.mathworks.com/access/helpdesk_ja_JP/help/toolbox/signal/pwelch.html

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

2010-03-30T22:15:00.001+09:00

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

pad_to: integerThe number of points to which the data segment is padded whenperforming the FFT. This can be different from NFFT, whichspecifies the number of data points used. While not increasingthe actual resolution of the psd (the minimum distance betweenresolvable peaks), this can give more points in the plot,allowing for more detail. This corresponds to the n parameterin the call to fft(). The default is None, which sets pad_toequal to NFFT

下がその例である。
差を良く見るために２年から７年の周期（周波数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()

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

2010-03-29T23:43:00.002+09:00

前回使った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章

Matplotlib for Python Developers

2010-03-28T13:30:00.000+09:00

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

ENSOの周期

2010-03-12T22:43:00.000+09:00

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

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

2010-03-06T16:20:00.002+09:00

寒暖の差大きく、大雪も　12～２月、気象庁まとめ　（日経新聞)
気象庁は１日、この冬（昨年12月～今年２月）の天候まとめを発表した。３カ月間の平均気温は全国的に高めとなったものの、暖冬をもたらす「エルニーニョ現象」と北極域の気圧が高くなる「北極振動」の影響で、例年以上に寒暖の差が大きくなった。　気象庁によると、今冬はエルニーニョ現象の影響で西高東低の冬型の気圧配置が長続きせず、平均気温は西日本で平年より1.0度、東日本で0.9度、北日本と沖縄・奄美は0.6度高くなった。　季節外れの暖かさとなったのは１月下旬と２月下旬で、２月25日には全国約160カ所で、２月の最高気温の記録を更新した。(01日　19:53)

エルニーニョを話題とする。エルニーニョを監視するために、東太平洋のエルニーニョ監視領域の海面温度がよく使われる。エルニーニョ監視領域
エルニーニョ監視領域の温度とその平年値からの偏差(anomaly)のデータはNOAAのclimate predictionセンターから手に入る。ここの情報によれば、平年値は1971年から2000年までの平均で定義される。

以下は、これらのエルニーニョ監視領域の絶対値と偏差をネットから取得するモジュールと、４つの監視領域の偏差をデータ全期間もしくは過去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-03-04T21:32:00.004+09:00

以前、串本と浦神の潮位差の月平均データを読んだが、今回は日データを読む。
データは気象庁のサイトで手に入れることができる(速報値）。
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()

串本・浦神の潮位差と黒潮流路緯度　Scipyでregression

2010-02-17T00:11:00.002+09:00

regressionはscipyでも行える。stats.linregress関数（実際にはmasked array用のstats.mstat.linregress)を使う。
参照　Cookbook/Linear Regression

#read module
from scipy import stats

from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

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

#calculate regression
(a_s,b_s,r,prob,stderr,c)=stats.mstats.linregress(kushimoto_uragami,kuroshio_lat)
print "slope=",a_s
print "intersect=",b_s
print "correlation=",r
print "p value=",prob
print "std. of residuals",stderr
print "(Var. Latitude)/(Var. Sea level)=",c

最後のcはマニュアルには無いが、ソースコードを見てみたところ(ipythonではstats.mstats.linregress??で見ることができる。scipy version 0.7.0）、それぞれの変数の分散の比になっていた。結果は

slope= 0.0666657
intersect= 31.9559306074
correlation= 0.618468
p value= 4.38452479663e-62
std. of residuals 0.62129592706
(Var. Latitude)/(Var. Sea level)= 0.011619

regressionの線も引いてみよう。

プログラムはmatplotlibを使って以下のようになる。

#plot
import matplotlib.pyplot as plt
import numpy as np

x=np.arange(-7.,35.,1.) # x data from -7 to 35 step 1
y=a_s*x+b_s         # regression line

plt.figure()
plt.plot(kushimoto_uragami,kuroshio_lat,'.') # plot data
plt.plot(x,y,'r-') # plot regression line
plt.xlabel("Kushimoto-Uragami")
plt.ylabel("Kuroshio Latitude")
plt.xlim((-7.,35.))
plt.ylim((30.,33.5))
plt.savefig("kurosho_reg_scipy.png")
plt.show()

今回はつかわなかったが、以下の情報もある。
Cookbook/Fit Statistics

Scikits.statsmodels

multiple regression

串本・浦神の潮位差と黒潮流路緯度　Rで線形モデル　（続き2）

2010-02-14T23:36:00.002+09:00

前回のさらに続き。推定値・予測値の信頼区間も描くことができる。

前々回もとめた係数をα、切片をβとすると、
(1) -7から35を1刻みで分割した串本・浦神の各点x0
(2) αx0+βの推定値、およびその99%信頼区間
(3) αx0+βの推定値、および値αx0+β+ε0（ε0はx0で混入する誤差）の予測値の99%信頼区間
(4) αx0+βの推定値および予測値の信頼区間を黒、推定値の信頼区間を青でプロット　（計6本の線。直線ではない）。
(5) グラフを重ねがきする。
(6) データ点をプロットする。

rscript="""
new<-data.frame(rlevel=seq(-7,35,1))  # (1)
png(\'kuroshio_reg_R3.png\')
dc<-predict(fm,new,interval="confidence",level=0.99) #(2)   
dp<-predict(fm,new,interval="prediction",level=0.99) #(3)
matplot(new$rlevel,cbind(dp,dc),lty=1,col=c("black","black","black","black","blue","blue"),type="l",xlab="",ylab="",xlim=c(-7, 35),ylim=c(30,33.5))　 #(4)
par(new=T) #(5)
plot(rlevel,rlat,xlab=\"Kushimoto-Uragami\",ylab=\"Kuroshio Latitude\",xlim=c(-7, 35),ylim=c(30,33.5))  #(6)
dev.off()     
"""
r(rscript)

図

参考:　やさしい医学統計手法　6.二標本の関連性を見る。

http://www3.ocn.ne.jp/~stat/medical/med_023.htm

串本・浦神の潮位差と黒潮流路緯度　Rで線形モデル　（続き）

2010-02-13T10:36:00.002+09:00

前回の続き。
以下のようにして残差(residuals)に関するプロットを作ることができる。
(1) 図を四分割し、上下左右の間隔を指定。現在の作図パラメータはoparに退避。
(2）4種類の診断図を一度に描く。
(3) 作図パラメータをもとにもどす。

r('png(\'kuroshio_reg_R2.png\')')
r('opar<-par(mfrow=c(2,2),oma=c(0,0,1.1,0))')　　#(1)
r('plot(fm)')                                       #(2)
r('par(opar)')                                      #(3)  
r('dev.off()')

図

左上：　予測値対残差。残差が1.96σを超える（正規分布で95%内に入っていない）データにはデータ番号が示されている。

右上：　標準化された残差に対する累積比率に対する正規Q-Qプロット。残差が正規分布しているか。

左下：　残差の標準化の絶対値の平方根を予測値に対してプロットしている。

右下：　leverage(梃子比: 回帰直線への影響度）と標準化された残差。

例えば、右上は黒潮の緯度が高い（日本に近い）場合には、回帰直線は過大評価になっている。これは陸地があるので、それ以上は高い位置にいけないからであろう。緯度が低い（沿岸から離れる）場合に値がばらついているのは、黒潮が十分に岸から離れている場合には串本・浦神の潮位差は0に近く、黒潮の離岸距離との関係がうすまっていると考えられる。

串本・浦神の潮位差と黒潮流路緯度　Rで線形モデル

2010-02-11T17:52:00.016+09:00

前回の続き。今回はrpy2を通してRを使い、串本・浦神の潮位差と黒潮流路緯度を線形モデルとして評価してみる。
(1) データをRに渡し、scatter plotを描く。
(2) 潮位差(rlevel)と黒潮緯度(rlat)の線形モデルを考える。
(3) 線形モデルをもとに、regression line の直線を引く。
(4) 線形モデルの結果を打ち出す。
missing valueもmasked arrayから自動的に選別して除いてくれているようだ。

#read module
#for data
from read_kuroshio_lat import read_kuroshio_lat
from read_kushimoto_uragami import read_kushimoto_uragami

#for R
import rpy2.robjects as robjects

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


#for R
rlat=robjects.FloatVector(kuroshio_lat)
rlevel=robjects.FloatVector(kushimoto_uragami)
robjects.globalEnv["rlat"] = rlat
robjects.globalEnv["rlevel"] = rlevel
r=robjects.r
r('png(\'kuroshio_reg_R.png\')')    
r('plot(rlevel,rlat,xlab=\"Kushimoto-Uragami\",ylab=\"Kuroshio Latitude\")') #scatter plot 
r('fm<-lm(rlat ~ rlevel)')  # linear model
r('abline(fm)')             # plot regression line
r('dev.off()')            
summary=r('sm<-summary(fm)')
print summary

図

summaryの出力

Call:
lm(formula = rlat ~ rlevel)

Residuals:
     Min       1Q   Median       3Q      Max
-3.17592 -0.32259 0.08407 0.40907 1.30406

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 31.955937   0.038044 839.98   <2e-16 ***
rlevel       0.066666   0.003672   18.15   <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.6225 on 532 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.3825,   Adjusted R-squared: 0.3813
F-statistic: 329.5 on 1 and 532 DF, p-value: < 2.2e-16

線形モデルの係数は0.066666、切片は31.955937。
ともにp値は小さく(2e-16)、0であるという帰無仮説は棄却できる。
線形モデルのR^2値は約38%で、前回求めた相関係数0.62の2乗になっている。

係数などは以下のようにとりだせる。

print summary.names

[1] "call"          "terms"         "residuals"     "coefficients"
[5] "aliased"       "sigma"         "df"            "r.squared"
[9] "adj.r.squared" "fstatistic"    "cov.unscaled" "na.action"

cf=r('cf<-sm$coefficients')
print cf.names

[[1]]
[1] "(Intercept)" "rlevel"

[[2]]
[1] "Estimate" "Std. Error" "t value" "Pr(>|t|)"

for i in range(8):
print cf[i]

31.9559366577
0.066665747972
0.0380437732903
0.00367237036394
839.978106637
18.1533291486
0.0
1.13483941203e-57

参考：単回帰http://www.geocities.jp/kashiwabara_fact/statistics-tankaiki.htm

関連ポスト
北極振動　Rで統計
北極振動　Rで統計グラフ編

串本と浦神の潮位差と黒潮流路位置の相関

2010-02-10T23:49:00.001+09:00

串本と浦神の潮位差が黒潮の流路のよい指標となることは経験的に知られている（参考:気象庁）。
前々回作った黒潮の流路の緯度を読むプログラムと、前回作った串本と浦神の潮位差を読むプログラムをモジュールとして用いて、潮位差を横軸に、流路緯度を縦軸にプロットしてみた。

#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

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

#plot
plt.figure()
plt.plot(kushimoto_uragami,kuroshio_lat,'.')
plt.xlabel("Kushimoto-Uragami")
plt.ylabel("Kuroshio Latitude")
plt.savefig("kurosho_corr.png")
plt.show()

潮位差が大きいと、黒潮の緯度は高くなる（岸に近い流路をとる）傾向があることがわかる。
相関をとってみよう。

c1=np.corrcoef([kushimoto_uragami,kuroshio_lat])
print c1
[[ 1. 0.11969213]
[ 0.11969213 1. ]]

相関係数は 0.11969213　!? 小さ！！！

numpyのmasked arrayでマスクした潮位差の欠損値のところでうまく計算できていないらしい。ちゃんと欠損値を抜いて計算してみよう。

mask=(kushimoto_uragami.mask==False)
c2=np.corrcoef([kushimoto_uragami[mask],kuroshio_lat[mask]])
print c2

[[ 1. 0.61846868]
[ 0.61846868 1. ]]

今度はそれなりの相関係数 0.61846868 が出た。

実は初めからmasked array用の相関係数の関数ma.corrcoefを使えば、そのままでも計算できる。

c3=np.ma.corrcoef([kushimoto_uragami,kuroshio_lat])
print c3

[[1.0 0.61846868176]
[0.61846868176 1.0]]

串本と浦神の潮位差

2010-02-09T22:32:00.002+09:00

串本と浦神の潮位差は黒潮の流路と密接な関係があることが知られている（参考:気象庁）。串本と浦神の潮位差の月平均値のデータは気象庁から手に入れることができる。http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/ksur.txt

以下はそれを読むプログラム。基本的には前回と同じだが、欠損値があるので一手間が必要。

read_kushimoto_uragami.py

import numpy as np
import scikits.timeseries as ts

def read_kushimoto_uragami():
    day=[]
    value=[]
    ds = np.DataSource(None)
    f=ds.open('http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/ksur.txt')
    for line in f:
         fields = line.strip().split()
         if len(fields)>=1: day.append(fields[0])
         if len(fields)==1:
             value.append(-999.)
         elif len(fields)==2:
             value.append(fields[1])
    f.close()
    valuea=np.asarray(value,dtype='float32')
    kushimoto_uragami=ts.time_series(valuea,start_date=day[0],freq='M',mask=(valuea<-990.))
    return kushimoto_uragami

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

    kushimoto_uragami=read_kushimoto_uragami()

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

黒潮の流路変化

2010-02-08T23:00:00.002+09:00

*
*

初ガツオ：水揚げ　勝浦漁港、活気づく　／千葉
勝浦市の勝浦漁港に２５日、初ガツオが水揚げされた。江戸中期の俳人・山口素堂が「目には青葉山ほととぎす初鰹」と詠む通り、黒潮に乗って日本列島に近づく初ガツオの本格シーズンは５月ごろ。冬の終わりを予感させる水揚げに、漁港は活気づいた。
（毎日新聞 2010/1/25）

黒潮の流れは、日本南岸を岸に近づいたり、離れたりして流れる。
参考：黒潮の数か月から十年規模の変動（流路）　気象庁

黒潮の流軸の緯度のデータは上の気象庁のサイトから入手できる。

以下はそのデータを読取りプロットするpythonプログラム。

import numpy as np
import scikits.timeseries as ts

def read_kuroshio_lat():
    ds = np.DataSource(None)
    f=ds.open('http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/kuro_slat.txt')
    data=np.genfromtxt(f,dtype=[('date','S7'),('lat','f4')])
    f.close()
#    ds.__del__
    kuroshio_lat=ts.time_series(data['lat'],start_date=data['date'][0],freq='M')
    return kuroshio_lat

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

    kuroshio_lat=read_kuroshio_lat()

    fig = tpl.tsfigure()
    fsp = fig.add_tsplot(111)
    fsp.tsplot(kuroshio_lat, '-')
    plt.xlabel("YEAR")
    plt.ylabel("Kurosho latitude")
    plt.savefig("kurosho_lat.png")
    plt.show()

北極振動時の日本付近の冬の気圧パターン

2010-02-08T12:05:00.001+09:00

*
*
12月の1000hpaのgeopotentail heightに前にもとめた北極振動のパターンを重ねてみる。北極振動指数が正だと西高東低が弱まり、負だと強まる。

北極振動指数0の時

北極振動+2の時

北極振動-2の時

以下、プログラム

from mpl_toolkits.basemap import Basemap # for Basemap toolkit
import matplotlib.pyplot as plt
import numpy as np

lat=np.load('lat.npy') # latitude
lon=np.load('lon.npy') # longitude
lat[-1]=lat[-1]-1.E-4 # lat=90 does not work

#average field
mhgt=np.load('mhgt.npy') # seasonal climatology
dec_average=mhgt[11,:,:] # december

#AO pattern
factor = np.sqrt(np.cos(np.pi*lat/180.))
eof1=np.load('eof1.npy')
AO_pattern=eof1/factor[:,np.newaxis]

#------------- plot geopotential height (AO index 0, 2,-2)
for index in [0.,2.,-2.]:
    data=dec_average+AO_pattern*index
    fig=plt.figure()
    ax = fig.add_axes([0.1,0.1,0.8,0.8])
    map = Basemap(width=15.e6,height=12.e6,\
            projection='gnom',lat_0=60.,lon_0=135.)
    map.drawcoastlines()

    lons,lats=np.meshgrid(lon[:],lat[:])
    x,y=map(lons,lats)
    cy=map.contourf(x,y,data,np.arange(-100,300,20),extend='both')
    cz=map.contour(x,y,data,np.arange(-100,300,20),colors='k')
    plt.clabel(cz,fmt='%3.0f',)
    plt.savefig('AO_patter_'+str(index)+'.png')
    plt.show()

北極振動指数の時系列　2/2

2010-02-07T16:23:00.000+09:00

*
*
前回はEOFと同時に得られる時系列をプロットした。今回はEOFの期間外の北極振動指数を求めてプロットする。例えば2008年から2009年までの値を求めてプロットすると次のとおり。赤が自分で求めた値で、黒がNOAAの発表している値。ほとんどあっている。

以下がプログラム。

(1) python interface for grads を使い、2008年から2009年の1000hpaのgeopotential heightを読み込む。データソースはNCEP/NCAR reanalysis。

(2) EOFを求めるときにセーブしておいた、12ヶ月分の1979-2000年気候値(mhgt)を読み込んで、差し引く。mhgtはnumpyのtile関数を用いて時間方向に2回(2年分)の繰り返しにする。

(3)　緯度のファクターを掛ける。緯度90度の取扱いはEOFを求める時とおなじようにしておく。

(4)EOFを求めるときにセーブしておいたEOF1のパターンにプロジェクションをおこなう。また、北極振動を規格化するための定数pfactor1で、EOFを規格化するのと、時系列を規格化するので2回割る。

from grads.ganum import GaNum # for python interface for grads

from mpl_toolkits.basemap import Basemap # for Basemap toolkit
import matplotlib.pyplot as plt
import numpy as np

import scikits.timeseries as ts #scikits timeseries module
import scikits.timeseries.lib.plotlib as tpl

#------------ data import using python interface for grads　(1)

ga = GaNum(Bin='grads')
fh=ga.open("http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/hgt.mon.mean.nc") # data open
#fh=ga.open("/home/tmiyama/DATA/Reanalysis/NCEP1/hgt.mon.mean.nc") # data open
ga("set lat 20 90")
ga("set time JAN2008 DEC2009") # use data from 2008 to 2009
ga("set lev 1000") # 1000 hpa
hgt = ga.exp("hgt")
del ga

#------------- Get AO index ----------------------------

#remove seasonal climatogical (2)
mhgt=np.load('mhgt.npy')
hgt=hgt-np.tile(mhgt,(2,1,1))

# multiply latitude factor　(3)　
lat=hgt.grid.lat
lat[-1]=lat[-1]+1.E4
factor = np.sqrt(np.cos(np.pi*(lat/180.)))
datain = hgt[:,:,:] * factor[np.newaxis,:,np.newaxis]

# projection to EOF1 (4)
nt,ny,nx=datain.shape
eof1=np.load('eof1.npy')
pc1=np.dot(np.reshape(datain.data,(nt,ny*nx)),np.reshape(eof1,(ny*nx,1)))
pfactor1=np.load('pfactor1.npy')
pc1=pc1/pfactor1/pfactor1

dbeginmonth=ts.Date('M', '2008-1')
AO_2008_2009_series=ts.time_series(pc1,start_date=dbeginmonth,freq='M') # make timeseries object

# plot
from readAO import readAO
AOindex_series=readAO() # import NOAA data

fig = tpl.tsfigure()
fsp = fig.add_tsplot(111)
l1=fsp.tsplot(AOindex_series, 'k-')
l2=fsp.tsplot(AO_2008_2009_series, 'r-')
istart=ts.Date('M','2008-1').value
iend=ts.Date('M','2009-12').value
fsp.set_xlim(istart,iend)
plt.xlabel("YEAR")
plt.ylabel("AO index")
plt.legend((l1,l2),("NOAA","What I got"))
plt.savefig("check_AOindex-2008-2009.png")
plt.show()

北極振動指数の時系列　1/2

2010-02-04T23:31:00.000+09:00

*
*
前回の続き。前回もとめたfirst EOFの時系列が、NOAA CPCの発表している北極振動の値と一致しているか確かめる。黒がNOAAのCPCの値で、赤が自分でもとめたもの。
ちょっとだけ値がずれているが、ほとんどあっているのでいいんじゃないかな。前回よくわからなかった北極のとりあつかいとかがずれる原因かもしれない。

以下がプログラム。以前作ったreadAOモジュールでNOAA CPCの北極振動を読む。時系列オブジェクトを作るのと、プロットにはscikits.timeseriesを用いた。

# plot pc1
import scikits.timeseries as ts
import scikits.timeseries.lib.plotlib as tpl

dbeginmonth=ts.Date('M', '1979-1')
pc1_series=ts.time_series(pc1,start_date=dbeginmonth,freq='M')

from readAO import readAO
AOindex_series=readAO()

fig = tpl.tsfigure()
fsp = fig.add_tsplot(111)
l1=fsp.tsplot(AOindex_series, 'k-')
l2=fsp.tsplot(pc1_series, 'r-')
istart=ts.Date('M','1979-1').value
iend=ts.Date('M','2000-12').value
fsp.set_xlim(istart,iend)
plt.xlabel("YEAR")
plt.ylabel("AO index")
plt.legend((l1,l2),("NOAA","What I got"))
plt.savefig("check_AOindex.png")
plt.show()

北極振動　パターン

2010-02-04T00:03:00.004+09:00

*
*
NOAA Climate Prediction Centerの記述　http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/history/method.shtml　を参考にEOF解析を用いて北極振動のパターンをもとめる。

用いる”python interface for grads”にはEOF関数があるのだが、計算がうまくいかなかったので(緯度が90度の点があるとうまくいかないようだ）、自分でつぎのようなモジュールを作っておく。
”python interface for grads”から

from grads.ganum import GaNum
ga = GaNum(Bin='grads')
ga.eof??

で関数の中身がわかるのでそれを参考に作る。
eof.py

def eof (data):
    """
        calculates Empirical
        Orthogonal Functions (EOFS) using Singular Value Decomposition (SVD).
            V, d, c = eof(data)
        where
            V --- eigenvectors
            d --- NumPy array with eigenvalues
            c --- NumPy array with principal components
    """
    import numpy as np
#       At least 2 time steps
#       ---------------------
    if data.shape[0] < 2:
        raise GrADSError, \
                  'need at least 2 time steps for EOFS but got nt=%d'% data.shap
e[0]
    nt, ny, nx = (data.shape[0], data.shape[1], data.shape[2])
#       Export N-dimensional array
#       --------------------------
    u = data
#       Reshape as 3D
#       -------------
    u = u.reshape(nt,ny,nx)
#       ------------------------------------------
    pc, d, u = np.linalg.svd(u.reshape(nt,ny*nx),
                                full_matrices=0)
#       Eigenvalues/coefficients
#       ------------------------
    d = d * d
    pc = pc.transpose()
#       Reshape eigenvectors as original array
#       --------------------------------------
    u = u.reshape(nt,ny,nx)
    return (u, d, pc)

モジュールが準備できたら、EOF解析を実行する。
使用するモジュールをインポートする。最後に、上のeofモジュールをインポートしている。

from grads.ganum import GaNum # for python interface for grads

from mpl_toolkits.basemap import Basemap # for Basemap toolkit
import matplotlib.pyplot as plt
import numpy as np

from eof import eof # import EOF

python interface for gradsを用いて、1971年か2000までのデータを読み込む。データはOPeNDAPを用いてNCEP/NCARのデータをNOAA ESRLのサイトから読み込んでいる。あるいはローカルにnetcdfを落としておいてそこから読んでも良い。1000hpaのGeopotential Heightを用いる。Bin GuanさんのGradsスクリプトを用いて季節変動(mhgt)を引いて、アノマリ（ahgt)だけを取り出す。

#------------ data import using python interface for grads

ga = GaNum(Bin='grads')
fh=ga.open("http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/hgt.mon.mean.nc") # data open
#fh=ga.open("hgt.mon.mean.nc") # data open
ga("set lat 20 90")
ga("set time JAN1979 DEC2000") # use data from 1979 to 2000
ga("set lev 1000") # 1000 hpa
ga("deseason hgt ahgt mhgt JAN1979 DEC2000") #use "deseason" at http://www.atmos.umd.edu/~bguan/grads/GrADS_Scripts.html agt is anomaly
ahgt = ga.exp("ahgt")
mhgt = ga.exp("mhgt")
del ga

EOF解析をおこなう。cos（緯度)^0.5のファクターをかけるのであるが、緯度90度の部分があると、svdがもとまらなかった。どうするのが正しいかよく分からないが、ちょっとだけ緯度をずらしてみた。北極振動指数は、1979年から2000年までの時系列の分散が1になるように設定されているので、そのように規格化する。また通常の定義にあうように符号も逆転させておく。

#------------- EOF ----------------------------
lat=ahgt.grid.lat # latitude
lon=ahgt.grid.lon
lat[-1]=lat[-1]-1.E-4 # lat=90 does not work
factor = np.sqrt(np.cos(np.pi*lat/180.))
datain = ahgt[:,:,:] * factor[np.newaxis,:,np.newaxis]
v,d,c=eof(datain)
pc1=c[0,:]*np.sqrt(d[0])
pfactor1=pc1.std(axis=0)
pc1=-pc1/pfactor1
eof1=-v[0,:]*pfactor1

あとはプロット。cos（緯度)^0.5ももとに戻しておく。後で使うかもしれないのでデータも保存しておく。

#------------- plot loading pattern
fig=plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8])
map = Basemap(projection='nplaea',boundinglat=20,lon_0=0.,resolution='i')
map.drawcoastlines()

lons,lats=np.meshgrid(lon[:],lat[:]+1.E-4)
x,y=map(lons,lats)
cz=map.contourf(x,y,eof1/factor[:,np.newaxis],[-40,-35,-30,-25,-20,-15,-10,-5,5,10,15,20],extend='both',cmap=plt.cm.RdBu_r)

# get axes position, add colorbar axes to right of this.
pos = ax.get_position()
l, b, w, h = pos.bounds
cax = plt.axes([l+w+0.05, b, 0.03, h]) # setup colorbar axes.
plt.colorbar(cax=cax) # draw colorbar
plt.savefig("AO_loading.png")
plt.show()

# save data for later use
np.save('eof1.npy', eof1.data)
np.save('pc1.npy', pc1.data)
np.save('mhgt.npy',mhgt[0:12,:,:].data)
np.save('pfactor1.npy',pfactor1)
np.save('lon.npy',lon)
np.save('lat.npy',lat)

できたプロット。NOAAのサイトにあるパターンと似たパターンを得ることができた。

このパターン（EOF第一モード）が説明する割合は19%。

d[0]/d.sum(axis=0)

0.19079784

2010/2/7 追記
緯度90度でうまく計算できないのは、計算精度により、cos90度　が負になり、その平方根をとっているからのようである。
2010/2/7 訂正
pfactor1はmasked arrayではないので、np.save('pfactor1.npy',pfactor1.data)からnp.save('pfactor1.npy',pfactor1)に修正。
2010/2/8 追記
緯度(lat)、経度(lon)も保存するようにした。
np.save('lon.npy',lon)
np.save('lat.npy',lat)

Scientific Peer Review, ca. 1945

2010-01-28T22:07:00.000+09:00

いろいろな気候インデックスとの相関

2010-01-28T17:20:00.005+09:00

*
*
NOAA　ESRLの
Linear Correlations in Atmospheric Seasonal/Monthly Averages
http://www.esrl.noaa.gov/psd/data/correlation/というページで各種気候インデックスといろいろな物理量の相関がプロットできる。

NINO3.4と表面温度の相関　(12月から2月）
日本は暖冬傾向だが、北日本でははっきりしない。
参考：気象庁

北極振動と表面温度の相関　(12月から2月）
日本は正の相関である。つまり北極振動が負だと寒冬の傾向がある。
参考：気象庁

北極振動　気圧場

2010-01-24T00:07:00.001+09:00

*
*
2008年と2009年の気圧場（1000hpa Geopotential height）を比較する。
2009年は北極振動のパターンと逆の符号（負のフェーズ）の傾向が強い。

以下プログラム。

を使用。
データはOpendapとい仕組みでNOAAからNCEP/NCARのデータを取得

from grads.ganum import GaNum # for python interface for grads

from mpl_toolkits.basemap import Basemap # for Basemap toolkit
import matplotlib.pyplot as plt
import numpy as np

for nyear in [2008,2009]:
# data import using python interface for grads
    ga = GaNum(Bin='grads')
    fh=ga.open("http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/pressure/hgt.mon.mean.nc") # data open
    ga("set time dec%(nyear)s" % locals())
    ga("set lat 0 90")
    ga("q file")
    ga("set lev 1000") # 1000 hpa
    hgt=ga.exp("hgt")
    del ga

    lat=hgt.grid.lat
    lon=hgt.grid.lon

# map
    fig=plt.figure()
    ax = fig.add_axes([0.1,0.1,0.8,0.8])

    map = Basemap(projection='nplaea',boundinglat=10,lon_0=135,resolution='i')
    map.drawcoastlines()
    map.drawmapboundary()
    lons,lats=np.meshgrid(lon[:],lat[:])
    x,y=map(lons,lats)
    map.fillcontinents()
    cz=map.contour(x,y,hgt[:,:],np.arange(-100,300,20))
    plt.clabel(cz,fmt="%3.0f")
    plt.title("Geopotential height (m): Dec %(nyear)s " % locals())
    plt.savefig("gpt_%(nyear)s.png" % locals())
    plt.show()

北極振動　Rで統計　グラフ編

2010-01-23T01:30:00.009+09:00

**

前回の続き

#boxplot
scripts="""
png(’boxplot.png’)
boxplot(rdatain)
dev.off()
"""
r(scripts)

去年の12月の値は「外れ値」にはなっていない

#hist
scripts="""
png('hist.png')
hist(rdatain,seq(-4, 4, 0.4),prob=TRUE)
x<-seq(-4, 4, 0.2)
lines(x, dnorm(x, mean=mean(rdatain), sd=sqrt(var(rdatain))), lty=3)
dev.off()
"""
a=r(scripts)

#cdf plot
scripts="""
png('cdf.png')
plot(ecdf(rdatain), do.points=FALSE, verticals=TRUE)
x<-seq(-4, 4, 0.01)
lines(x, pnorm(x, mean=mean(rdatain), sd=sqrt(var(rdatain))), lty=3)
dev.off()
"""
r(scripts)

#qqplot
scripts="""
png('qqplot.png')
qqnorm(rdatain)
qqline(rdatain)
dev.off()
"""
r(scripts)

正規分布より分布がせまい。

北極振動　Rで統計

2010-01-21T13:24:00.008+09:00

北極振動　12月だけを抜き出しの続き

12月のデータだけ抜き出して、Rで統計処理する。 Rpy2を使う。
まずは以下のように準備。readAOとpickup_monthは以前作ったもの。

AO_rpy2.py

from readAO import readAO # read AO data
from pick_month import pick_month # pick month
#
import scikits.timeseries as ts
#
import rpy2.robjects as robjects

# read AO data
AOindex_series=readAO()

# pickup December
AOindex_dec=pick_month(AOindex_series,month=12)

#remove missig value
data=AOindex_dec[AOindex_dec.mask==False].data

# for R
rdata=robjects.FloatVector(data)
robjects.globalEnv["rdatain"] = rdata
r=robjects.r

それで、

a=r('summary(rdatain)')

または

a=r.summary(rdata)

または

summary=r['summary']
a=summary(rdata)

いずれも同じ結果

print a
Min. 1st Qu. Median Mean 3rd Qu. Max.
-3.41300 -1.24200 -0.08762 -0.19740 0.82490 2.28200

print a[0]
-3.413

print a.names
[1] "Min." "1st Qu." "Median" "Mean" "3rd Qu." "Max."

b=a.r["Mean"]
print b
Mean
-0.1974

b[0]
-0.19739999999999999

print a.subset(1)
Min.
-3.413

参考
Rで統計: データ集合中の最大、最小、平均、中央値 - summary()関数http://www.yukun.info/blog/2008/09/r-summary-mean-median.html

北極振動　12月だけを抜き出し

2010-01-19T22:39:00.001+09:00

続き
北極振動
北極振動　年平均
AO 各年ごとにプロット　および　気候学平均

北極振動の時系列から12月だけを抜き出す。

以下、プログラム。
pick_month.py

from readAO import readAO
#
import numpy as np
import matplotlib.pyplot as plt
import scikits.timeseries as ts
import scikits.timeseries.lib.plotlib as tpl

def pick_month(series,month):
    """ pick time series only month==month
    """
    mask=(series.dates.month==month)
    out=series[mask].asfreq("A")
    return out

if __name__ == '__main__':

    AOindex_series=readAO()

# pickup December
    AOindex_dec=pick_month(AOindex_series,month=12)

#plot
    fig = tpl.tsfigure()
    fsp = fig.add_tsplot(111)
    fsp.tsplot(AOindex_dec, 'r-')
    plt.xlabel("YEAR")
    plt.ylabel("AO index")
    plt.savefig("AOindex_dec.png")
    plt.show()

AO 各年ごとにプロット　および　気候学平均

2010-01-19T00:39:00.002+09:00

続き
北極振動
北極振動　年平均

AOを各年毎にプロット（黒線）
2009年だけ赤線
気候学的平均（各月ごと、青線）

2009年の12月は異常に値が小さかったことがわかる。

以下プログラム。

plot_AO_season.py
from readAO import readAO
#
import numpy as np
import matplotlib.pyplot as plt
import scikits.timeseries as ts
import scikits.timeseries.lib.plotlib as tpl

AOindex_series=readAO()
# to seasonal
AOindex_season=ts.extras.convert_to_annual(AOindex_series)

#plot each year
fig = tpl.tsfigure()
fsp = fig.add_tsplot(111)

for n in range(AOindex_season.shape[0]):
    series = ts.time_series(AOindex_season[n,:], start_date=ts.Date(freq="m",year=1990,month=1),length=12)
    lstyle='k-'
    if AOindex_season.dates.year[n] == 2009:
        lstyle='r-'
    fsp.tsplot(series, lstyle)

# climatological average

AOindex_average = ts.time_series(np.mean(AOindex_season,axis=0), start_date=ts.Date(freq="m",year=1990,month=1),length=12)
fsp.tsplot(AOindex_average, 'b-',linewidth=2)
fsp.xaxis.set_ticklabels(["Jan"]) # dirty trick for plot
plt.savefig("AO_season.png")
plt.show()

Sageで数学　素数

2010-01-18T22:00:00.054+09:00

Sageで素数を遊んでみる。

sage: A=Primes()
sage: A.first()
2
sage: A.next(2)
3
sage: A.next(3)
5
sage: A.next(5)
7
sage: A.next(4)
5
sage: A.unrank(0)
2
sage: A.unrank(4)　#Returns the 4-th prime number.
11

sage: A.unrank(100)

547
sage: 17 in A
True
sage: A.cardinality() # Counts the elements of the enumerated set.
+Infinity

sage: prime_pi(548) # Return the number of primes <=549
101
sage: plot(prime_pi,1,600)

2010/9/17 追記　素数定理
x 以下の素数の数（下、赤線）は

（緑線）に近い。
あるいは、その近似である

（青線）に近い。

a=plot(prime_pi,2,10000,color="red")
b=plot(Li,2,10000,color="green")
c=plot(x/ln(x),2,10000,color="blue")
a+b+c

参考: Wikipedia 素数定理

算数の冒険(川端裕人著）

北極振動　年平均

2010-01-13T11:35:00.000+09:00

続き
北極振動

年平均をつくる。

前回のreadAO.pyを用いる。

from readAO import readAO
#
import numpy as np
import matplotlib.pyplot as plt
import scikits.timeseries as ts
import scikits.timeseries.lib.plotlib as tpl

AOindex_series=readAO()
# to annual mean
AOindex_annual=ts.convert(AOindex_series,'A',func=np.ma.mean)

#plot
fig = tpl.tsfigure()
fsp = fig.add_tsplot(111)
fsp.tsplot(AOindex_annual, 'r-+')
istart=fsp.get_xlim()[0]
iend=ts.Date('A','2010').value
fsp.set_xlim(istart,iend)
plt.xlabel("YEAR")
plt.ylabel("AO index")
plt.savefig("AOindex_anu.png")
plt.show()

北極振動

2010-01-12T16:02:00.004+09:00

最近、寒いのは北極振動が負になっているのが関係しているらしい。
北極振動　wikipedia
http://ja.wikipedia.org/wiki/%E5%8C%97%E6%A5%B5%E6%8C%AF%E5%8B%95
http://en.wikipedia.org/wiki/Arctic_oscillation

AOのデータは以下から入手できる。
http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/ao.shtml

プロットしてみる

以下プログラム(python)。scikits.timeseriesを使っている。

readAO.py

def readAO():
    "import monthly AO data"

    import numpy as np
    import scikits.timeseries as ts

# data is downloaded from http://www.cpc.ncep.noaa.gov/products/precip/CWlink/daily_ao_index/monthly.ao.index.b50.current.ascii
    f=open("monthly.ao.index.b50.current.ascii",'r')
    data=np.loadtxt(f)
    f.close()

    AOindex=data[:,2]
    year=data[:,0]
    month=data[:,1]

# make time series
    tbeginmonth=str(int(year[0]))+"-"+str(int(month[0]))
    dbeginmonth=ts.Date('M', tbeginmonth)
    AOindex_series=ts.time_series(AOindex,start_date=dbeginmonth,freq='M',mask=(AOindex<-99.))

    return AOindex_series
##################################

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

    AOindex_series=readAO()

    fig = tpl.tsfigure()
    fsp = fig.add_tsplot(111)
    fsp.tsplot(AOindex_series, '-')

    istart=fsp.get_xlim()[0]
    iend=ts.Date('M','2010-1').value
    fsp.set_xlim(istart,iend)
    plt.xlabel("YEAR")
    plt.ylabel("AO index")
    plt.savefig("AOindex.png")
    plt.show()