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

0 件のコメント:
コメントを投稿