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
これで図を作ったところ、従来横線が入っていた場合に、横線が入らなくなった。常にそうであるかは要確認。
2012年2月6日月曜日
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ヶ月移動平均にした場合。
参考） Scipy Cookbook FIRFilter
2011年5月6日金曜日
気象庁全球数値予報モデル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)
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.ｆ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
2011年2月25日金曜日
Gradsでregrid
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"
2011年2月18日金曜日
GaLab in Python Interface to GrADS
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")
2010年12月17日金曜日
気象庁全国合成レーダー降雨
気象庁の全国合成レーダー降雨データがこちらで配布されている。
データは１０分ごとにある。
データは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")
登録:
投稿 (Atom)