pythonでssh

pythonで外部のサーバーにsshでアクセスして、そこでコマンドを実行する。
そのためにparamikoというパッケージを用いる。

以下はその例。
外部サーバにアクセス(ホスト名***, user名xxx, password +++）にアクセスし、コマンド"du . | sort -n -r"（どのディレクトリが容量を使っているかを調べる）を実行して、標準出力に表示する。

#!/usr/bin/env python

import paramiko

host="***"
port=22
user="xxx"
passwd="+++"

if __name__ == "__main__":
    paramiko.util.log_to_file('paramiko.log')
    s=paramiko.SSHClient()
    s.load_system_host_keys()
    s.connect(host,port,user,passwd)
    stdi,stdo,stde = s.exec_command("du . | sort -n -r")
    print stdo.read()
    s.close()

参考文献
Python ポケットリファレンス (Pocket Reference)　３１７ページ

(追記）
私がparamikoを使っている時、時々timeoutの問題が起こることがあった。
その場合
http://www.stillhq.com/commentform.cgi?post=python/paramiko/000004
などを参考にして、回避する。

ENDNOTE: MS Wordの文献リストのhanging indent

MS WORDでENDNOTEを使ってHanging indentをするやり方がずっとわからなかったのだが、やっとわかった。

Hanging indentというのはこういうの

設定はWordのEndnoteタブのBibliographyのLayoutでHanging indentの大きさを指定するとともに

Endnote本体のStyle ManagerのEditorで、BibliographyのLayoutの右隅（これに気が付かなかった）のHanging indentで"All paragraphs"など選んでやれば良い

黒潮の流路をpandasで

以前、気象庁黒潮の流軸の緯度のデータをネットから読み込んでプロットしてみた。

同じことをpandas (Python Data Analysis Library)で行ってみる。

import numpy as np
from  pandas import *
import matplotlib.pyplot as plt
ds = np.DataSource(None)
f=ds.open('http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/kuro_slat.txt')
data=read_csv(f,sep=" ",names=("space","val"),parse_dates=True)
f.close()
data["val"].plot()
plt.xlabel("YEAR")
plt.ylabel("Kurosho latitude")
plt.savefig("pandas_example.png")

(注記）　上のデータ読み込みでは、元データがスペースを使っている関係で、"space"という余計なデータ列(column)を作ってしまう。もっとスマートな方法があるかもしれない。
(追記:2012/7/2３)　後日、改良版を書いた。

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

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
これで図を作ったところ、従来横線が入っていた場合に、横線が入らなくなった。常にそうであるかは要確認。

フィルタの応答関数

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

以下の例は、京都大学の講義「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)

気象庁全球数値予報モデル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

気象庁の海面水温データ（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