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