Google Earthでアニメーション

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

手順１
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>

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

・GroundOverlayから/GraoundOverlayの部分を図の数だけ繰り返す。
・それぞれ図の名前を指定する。
・図がどの時間を占めるべきかTimeSpan /TimeSpanの部分で指定する。なお、今の場合、図は8月だけの図であるが、アニメーションで図が続く様に１年中同じ図が表示されるよう、１月１日から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で開くとこんな感じ。

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が出てくるようになった。
あとは再起動やアップグレードをしているうちに、ちゃんと正常に動くようになった。
良かった、良かった。

富士山の初冠雪日(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)が外れ値、すなわち１９１４年（ayear[21-1]=1914)と２００８年（ayear[115-1]=2008）が外れ値になるのは良いとして、weightが１に近い点は13だけで、あとの外れ値を除く１０２点はweightが１以下に落とされている。

これを見るために、回帰直線の係数、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()

ちなみに今年２０１０年の重みはweight[-1]=0.86999219949571671
このように多数の点の重みを落とすのは統計的に良くても、物理的な理由は自明ではない。
今回、この方法を用いたのは適当ではなかったかもしれない。

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

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

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

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

前回作ったモジュール(read_FirstSnowFuji)で初冠雪日を読み、年初からの日にち（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%でトレンドがない可能性を否定できない。データのばらつきが大きいことから、さもありなんといったところである。

富士山の初冠雪日(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()