2010年12月17日金曜日

気象庁全国合成レーダー降雨



気象庁の全国合成レーダー降雨データがこちらで配布されている。
データは10分ごとにある。
データは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時00分世界標準時(日本時間午後5時)の降水を東経125-135度北緯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")


2010年11月7日日曜日

Google Graphで棒グラフ



GE-GraphというGoogle Graph上に棒グラフをプロットできるフリーソフトがある。 http://www.sgrillo.net/googleearth/gegraph.htm
これを使ってみた。

結果から。下の図は正味海洋-大気熱フラックス(海から大気)の1月長期平均をGoogle Earth上で棒グラフにしたものである。正値のみ使っている。
このkmzファイルはこちら

GE-Graphに入力するデータとして
場所名,緯度,経度,値
で構成されるcsvファイルを用意する必要がある。
上の例のcsvファイルはコレ(その1,その2
GE-Graphが5000地点以上を受け付けないのでファイルを2つに分けている。場所名は特に意味は無いので全て0を入れている。

上のデータをGE-Graphに読み込ませる。設定は以下の図のようにした。
それぞれのデータで"Run"を行ってkmzファイルを作った。Google Earth上で二つのファイルをまとめて、saveしたのが上のkmzファイルである。

最後にcsvファイルを作るpythonプログラムを添付する。
NCEP/NCARのデータをESRLのサイトからPython Interface to GrADS を使って読み込む。
latent/sensible heat flux, longwave/shortwave radiationを全て足して正味の熱フラックスを計算する。
陸・海マスクを読み、海のデータだけを使う。また正(海から大気)のデータだけ残している。データは4000点づつセーブする。
from grads.ganum import GaNum  # for python interface for grads
import numpy as np

#====== DATA read
ga = GaNum(Bin='grads')
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/lhtfl.mon.ltm.nc") #net latent heat flux
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/shtfl.mon.ltm.nc") #net sensible heart flux
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/nswrs.mon.ltm.nc") #net shortwave radiation
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis.derived/surface_gauss/nlwrs.mon.ltm.nc") #net longwave radiation
ga("set t 1")  # January
ga("netheat=lhtfl.1+shtfl.2+ nswrs.3+nlwrs.4") # all net heat
netheat=ga.exp("netheat")
ga("close 4")
ga("close 3")
ga("close 2")
ga("close 1")
ga("sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/ncep.reanalysis/surface_gauss/land.sfc.gauss.nc")  #land sea mask
land=ga.exp("land")
del ga

lon=netheat.grid.lon
lat=netheat.grid.lat

#===== DATA manupilation
data=netheat.data
data[land>0]=-999.  # remove values on land

x,y=np.meshgrid(lon,lat)

vl=0.               # data criteria 

xout=x[data>vl].flatten()      # only  values >vl  are used
xout[xout>180]=xout[xout>180]-360. #   -180.=<lon<=180.  
yout=y[data>vl].flatten()      # only  values >vl  are used
zout=data[data>vl].flatten()   # only  values >vl  are used
points=np.zeros(data[data>vl].size) # position number (anythin is ok)

#==== DATA output
nn=zout.size/4000+1
for n in range(nn):
    st=4000*n
    en=4000*(n+1)
    np.savetxt("data%(n)s.txt" %locals(),np.transpose((points[st:en],yout[st:en],xout[st:en],zout[st:en])),fmt='%4i,%7.3f,%7.3f,%7.3f')


2010年10月29日金曜日

Google Earthでアニメーション



以前、Google Earthで図を描いたが、今回はアニメーションを作ってみる。図が3枚だけのアニメーションだ。

手順1
Gradsで図を3枚つくる。GradsでGoogle Earth向きの図をつくるやり方は以前の記事を参照のこと。下の例はNCEP/NCARの最解析データを使って、2008年から2010年の8月の表面大気温度のアノマリをプロットしている。colorbarはxcbarを使って図の上に一緒に乗せた。
"reinit"
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/air.mon.mean.nc"
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/air.mon.ltm.nc"
year=2008
while (year<=2010)
"set time aug"year
"set x 1 144"
"set y 1 73"
"set parea 0 11 0 8.5"
"set grid off"
"set grads off"
"set mproj scaled"
"set mpdraw off"
"set gxout shaded"
"color -5 5 0.5 -kind blue->gray->red"
"d air.1-air.2(t=8)"
"xcbar  4.5 6.5 4 4.2 -edge triangle -fstep 10 "
"printim  airtemp_"year"aug.png x2048 y1024 -t 0"
"c"
year=year+1
endwhile

こんな感じの図ができる。


手順2
基になるkmlファイルをGradsで作成する。やり方は以前の記事と同じだ。
*************KML
"reinit"
"sdfopen http://www.esrl.noaa.gov/psd/thredds/dodsC//Datasets/ncep.reanalysis.de
rived/surface/air.mon.mean.nc"
"set time aug2010"
"set x 1 144"
"set y 1 73"
"set kml airtemp"
"set gxout  kml"
"d air"

以下のようなkmlファイルが出来る。
<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
  <GroundOverlay>
    <name>Monthly Mean Air Temperature</name>
    <Icon>
      <href>airtemp.tif</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>
</kml>

手順3
kmlファイルを以下の様に書き換える。
ポイントは

GroundOverlayから/GraoundOverlayの部分を図の数だけ繰り返す。
・それぞれ図の名前を指定する。
・図がどの時間を占めるべきかTimeSpan /TimeSpanの部分で指定する。なお、今の場合、図は8月だけの図であるが、アニメーションで図が続く様に1年中同じ図が表示されるよう、1月1日から12月31日をTimeSpanで設定している。
・繰り返しの全体をFolder/Folderではさむ。 

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2">
<Folder>
  <GroundOverlay>
    <name>Monthly Longterm Mean of Air Temperature</name>
    <TimeSpan>
      <begin>2008-01-01</begin>
      <end>  2008-12-31</end>
    </TimeSpan>
    <Icon>
      <href>airtemp_2008aug.png</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>

  <GroundOverlay>
    <name>Monthly Longterm Mean of Air Temperature</name>
    <TimeSpan>
      <begin>2009-01-01</begin>
      <end>  2009-12-31</end>
    </TimeSpan>
    <Icon>
      <href>airtemp_2009aug.png</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>

  <GroundOverlay>
    <name>Monthly Longterm Mean of Air Temperature</name>
    <TimeSpan>
      <begin>2010-01-01</begin>
      <end>  2010-12-31</end>
    </TimeSpan>
    <Icon>
      <href>airtemp_2010aug.png</href>
    </Icon>
    <LatLonBox>
      <west>     -1.25</west>
      <east>    358.75</east>
      <south>    -91.25</south>
      <north>     91.25</north>
      <rotation>0.0</rotation>
    </LatLonBox>
  </GroundOverlay>
</Folder>
</kml>

出来たファイルをkmzにまとめたものをこちらに置いた。kmzファイルはkmlと図をまとめたもので、Google Earthで一度kmlファイルを開いて、「名前をつけて保存」で作成できる。

Google Earthで開くとこんな感じ。


2010年10月25日月曜日

Ubuntu 10.10にアップグレードメモ



Ubuntuを10.10にアップグレードした。
順調にアップグレードしてくれると思ったらあまかった。
再起動後GUIの画面が始まらず、CUIでログインできるだけ。
リカバリーモードさえうまくいかない。
その夜はそれ以上進まず、諦めて寝ることにした。
次の日ネットで調べたら以下の情報を見つけた。

After upgrade to 10.10 no GUI just command line
http://ubuntuforums.org/showthread.php?t=1575262

これに書いてある中で、 fglrxを消すというのが効果があった。
これで、リカバリーモードのSafe Graphics ModeでGUIが出てくるようになった。
あとは再起動やアップグレードをしているうちに、ちゃんと正常に動くようになった。
良かった、良かった。

2010年10月7日木曜日

富士山の初冠雪日(2) 続・トレンド



前回の続きとして、私自身あまりよく理解しているとは言えないが、ココを参考にして、ロバスト推定でも線形回帰をとってみた。

参考:
ロバスト推定(統計数理研究所)
http://www.ism.ac.jp/~fujisawa/research/robust.html
ロバスト推定(朱鷺の社)
http://ibisforest.org/index.php?%E3%83%AD%E3%83%90%E3%82%B9%E3%83%88%E6%8E%A8%E5%AE%9A

Rのパッケージrobustbaseを使うために
install.packages("robustbase")
を行っておく。
以下がプログラム。
import numpy as np
import scikits.timeseries as ts
#module
from read_FirstSnowFuji import *

#for R
import rpy2.robjects as robjects

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)


#for R
ryear=robjects.FloatVector(ayear)
rdays=robjects.FloatVector(days)
robjects.globalEnv["ryear"] = ryear
robjects.globalEnv["rdays"] = rdays

r=robjects.r

scripts="""
library("robustbase")   #import robustbase
png('FujiSnow_trend_robust.png')
plot(ryear,rdays,xlab="Year",ylab="Day of year") #scatter plot 
fm<-lmrob(rdays ~ ryear)  # robust linear model
abline(fm)  #plot regression line
abline(lm(rdays ~ ryear),lty=2) # usual linear model
dev.off()
"""
r(scripts)

# summary
summary=r('sm<-summary(fm)')
print summary


以下が図。点線が前回の回帰直線、実線がrobustbaseで得られた直線である。


出力(print summary)は


Call:
lmrob(formula = rdays ~ ryear)

Weighted Residuals:
    Min      1Q  Median      3Q     Max
-57.989  -7.872  -0.125   5.986  25.291

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept) 76.10515   51.85081   1.468 0.144896   
ryear        0.10154    0.02673   3.798 0.000234 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Robust residual standard error: 10.03
Convergence in 12 IRWLS iterations

Robustness weights:
 2 observations c(21,115) are outliers with |weight| <= 0.00057 ( < 0.00085);
 13 weights are ~= 1. The remaining 102 ones are summarized as
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
0.00365 0.88600 0.94480 0.88410 0.98680 0.99880
Algorithmic parameters:
tuning.chi         bb tuning.psi refine.tol    rel.tol
 1.5476400  0.5000000  4.6850610  0.0000001  0.0000001
 nResample     max.it     groups    n.group   best.r.s   k.fast.s      k.max
       500         50          5        400          2          1        200
 trace.lev compute.rd
         0          0
seed : int(0) 

これによると回帰係数は0.10154(100年で約10日冠雪日が遅くなる)で、しかもp値が非常に小さく統計的に有意である。
ただし、c(21,115)が外れ値、すなわち1914年(ayear[21-1]=1914)と2008年(ayear[115-1]=2008)が外れ値になるのは良いとして、weightが1に近い点は13だけで、あとの外れ値を除く102点はweightが1以下に落とされている。

これを見るために、回帰直線の係数、weightなどの数値をpythonに戻し、プロットしてみる(同じことをRでもできるはずだが、pythonのmatplotlibのほうがなれているので)。
以下プログラム。
robust推定で得られた回帰直線を緑線で、
weightが1に近い点を青○で、
weightが0.5から0.9988の間にある点を黒点、
weightが0.00085から0.5の間にある点を赤点、
外れ値を赤xでプロットした。
#weights
weights=r('weights<-sm$weights')
weight=np.array(weights)

#plot
import matplotlib.pyplot as plt
plt.plot(ayear,ayear*a+b,'g-')
range=(weight>0.9988)
plt.plot(ayear[range],days[range],"bo")
range=np.logical_and(weight<=0.9988,weight>=0.5)
plt.plot(ayear[range],days[range],"k.")
range=np.logical_and(weight<0.5,weight>=0.00085)
plt.plot(ayear[range],days[range],"r.")
range=(weight<0.00085)
plt.plot(ayear[range],days[range],"rx")
plt.xlim(ayear[0],ayear[-1]) 
plt.xlabel("year")
plt.ylabel("day of year")
plt.savefig("FujiSnow_trend_robust_2.png")
plt.show()
ちなみに今年2010年の重みはweight[-1]=0.86999219949571671
このように多数の点の重みを落とすのは統計的に良くても、物理的な理由は自明ではない。
今回、この方法を用いたのは適当ではなかったかもしれない。

前回と今回の解析で、初冠雪日にトレンドがあるかは統計的にははっきりしなかったが、初冠雪日が遅くなっていっているかもしれないという作業仮説自体は興味深いものがある。これから、将来、どのような傾向があらわれてくるだろうか?

参考
異常気象を追う 初冠雪
http://www.bioweather.net/column/essay2/aw15.htm

縮小する富士山の永久凍土
http://www.fujisan-net.jp/data/article/1583.html

2010年10月6日水曜日

富士山の初冠雪日(2) トレンド



前回作ったモジュール(read_FirstSnowFuji)で初冠雪日を読み、年初からの日にち(1月1日を1とする)に変換したうえで、線形回帰を当てはめてみる。
以前やったように、rpy2を用いて、Rにデータを渡してやる。

import numpy as np
import matplotlib.pyplot as plt
import scikits.timeseries as ts
#module
from read_FirstSnowFuji import *

#for R
import rpy2.robjects as robjects

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)


#for R
ryear=robjects.FloatVector(ayear)
rdays=robjects.FloatVector(days)
robjects.globalEnv["ryear"] = ryear
robjects.globalEnv["rdays"] = rdays
r=robjects.r

scripts="""
png('FujiSnow_trend.png')
plot(ryear,rdays,xlab="Year",ylab="Day of year") # scatter plot
fm<-lm(rdays ~ ryear) # linear model
abline(fm) # plot regression line                 
dev.off()             
"""
r(scripts)
summary=r('sm<-summary(fm)')
print summary
出力(print summary)は
Call:
lm(formula = rdays ~ ryear)

Residuals:
     Min       1Q   Median       3Q      Max
-54.2679  -5.8555   0.9491   8.1001  27.3063

Coefficients:
             Estimate Std. Error t value Pr(>|t|) 
(Intercept) 138.25153   71.89542   1.923   0.0570 .
ryear         0.06873    0.03683   1.866   0.0645 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 13.45 on 115 degrees of freedom
Multiple R-squared: 0.0294,    Adjusted R-squared: 0.02096
F-statistic: 3.484 on 1 and 115 DF,  p-value: 0.06453

係数が0.06873/yearだから、100年に一週間ほどのトレンド(初冠雪日が遅くなる)となる。ただし、p値が約0.0645であるから危険率5%でトレンドがない可能性を否定できない。データのばらつきが大きいことから、さもありなんといったところである。

2010年10月1日金曜日

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



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

甲府気象台のwebページ表から、富士山初冠雪日を取得するプログラムを作る。
以下がプログラム。
関数 read_FirstSnowFuji()が初冠雪日を取得するモジュールで、年、月、日をnumpyの配列として返す。HTMLから表を読むのにtable parserを用いている。
メインプログラムでは、日付をscikits.timeseriesで年初からの日にち(1月1日を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()