富士山の初冠雪日(1) データ取得

富士山　一気に冬の顔
初冠雪、平年より６日早く
山頂から８合目付近にかけてうっすらと雪をかぶった富士山＝富士吉田市上吉田 甲府地方気象台は２５日、富士山が初冠雪したと発表した。平年より６日、昨年より１２日早い。独自に観測している富士吉田市も同日、昨年より１５日早く「富士山初雪化粧」を宣言した。山梨県内は２２日に甲府など３地点で最高気温が３５度以上の猛暑日を記録したばかり。富士山は早くも冬に向けての“衣替え”が始まった。
2010年09月26日 山梨日日新聞　web魚拓

甲府気象台のwebページ表から、富士山初冠雪日を取得するプログラムを作る。
以下がプログラム。
関数 read_FirstSnowFuji()が初冠雪日を取得するモジュールで、年、月、日をnumpyの配列として返す。HTMLから表を読むのにtable parserを用いている。
メインプログラムでは、日付をscikits.timeseriesで年初からの日にち（1月１日を１とする）に変換し、グラフにしている。平均値を黒い実線で、平均値から標準偏差（約14日）の2倍離れた値に黒点線を引いた。ばらつきがかなり大きい印象だ。今年は少しだけ平年より早い。

read_FirstSnowFuji.py

import numpy as np
# table parser from http://simbot.wordpress.com/2006/05/17/html-table-parser-using-python/
from table_parser import *

def read_FirstSnowFuji():
    import urllib
    """
   read first snow cover day of the FUJI mountain from KOFU local meteorological observatory
   usage:
      year,month,day=read_FirstSnowFuji()
   input
      none
   output
      year,month,day:  first snow cover day (numpy array)
   
   """

    f = urllib.urlopen('http://www.jma-net.go.jp/kofu/menu/siryo_hatsukansetsu.dat.html')
    p = TableParser()
    p.feed(f.read())
    f.close()
    data=p.doc # Get to the data

    #get 
    year=[[],[],[]] # 3 column
    month=[[],[],[]] # 3 column
    day=[[],[],[]] # 3 column
    for d in data[0][1:]:
        for n,c in enumerate([0,4,8]):
            if not d[c]=="":
                year[n].append(int(d[c]))
                mi,di=d[c+1].split("/")
                month[n].append(int(mi))
                day[n].append(int(di))
    # to numpy array
    ayear=np.array([],dtype="i4")
    amonth=np.array([],dtype="i4")
    aday=np.array([],dtype="i4")
    for n in range(3):
        ayear=np.append(ayear,year[n])
        amonth=np.append(amonth,month[n])
        aday=np.append(aday,day[n])

    return ayear,amonth,aday

############################################
if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import scikits.timeseries as ts

    ayear,amonth,aday= read_FirstSnowFuji()

    days=np.array([],dtype="i4")
    for  n in range(ayear.size):
        td=ts.Date("D",year=ayear[n],month=amonth[n],day=aday[n])
        days=np.append(days,td.day_of_year)

    #plot
    plt.plot(ayear,days,"b-")
    mean=days.mean()
    stdp=days.mean()+2*days.std()
    stdm=days.mean()-2*days.std()
    plt.plot([ayear[0],ayear[-1]],[mean,mean],"k-")
    plt.plot([ayear[0],ayear[-1]],[stdp,stdp],"k--")
    plt.plot([ayear[0],ayear[-1]],[stdm,stdm],"k--")
    plt.xlim(ayear[0],ayear[-1])
    plt.xlabel("year")
    plt.ylabel("day of year")
    plt.savefig("FujiFirstSnow.png")
    plt.show()