2012年9月28日金曜日

Sageで方程式

Sageで200万が5年後に202万になる利息を計算する。

x=var(x)
q=solve(200*(1+x)**5==202,x)
for n in range(len(q)):
    print N(q[n].rhs())


結果

-0.690367429042489 + 0.952951066209182*I
-1.81062859479078 + 0.588956148532726*I
-1.81062859479078 - 0.588956148532726*I
-0.690367429042489 - 0.952951066209181*I
0.00199204766653335

したがって0.00199204766653335=約0.199%



2012年7月26日木曜日

NCEP Climate Forecast System Reanalysis




NCEP Climate Forecast System Reanalysisというデータがある。
http://cfs.ncep.noaa.gov/cfsr/

大気海洋モデル結合モデルに対してデータ同化が行われているのが特色だ。
分解能は従来のNCEPの再解析よりも高い。

大気: 水平分解能~38 km (T382)、 鉛直64 levels
海洋: 熱帯で0.25度、それ以外で0.5度の水平分解能、鉛直40 levels

また変数によっては1時間間隔のデータがある。
もっとも、データ同化が行われるのは6時間間隔だから、その間は予報データでうめてある。

このデータを見てみた。

UCARのサイトの
http://rda.ucar.edu/pub/cfsr.html
NCEP Climate Forecast System Reanalysis (CFSR) Selected Hourly Time-Series Products, January 1979 to December 2010http://rda.ucar.edu/datasets/ds093.1/
から、2001年の12月のSSTデータ(中身は海面下5mの温度)
ocnsst.gdas.200112.grb2
をダウンロードした。

これはgrib2形式なので、これをwgrib2でgrads形式に変換する。

g2ctl ocnsst.gdas.200112.grb2 > sst1hr.ctl
gribmap -i sst1hr.ctl

で、grads ctlファイルとindexファイル

sst1hr.ctl
ocnsst.gdas.200112.grb2.idx

ができる。水平分解能は0.5x0.5になっていた。

gradsで最初のステップをプロットすると、こんなかんじ


0N, 180Eの最初の48ステップをプロットすると
6時間ごとに値が飛んでいるように見えるのは、6時間ごとに同化(3D-VAR)で解析値に修正していることにともなうものであろう。

grads形式に変換する時に、解析値だけ取り出すには"-0"(ハイフンゼロ)をつけて

g2ctl -0 ocnsst.gdas.200112.grb2 > sst1hr.ctl
gribmap -0 -i sst1hr.ctl


この場合、0N, 180Eの最初の8ステップ(6時間ごとなので、48時間分)をプロットすると


2012年7月24日火曜日

Wvisで気象情報を3次元可視化




Wvisという、気象庁のメソモデル(MSM)を3次元可視化してくれるソフトがある。
今のところMSMしか可視化できないが、MSMなら手軽に3次元可視化できる。

Wvisは以下からダウンロードできる。
http://www.cybernet.co.jp/avs/download/wvis.html

Windowsのみのソフトだが、私の試した限り、Linux上でも、wineを用いることで問題なく使えた

Wvisについては、気象学会の「天気」にも解説がある。
http://www.metsoc.jp/tenki/pdf/2011/2011_09_0073.pdf

使い方については、Wvisのマニュアルまたは、「天気」の記事を参照されたい。

MSMのデータは
http://database.rish.kyoto-u.ac.jp/arch/jmadata/gpv-original.html
または
http://gpvjma.ccs.hpcc.jp/~gpvjma/
でgrib2形式のデータが入手可能である。

例として、下の動画は2012年7月11-14日にかけて九州に豪雨をもたらした暖かく湿った空気(相当温位355K)をプロットしたものである。説明についてはYoutubeのキャプションの欄を参照されたい。http://www.youtube.com/watch?v=GO0ZrA3wx5s


浦神と串本の日平均潮位差をpandasで



以前書いたpandasでの時系列処理の続き。バージョンアップ(これを書いてる時点でversion 0.8.1)により、さらに便利になった。

浦神と串本の日平均潮位差をプロットしてみる。前回のようにwebから直接読むこともできるが、今回はhttp://www.data.kishou.go.jp/db/kobe/kuroshio/chouisa/chouisa.datからchouisa.datをダウンロードしてそれを読むことにする。

下がプログラム。
read_psvでデータを読んでいる。
前回は日付とデータの間にあるスペースをスマートに処理できなかったが、sep='\s*' (任意の数のスペース)を使えばよかったらしい。
column=0(右端)をindexとして用い(index_col=0)、日付データとして解釈する(parse_dates=Trule)。
”-”を欠損値と解釈する。

version 0.8以降、月平均を計算するのが簡単になり、data.resample("M",how='mean')で良い。


import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# read data
data=pd.read_csv("chouisa.dat",sep='\s*',skiprows=1,header=None,index_col=0,parse_dates=True,names=("date","chouisa"),na_values="-")

#monthly value
data2=data.resample("M",how='mean')

#plot
plt.subplot(111)
data["chouisa"].plot(style='b-',label="Daily")
data2["chouisa"].plot(style='r-',label="Monthly",linewidth=2)
plt.legend(loc='best')
plt.savefig("choisa.png")
plt.show()



2012年4月28日土曜日

pythonでssh



pythonで外部のサーバーにsshでアクセスして、そこでコマンドを実行する。
そのためにparamikoというパッケージを用いる。

以下はその例。
外部サーバにアクセス(ホスト名***, user名xxx, password +++)にアクセスし、コマンド"du . | sort -n -r"(どのディレクトリが容量を使っているかを調べる)を実行して、標準出力に表示する。

#!/usr/bin/env python

import paramiko

host="***"
port=22
user="xxx"
passwd="+++"

if __name__ == "__main__":
    paramiko.util.log_to_file('paramiko.log')
    s=paramiko.SSHClient()
    s.load_system_host_keys()
    s.connect(host,port,user,passwd)
    stdi,stdo,stde = s.exec_command("du . | sort -n -r")
    print stdo.read()
    s.close()

参考文献
Python ポケットリファレンス (Pocket Reference) 317ページ

(追記)
私がparamikoを使っている時、時々timeoutの問題が起こることがあった。
その場合
http://www.stillhq.com/commentform.cgi?post=python/paramiko/000004
などを参考にして、回避する。

ENDNOTE: MS Wordの文献リストのhanging indent



MS WORDでENDNOTEを使ってHanging indentをするやり方がずっとわからなかったのだが、やっとわかった。

Hanging indentというのはこういうの


設定はWordのEndnoteタブのBibliographyのLayoutでHanging indentの大きさを指定するとともに

Endnote本体のStyle ManagerのEditorで、BibliographyのLayoutの右隅(これに気が付かなかった)のHanging indentで"All paragraphs"など選んでやれば良い





2012年3月16日金曜日

黒潮の流路をpandasで



以前、気象庁黒潮の流軸の緯度のデータをネットから読み込んでプロットしてみた。

同じことをpandas (Python Data Analysis Library)で行ってみる。

import numpy as np
from  pandas import *
import matplotlib.pyplot as plt
ds = np.DataSource(None)
f=ds.open('http://www.data.kishou.go.jp/kaiyou/shindan/b_2/kuroshio_stream/kuro_slat.txt')
data=read_csv(f,sep=" ",names=("space","val"),parse_dates=True)
f.close()
data["val"].plot()
plt.xlabel("YEAR")
plt.ylabel("Kurosho latitude")
plt.savefig("pandas_example.png")

(注記) 上のデータ読み込みでは、元データがスペースを使っている関係で、"space"という余計なデータ列(column)を作ってしまう。もっとスマートな方法があるかもしれない。
(追記:2012/7/23) 後日、改良版を書いた。

2012年2月6日月曜日

GrADSで作成するpostscriptの不要な横線

GrADSでpostscript形式の図を作ると不要な横線が入ってしまうことがある。

従来の対策は、例えば
http://www.sci.hokudai.ac.jp/~sasakiyo/tips_grads.html

GrADS 2.0.0から従来のgxout shaded に加えて、gxout shade2が加わっている。
http://www.iges.org/grads/gadoc/gradcomdsetgxout.html
これで図を作ったところ、従来横線が入っていた場合に、横線が入らなくなった。常にそうであるかは要確認。