2011年10月17日月曜日

フィルタの応答関数



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

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

2011年5月6日金曜日

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



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

以下がそのモジュール。サンプルのメインプログラムでは、2011年1月1日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)



2011年3月5日土曜日

気象庁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.f90
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

2011年2月25日金曜日

Gradsでregrid



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

OpenGrADSre関数を用いた。
下のソース例、
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)という二つの1月の海面水温データセットから2011年1月の値をそれぞれプロットし、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"




2011年2月18日金曜日

GaLab in Python Interface to GrADS



Python Interface to GrADSGaLabを用いてプロットしてみた。
プロットしたデータはEarth System Research Laboratoryで配布されているNCEP/NCAR再解析の2011年1月の海水面気圧を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よりも手軽にプロットできる半面、緯度・経度線を取り除いたりとかの自由はききにくいようだ。