Wavelet解析

Torrence and Compo [1998]  のwavelet解析手法と、そのツールは我々の分野でよく使われる。私もMiyama and Miyazawa [2014]をで使用している。その関連した発表をAOGS 2014でした時に、Torrence and Compo [1998] はスケールが大きいほうを過大評価するバイアスがあり、それに関する論文（Liu et al. 2007と関連するweb page)を教えていただいた。

下の図はLiu et al. 2007のFig.2に対応する、1,8,32,128,365日のsineカーブの単純な重ね合わせ(a)のシグナルにwavelet解析をかけたものである。Torrence and Compo [1998] の手法(b,c)ではスケールが大きい（長周期）のほうがシグナルが大きいと解析されてしまう。一方、Liu et al. [2007]はスペクトルをスケールで割ることを提案しており、これであれば(d,e)、各周期が同じくらいの強さであるという合理的な結果が出る。





上記の図を作るwavelet解析のツールをpythonに翻訳し、IPython Notebookにしたものはこちら。
http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/u0wb931fs4qfrpn/wavelet_test_sine.ipynb

さらに、NINO SST3のwavelet解析した(Liu et al. 2007のFig 4に対応）物をIPython Notebookしたものはこちら。
http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/400j051n0sustcy/wavelet_test_ElNino3_Liu.ipynb

Surface velocity data of MOVE/MRI.com at NEAR-GOOS RRTDB

NEAR-GOOSのRegional Real Time Databaseのサイトは、最近更新されて、データのダウンロードがregistration無しでダウンロードできるようになった（自動ダウンロードを行うにはregistrationが必要）。

ここでは、その中から、気象庁(JMA)のDaily Surface Currents in the seas adjacent to Japanをプロットして見る。

2014年11月6日のデータ "move_subc_jpn_D20141106.txt"をダウンロードして置き、CartopyでプロットするIPython Notebookはこちら。
https://www.dropbox.com/s/bpow09n6x4rm9sy/MOVE_velocity.ipynb?dl=0

Colors (diverging data)

前回はsequential data用のカラーマップについてテストしたが、今回は正負のアノマリなどをプロットする diverging カラーマップについてテストした。

ipython notebookはこちら
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/t8um9tdmauu8qpv/divergent_color.ipynb?dl=0

テストデータのしてはNOAA OISST　version 2の海面温度アノマリを使った。データ読み込みとプロットにはIrisを使った。

Irisはイギリスのソフトのためか、日付変更線をまたぐ場合には少し癖があり、そのデモでもある。

Colors (sequential data)

最近、カラーマップに関する考察をいくつか読んだ。

• THE RAINBOW IS DEAD…LONG LIVE THE RAINBOW!
   http://mycarta.wordpress.com/2012/05/29/the-rainbow-is-dead-long-live-the-rainbow-series-outline/
• Subtleties of Color
  http://earthobservatory.nasa.gov/blogs/elegantfigures/2013/08/05/subtleties-of-color-part-1-of-6/
• Stauffer, R., G. J. Mayr, M. Dabernig, and A. Zeileis, 2014: Somewhere over the rainbow: How to make effective use of colors in meteorological visualizations. Bulletin of the American Meteorological Society, 140710055335002. doi:10.1175/bams-d-13-00155.1

いずれも"rainbow"（またの名を"jet")と呼ばれるカラーマップの弊害を述べている。

それに対する代替はいくつか提案されている。

以下のipython notebookでは、それらを試してみた。
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/5t9jrhr7va3wscj/colors.ipynb?dl=0

データはGHRSST OSTIA SSTという海面温度で、2014年10月15日の117-150E, 17-40Nの範囲をIRISを使ってIPRC APDRCから読み込み、プロットさせている。

下は"cube1"と呼ばれるカラーマップでプロットした例。

ちなみに、この日の海面温度には台風19号と台風20号の通過に伴い、海水の低下した場所が見られる。

Reference

Kuroshio Oyashio Watch 2014/10/20

http://www.jamstec.go.jp/jcope/htdocs/e/kow/141020/KOWatch_20141020.html

等値線の位置を求める

Qiu and Chen (2005)などでは、海面高度データのある等値線の位置を黒潮続流の位置とし、Decadal Variationの議論している。

等値線の位置を決めるのは、渦などの紛らわしい等値線があり、目でみていちいち確認しないで自動的に欲しい等値線の位置を決めるのは難しい。

ここでは、図作成ソフトmatplotlibのcontourルチーンから情報を取り出して、等値線の位置を決める工夫をする。

参考にしたのはstackoverflowの

• How can I get the (x,y) values of the line that is ploted by a contour plot (matplotlib)? http://stackoverflow.com/questions/1560424/how-can-i-get-the-x-y-values-of-the-line-that-is-ploted-by-a-contour-plot-mat  
• matplotlib - extracting data from contour lines http://stackoverflow.com/questions/5666056/matplotlib-extracting-data-from-contour-lines  
• Python: find contour lines from matplotlib.pyplot.contour() http://stackoverflow.com/questions/18304722/python-find-contour-lines-from-matplotlib-pyplot-contour/18309914#18309914

ipythonノートブックは
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/qt804lsufj4whyx/Kuroshio_extension_path.ipynb?dl=0

Step1 力学高度のデータをferretを使い、読み込み、図示する

データはAVISOの絶対力学高度の2009年1月21日のデータを用いる。データはAVISOからダウンロードしておく。

データの読み込みと図示には、前回紹介したferret magicを使う。

このデータでは、0.8mの等値線が黒潮続流位置としてはよさそうである。本当はheat fluxによるsteric heightの変動を考慮にいれないとならないが(Qiu and Chen, 2005)、ここではその手順は省略。
Qiu and Chen (2005)では東経141から153度のデータで議論しているので、その範囲だけ抜き出しておく。

等値線は、渦や蛇行のなどのため等値線が４つあるが、黒潮続流にあたる本来欲しい線（東西にのびる線）だけを取り出すのが課題である。

Step 2 Ferretからpythonにデータを読み込む

matlabの機能を使うために、ferretからpythonへとデータを読み込む。これもferret magicを用いる。データがちゃんと読めているかを確認もしておく。

Step ３ 等値線の位置を求める

matlabのcontour routineの機能を使って、等値線の位置を求める。等値線の位置は４つ求まっている。

Step 4 等値線の位置情報から、もっとも長い等値線を選ぶ。

等値線の位置が数値的に求まったので、この中から、欲しいpathだけを選ぶ。最も長いpathが黒潮続流にあたる欲しいpathだと考えられる（渦などを避けるために、始点と終点がもっとも離れているpathという条件でも良いかもしれない。）。

まず、pathの長さを求める関数を作成する。

get_distance　は２点間の距離を km で求める関数。

get_pathlength　はget_distanceを用いてpath位置の情報の配列からpathの長さを求める関数。

上の４つのpathの長さを求め、path #2が最も長いことがわかる(約１432km)。

path #= 0 length= 415.785621835 km
path #= 1 length= 527.959783762 km
path #= 2 length= 1432.08536051 km
path #= 3 length= 563.215908518 km

これが、求めるpathである。

Reference:

Qiu, B., and S. Chen, 2005: Variability of the Kuroshio Extension Jet, Recirculation Gyre, and Mesoscale Eddies on Decadal Time Scales. Journal of Physical Oceanography, 35, 2090-2103. doi:10.1175/JPO2807.1

pyferret, ipython ferret magic

pyferretは、大気海洋分野のお絵描きソフトferretの、pythonとの結合を深めたバージョンである。

さらにipython-ferretmagicはipythonの中でferretを使うのを可能にしてくれる。

サンプルのノートブックはこちら。
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/avefneoyxdxf7ov/ferret_magic_example.ipynb?dl=0

この中で、Levitusのデータを読み込んで、9月の温度26度の等温線を求めてプロットしている。
ferretからpythonにデータを読み込むことも可能である。

surface drifter

NOAA, Atlantic Oceanographic and Meteorological Laboratory (AOML)のGlobal Drifter Program (GDP), Drifter Data Assembly Center (DAC)には表層漂流ブイのアーカイブがある。このデータを読み込んで図示してみる。

サイトからftpでデータをあらかじめダウンロードしておく。このうち、buoydata_10001_mar14.dat　(ブイ1万個目以降から2014年3月までのデータ）を例として扱う。

IPython Notebookはこちら。
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/uzv1bktgazcf821/Test_buoys.ipynb?dl=0

まず、pandasでデータを読み込む。データが多量なので、読み込むまで時間がかかるかもしれない。

データ(data)の内で、一部に緯度が北緯25から45度、東経120から160度が含まれるデータを選びだす(selecteddata)。そのうちで、ユニークなid番号を持つもののは851個ある。

あらためて、もとのデータ(data)から、そのid番号を持つデータを抜き出し(wanteddata)、pickle形式で保存する。

次に、そのデータの中から、例として最初の10個だけcartopyを用いてプロットした。

参考文献

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

IrisとETOPO1

Irisはイギリスで開発されているpythonベースの気象・海洋データ解析パッケージである。

Irisを使って地形データETOPO1を読み込んで、プロットしたり、補間してみた。

ETOPO1は1/60°の海底地形で、NOAA NGDCのサイトからETOPO1 Bedrock（氷無しの地形）、cell-registered（グリッドセルの位置で定義される）の中から、netcdfファイル　ETOPO1_Bed_c_gmt4.grdをダウンロードした。

IPython Notebookはこちらである。
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/67khdbgzc2bf0sp/Read_ETOPO1.ipynb

下の図はプロットの一例。

PDO index

Pacific Decadal Oscillation (PDO) index をJISAOのPDOページからデータを取得し、ある月から最新のデータまで、時系列をプロットする。

ipython notobookはこちら
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/y1lol5bl7h4xqdi/PDO.ipynb

bokehを使って、串本浦神の潮位差をプロット

串本と浦神の潮位差のデータをbokehを使ってプロットしてみる。bokehはpythonのプロットライブラリの一つで、ズームなどが可能なインタラクティブなプロットをweb上につくることができる。（bokehは日本語のボケから来てるらしい。何故そんな名前にしたんだか。）

ipython notebookはこちら。
http://nbviewer.ipython.org/urls/dl.dropbox.com/s/rry6nt3xs49o6s6/KushimotoUragami.ipynb

ズームなどができるのを確認してみてください。

データは以前のようにpandasで読み込んでいる。
daily dataの黒線で、30日の移動平均を赤線でプロットしている。ここでの30日平均は基準日を中心とした30日平均ではなく、基準日をさかのぼっての30日平均であることに注意。

Numba

Pythonでforをまわすような数値計算を行うとmatlabなどと同じで非常に遅い。
それを克服するためにpypyやcythonなどがある。しかし基本的にコードを書き換える必要があるので、もっと楽をしたい。

最近numbaというものを知った。ここのサンプルにある例だと、

from numba import autojit

@autojit

のおまじないだけで関数が劇的に早くなる。

少し、自分で例で実験してみた。
全球のSSTに対して、NINO3の相関をとるということをやってみる（これはnumbaの実験のために作った例であり、相関をとるだけならもっとエレガントな方法があるだろう）。

そのノートブックが
http://nbviewer.ipython.org/5396344

説明
(1) 図をプロットするためのおまじない。

(2) SSTアノマリとNINO3を読む関数を定義。
　　元データは　NOAA　OISST v2でmonthly dataのsst.mnmean.ncを使う。
　　データはPython Interface to GrADSを使って読む。
　　1982から2010年のデータから気候値からのanomaly（Bin Guan's GrADS Script Libraryのdeseasonを使った。古いバージョンを使ったので、今は使い方が違うかも知れない）とNINO3 SSTアノマリを定義し、返す。
　ここらへんは以前の記事も参照。
　/usr/local/lib/python2.7　なんちゃらというメッセージは私の環境でインストール上のゴミなので気にしないで。

(3) 上で定義した関数を呼んでいる。

(4)  時系列xと二次元時系列（水平2次元と時間の3次元）から2次元の相関係数を返す関数。

from numba import autojit
@autojit
の部分がnumba.

(5) ベンチマーク
　　この場合、5.7秒かかっている。
   
　　@autojitの代わりに@jit('f8[:,:](f8[:],f8[:,:,:])')としても速度はかわらなかった。
　　
　　もし上で@autojitを付けない場合、7.34秒だった。
　
　　つまりnumbaを使っていることで、確かに速くなっているが、劇的にではなかった。

　　これは、forでまわしている中身がプリミティブな式ではなくて、numpyの関数corrcoefであるからだろうか。

　　ちなみに、欠損値処理は何もやっていないので、warningが出ている。この場合はプロットするときに陸地のデータはマスクするので、これで構わない。

　　欠損値処理もできるma.corrcoefを使うとかなり遅くなる。
　　numbaを使わない場合で50.2秒、numbaを47.5秒だった。

(6)　実際に相関係数を求める。

(7)  land sea maskを読む関数

(8) land sea maskを読む

(9) プロット 　

　　エルニーニョに典型的な馬蹄形が見られる。インド洋ダイポールへの相関も見られる。


結論
今回試したケースでは、劇的というほどではないが、ちょっとした追加で、確かに速くはなった。場合によっては劇的に速くなる場合もあるだろう。。今後の進展に期待したい技術である。


　　
　　

Lanczos filter と Wakari

Lanczos filterによるlow-pass filterを理解するためにIPythonノートブックを作成した。
https://www.wakari.io/nb/tmiyama/Lanczos_test

結論は、少ない点でもadjustすれば低周波で応答関数が１に近づくようにできるが、基本は点数を大きくする必要があるということだ。

前回はIPython Viewerを使ってノートブックを公開したが、今回はWakariを使用してみた。

Wakariは端末がなくとも、クラウドでいろいろ計算できるので、なかなかおもしろい。

参考文献
Duchon, C. E., 1979: Lanczos Filtering in One and Two Dimensions. Journal of Applied Meteorology, 18, 1016–1022.

気象研究ノート 221号　気象学と海洋物理学で用いられるデータ解析法　伊藤　久徳・見延　庄士郎/著　日本気象学会/刊　A4判　253p

GFD電脳Ruby小物置き場 (Library) Lanczosフィルター
http://davis.gfd-dennou.org/rubygadgets/ja/?(Library)+Lanczos%A5%D5%A5%A3%A5%EB%A5%BF%A1%BC

Torrence and Compo (1998)に出てくる カイ２条分布を理解するためのiPython Notebookを作った。
http://nbviewer.ipython.org/5311742

(追記するかも）



参考文献
Torrence, C., and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bulletin of the American Meteorological Society, 79, 61–78.　http://journals.ametsoc.org/doi/abs/10.1175/1520-0477%25281998%2529079%253C0061%253AAPGTWA%253E2.0.CO%253B2

IPython Notebook

IPython　Notebookで解析して、それをパブリックに表示するデモンストレーション。

作ったnotebookはIPytnon Notebook viewerを使って、ここで見られる。
http://nbviewer.ipython.org/urls/dl.dropbox.com/u/439394/ipyhon/seaice.ipynb

1. データとしては、北極海海氷面積データのcsvでダウンロード
2. pandasでデータを読み、9月の平均データを作る。
3. IPythonのRmagicを使い、線形回帰を行いトレンドラインを引く

参考文献

Python for Data Analysis

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%

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時間間隔のデータがある。
もっとも、データ同化が行われるのは６時間間隔だから、その間は予報データでうめてある。

このデータを見てみた。

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ステップをプロットすると

６時間ごとに値が飛んでいるように見えるのは、６時間ごとに同化（３D-VAR）で解析値に修正していることにともなうものであろう。

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

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

gribmap -0 -i sst1hr.ctl

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

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=０（右端)を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()

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)　３１７ページ

(追記）
私が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"など選んでやれば良い