tag:blogger.com,1999:blog-49183906558188643742024-02-08T13:03:13.090+09:00Ocean Science Hacktmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.comBlogger83125tag:blogger.com,1999:blog-4918390655818864374.post-38330035354719556582014-11-20T10:42:00.000+09:002016-08-20T19:54:59.401+09:00Wavelet解析<br />
<br />
<a href="http://dx.doi.org/10.1175/1520-0477%281998%29079%3C0061%3AAPGTWA%3E2.0.CO%3B2" target="_blank">Torrence and Compo [1998] </a> のwavelet解析手法と、<a href="http://paos.colorado.edu/research/wavelets/" target="_blank">そのツール</a>は我々の分野でよく使われる。私も<a href="http://10.0.3.239/s10236-014-0701-1" target="_blank">Miyama and Miyazawa [2014]</a>をで使用している。その関連した発表をAOGS 2014でした時に、<a href="http://dx.doi.org/10.1175/1520-0477%281998%29079%3C0061%3AAPGTWA%3E2.0.CO%3B2" target="_blank">Torrence and Compo [1998] </a>はスケールが大きいほうを過大評価するバイアスがあり、それに関する論文(<a href="http://dx.doi.org/10.1175/2007JTECHO511.1" target="_blank">Liu et al. 2007</a>と<a href="http://ocgweb.marine.usf.edu/~liu/wavelet.html" target="_blank">関連するweb page</a>)を教えていただいた。<br />
<br />
下の図は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)、各周期が同じくらいの強さであるという合理的な結果が出る。<br />
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEib1tfVs3uOL7OAt-AxfV0Mh7CThPhCpcbvBQ1hv7KPbTR9Ux61at_C1qVW_ldJm1PiTnkGndtfiYri7qbObWppzHjTRglt9z4kZjteWX5h7nijHrLQHP4BZGMLwEmC2RXfO-24Ta2PgAA/s1600/wavelet_test_sine.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em; text-align: center;"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEib1tfVs3uOL7OAt-AxfV0Mh7CThPhCpcbvBQ1hv7KPbTR9Ux61at_C1qVW_ldJm1PiTnkGndtfiYri7qbObWppzHjTRglt9z4kZjteWX5h7nijHrLQHP4BZGMLwEmC2RXfO-24Ta2PgAA/s1600/wavelet_test_sine.png" width="640" /></a><br />
<br />
<br />
<br />
上記の図を作るwavelet解析のツールをpythonに翻訳し、IPython Notebookにしたものはこちら。<br />
<a href="http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/u0wb931fs4qfrpn/wavelet_test_sine.ipynb" target="_blank">http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/u0wb931fs4qfrpn/wavelet_test_sine.ipynb</a><br />
<br />
さらに、NINO SST3のwavelet解析した(Liu et al. 2007のFig 4に対応)物をIPython Notebookしたものはこちら。<br />
<a href="http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/400j051n0sustcy/wavelet_test_ElNino3_Liu.ipynb" target="_blank">http://nbviewer.jupyter.org/urls/dl.dropbox.com/s/400j051n0sustcy/wavelet_test_ElNino3_Liu.ipynb</a><br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh5cAaW7KPU3Yhbp1QvxRmQhqvqdw1mbA1px90Pq3JxldsOJtAmZqWaHP2KOae2KftO2-mx-e1GWJZI_rgT8EWKeKpiDlBbD_SU5g2MXMh6hX0fO_8GI3957897Qn_TANS-G3pyPRwfd18/s1600/nino3_liu.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh5cAaW7KPU3Yhbp1QvxRmQhqvqdw1mbA1px90Pq3JxldsOJtAmZqWaHP2KOae2KftO2-mx-e1GWJZI_rgT8EWKeKpiDlBbD_SU5g2MXMh6hX0fO_8GI3957897Qn_TANS-G3pyPRwfd18/s1600/nino3_liu.png" width="640" /></a></div>
<br />
<br />
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com1tag:blogger.com,1999:blog-4918390655818864374.post-68645838329414739602014-11-11T11:58:00.000+09:002014-11-11T11:58:12.308+09:00Surface velocity data of MOVE/MRI.com at NEAR-GOOS RRTDB<br />
<br />
<a href="http://ds.data.jma.go.jp/gmd/goos/data/database.html" target="_blank">NEAR-GOOSのRegional Real Time Databaseのサイト</a>は、最近更新されて、データのダウンロードがregistration無しでダウンロードできるようになった(自動ダウンロードを行うにはregistrationが必要)。<br />
<br />
ここでは、その中から、気象庁(JMA)の<a href="http://ds.data.jma.go.jp/gmd/goos/data/rrtdb/jma-pro/move_subc_jpn_D.html" target="_blank">Daily Surface Currents in the seas adjacent to Japan</a>をプロットして見る。<br />
<br />
2014年11月6日のデータ "move_subc_jpn_D20141106.txt"をダウンロードして置き、<a href="http://scitools.org.uk/cartopy/" target="_blank">Cartopy</a>でプロットする<a href="http://www.nature.com/news/interactive-notebooks-sharing-the-code-1.16261?WT.mc_id=TWT_NatureNews" target="_blank">IPython Notebook</a>はこちら。<br />
<a href="https://www.dropbox.com/s/bpow09n6x4rm9sy/MOVE_velocity.ipynb?dl=0" target="_blank">https://www.dropbox.com/s/bpow09n6x4rm9sy/MOVE_velocity.ipynb?dl=0</a><br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjOLzWiSXSO2wzk3a-WYezZbwYWM0TY5lzzJfvhOl9X8z2janmDKY4aj9217nbfedGX34-axf2nPadCmwS2dINDgnQYX19XhWXKj3koBjxBLm3bwt3n2myWsyfVGxIopJJdzvAriNFUX1k/s1600/move_vel.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="640" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjOLzWiSXSO2wzk3a-WYezZbwYWM0TY5lzzJfvhOl9X8z2janmDKY4aj9217nbfedGX34-axf2nPadCmwS2dINDgnQYX19XhWXKj3koBjxBLm3bwt3n2myWsyfVGxIopJJdzvAriNFUX1k/s640/move_vel.png" width="640" /></a></div>
<br />
<br />
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-49654527934115673922014-10-30T00:46:00.001+09:002014-10-30T00:46:36.309+09:00Colors (diverging data)<br />
<br />
前回はsequential data用のカラーマップについてテストしたが、今回は正負のアノマリなどをプロットする diverging カラーマップについてテストした。<br />
<br />
ipython notebookはこちら<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/t8um9tdmauu8qpv/divergent_color.ipynb?dl=0" target="_blank">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/t8um9tdmauu8qpv/divergent_color.ipynb?dl=0</a><br />
<br />
テストデータのしては<a href="http://www.ncdc.noaa.gov/sst/" target="_blank">NOAA OISST version 2</a>の海面温度アノマリを使った。データ読み込みとプロットには<a href="http://scitools.org.uk/iris/" target="_blank">Iris</a>を使った。<br />
<br />
Irisはイギリスのソフトのためか、日付変更線をまたぐ場合には少し癖があり、そのデモでもある。<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgVCi6vDvVjcWVbNVcJsShbdimSWi-q6cvq6Om-O03RLijI562LJAFuU6kUAR8mR8ZVWApT8dFxotpMeZjd34-pOCrks05mAmApLTR26VylYowmt58rXJ5-SMU1XwCFbvlNgL7Cxx3zcsQ/s1600/coolwarm.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgVCi6vDvVjcWVbNVcJsShbdimSWi-q6cvq6Om-O03RLijI562LJAFuU6kUAR8mR8ZVWApT8dFxotpMeZjd34-pOCrks05mAmApLTR26VylYowmt58rXJ5-SMU1XwCFbvlNgL7Cxx3zcsQ/s1600/coolwarm.png" height="353" width="640" /></a></div>
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-83830336572401879282014-10-28T15:14:00.001+09:002014-11-01T21:57:14.257+09:00Colors (sequential data)<br />
<br />
最近、カラーマップに関する考察をいくつか読んだ。<br />
<br />
<ul>
<li>THE RAINBOW IS DEAD…LONG LIVE THE RAINBOW!<br /> <a href="http://mycarta.wordpress.com/2012/05/29/the-rainbow-is-dead-long-live-the-rainbow-series-outline/">http://mycarta.wordpress.com/2012/05/29/the-rainbow-is-dead-long-live-the-rainbow-series-outline/</a></li>
<li>Subtleties of Color<br /><a href="http://earthobservatory.nasa.gov/blogs/elegantfigures/2013/08/05/subtleties-of-color-part-1-of-6/">http://earthobservatory.nasa.gov/blogs/elegantfigures/2013/08/05/subtleties-of-color-part-1-of-6/</a></li>
<li>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:<a href="http://dx.doi.org/10.1175/bams-d-13-00155.1" target="_blank">10.1175/bams-d-13-00155.1</a></li>
</ul>
いずれも"rainbow"(またの名を"jet")と呼ばれるカラーマップの弊害を述べている。<br />
<div>
<br /></div>
<div>
それに対する代替はいくつか提案されている。</div>
<div>
<br /></div>
<div>
以下のipython notebookでは、それらを試してみた。<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/5t9jrhr7va3wscj/colors.ipynb?dl=0" target="_blank">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/5t9jrhr7va3wscj/colors.ipynb?dl=0</a></div>
<div>
<br /></div>
<div>
データは<a href="http://apdrc.soest.hawaii.edu/datadoc/ghrsst.php" target="_blank">GHRSST OSTIA SST</a><span style="background-color: white; font-family: 'Helvetica Neue', Helvetica, Arial, sans-serif; font-size: 14px; line-height: 20px; text-align: justify;">という海面温度で、2014年10月15日の</span><span style="font-family: Helvetica Neue, Helvetica, Arial, sans-serif;"><span style="font-size: 14px; line-height: 20px;">117-150E, 17-40Nの範囲を<a href="http://scitools.org.uk/iris/docs/latest/index.html" target="_blank">IRIS</a>を使って<a href="http://apdrc.soest.hawaii.edu/index.php" target="_blank">IPRC APDRC</a>から読み込み、プロットさせている。</span></span></div>
<div>
<span style="font-family: Helvetica Neue, Helvetica, Arial, sans-serif;"><span style="font-size: 14px; line-height: 20px;"><br /></span></span></div>
<div>
<span style="font-family: Helvetica Neue, Helvetica, Arial, sans-serif;"><span style="font-size: 14px; line-height: 20px;">下は"cube1"と呼ばれるカラーマップでプロットした例。</span></span></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhbz9NmRNAZYhyphenhyphensKuPpP1N01UspWD9FxMFANdjdxSlqF6UMY7-bqZasKD_u-Hr-uUpa3TxloO3kvU_sQJOfn9nB803P_c70qmydAxPxtxtL9xRnJkWRB8P3LxqTw3LcRU6TxcSUJUfx8JA/s1600/example.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhbz9NmRNAZYhyphenhyphensKuPpP1N01UspWD9FxMFANdjdxSlqF6UMY7-bqZasKD_u-Hr-uUpa3TxloO3kvU_sQJOfn9nB803P_c70qmydAxPxtxtL9xRnJkWRB8P3LxqTw3LcRU6TxcSUJUfx8JA/s1600/example.png" height="377" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
ちなみに、この日の海面温度には台風19号と台風20号の通過に伴い、海水の低下した場所が見られる。</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjWHnqRY6nqzckVNFWWVE6s4rhEcpvhTju80z9UZgJ19Ne0qV8XM45aIfCY6e_k8qCMQ_fWnazgzy8trKRvRZQA6QKljrP3KkGTZObbttnGQm0o25L4GfJMJD7Flc5MZEGPiUCJ_A5zfzE/s1600/example2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjWHnqRY6nqzckVNFWWVE6s4rhEcpvhTju80z9UZgJ19Ne0qV8XM45aIfCY6e_k8qCMQ_fWnazgzy8trKRvRZQA6QKljrP3KkGTZObbttnGQm0o25L4GfJMJD7Flc5MZEGPiUCJ_A5zfzE/s1600/example2.png" height="377" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<b>Reference</b></div>
<div class="separator" style="clear: both; text-align: left;">
Kuroshio Oyashio Watch 2014/10/20</div>
<div class="separator" style="clear: both; text-align: left;">
<a href="http://www.jamstec.go.jp/jcope/htdocs/e/kow/141020/KOWatch_20141020.html">http://www.jamstec.go.jp/jcope/htdocs/e/kow/141020/KOWatch_20141020.html</a></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div>
<span style="font-family: Helvetica Neue, Helvetica, Arial, sans-serif;"><span style="font-size: 14px; line-height: 20px;"><br /></span></span></div>
<div>
<br /></div>
tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-42239788455292084012014-10-26T08:38:00.000+09:002014-10-26T08:38:15.480+09:00等値線の位置を求める<br />
<br />
<a href="http://dx.doi.org/10.1175/JPO2807.1" target="_blank">Qiu and Chen</a> (2005)などでは、海面高度データのある等値線の位置を黒潮続流の位置とし、Decadal Variationの議論している。<br />
<br />
等値線の位置を決めるのは、渦などの紛らわしい等値線があり、目でみていちいち確認しないで自動的に欲しい等値線の位置を決めるのは難しい。<br />
<br />
ここでは、図作成ソフト<a href="http://matplotlib.org/" target="_blank">matplotlib</a>のcontourルチーンから情報を取り出して、等値線の位置を決める工夫をする。<br />
<br />
参考にしたのはstackoverflowの<br />
<br />
<ul>
<li>How can I get the (x,y) values of the line that is ploted by a contour plot (matplotlib)? <span class="Apple-tab-span" style="white-space: pre;"> </span><a href="http://stackoverflow.com/questions/1560424/how-can-i-get-the-x-y-values-of-the-line-that-is-ploted-by-a-contour-plot-mat">http://stackoverflow.com/questions/1560424/how-can-i-get-the-x-y-values-of-the-line-that-is-ploted-by-a-contour-plot-mat</a> <span class="Apple-tab-span" style="white-space: pre;"> </span> </li>
<li>matplotlib - extracting data from contour lines <span class="Apple-tab-span" style="white-space: pre;"> </span><a href="http://stackoverflow.com/questions/5666056/matplotlib-extracting-data-from-contour-lines">http://stackoverflow.com/questions/5666056/matplotlib-extracting-data-from-contour-lines</a> <span class="Apple-tab-span" style="white-space: pre;"> </span> </li>
<li>Python: find contour lines from matplotlib.pyplot.contour() <span class="Apple-tab-span" style="white-space: pre;"> </span><a href="http://stackoverflow.com/questions/18304722/python-find-contour-lines-from-matplotlib-pyplot-contour/18309914#18309914">http://stackoverflow.com/questions/18304722/python-find-contour-lines-from-matplotlib-pyplot-contour/18309914#18309914</a></li>
</ul>
<span style="white-space: pre;"><span style="white-space: pre;"><br /></span>
ipythonノートブックは</span><br />
<span style="white-space: pre;"><a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/qt804lsufj4whyx/Kuroshio_extension_path.ipynb?dl=0">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/qt804lsufj4whyx/Kuroshio_extension_path.ipynb?dl=0</a></span><br />
<br />
<h3>
Step1 力学高度のデータをferretを使い、読み込み、図示する</h3>
<div>
データはAVISOの絶対力学高度の2009年1月21日のデータを用いる。データは<a href="http://www.aviso.altimetry.fr/en/home.html" target="_blank">AVISO</a>からダウンロードしておく。</div>
<div>
データの読み込みと図示には、<a href="http://oceansciencehack.blogspot.jp/2014/10/pyferret-ipython-ferret-magic.html" target="_blank">前回紹介したferret magic</a>を使う。</div>
<div>
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiLnw2Aid3moNs_y15jIp3ivaFGke0RyClP8U6ooZD4hjmd0LnXofD1qqR7oMs_KZ_UaBzmUFnDMhj5f_A8dohhzSdYdupbfko9shVL5yBhNZa92MkTuoRWoOHRQ2MXavbqsC68dZesJeY/s1600/adt.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiLnw2Aid3moNs_y15jIp3ivaFGke0RyClP8U6ooZD4hjmd0LnXofD1qqR7oMs_KZ_UaBzmUFnDMhj5f_A8dohhzSdYdupbfko9shVL5yBhNZa92MkTuoRWoOHRQ2MXavbqsC68dZesJeY/s1600/adt.png" height="320" width="640" /></a></div>
<div>
<br /></div>
<div>
<br /></div>
このデータでは、0.8mの等値線が黒潮続流位置としてはよさそうである。本当はheat fluxによるsteric heightの変動を考慮にいれないとならないが(Qiu and Chen, 2005)、ここではその手順は省略。<br />
Qiu and Chen (2005)では東経141から153度のデータで議論しているので、その範囲だけ抜き出しておく。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBWt2wOWcifRuZcCatz-VyYRlcviUeeSeykaPzCQdnhrhfsmFm0JcH3TrLJWZvGl1ZIlGbc-bjsohKvILX6lJJm3wLgJsZ35AfH4J0Njgt4MzMFUbWFL-GQinUhB4U-pFbpfrbC10kBTU/s1600/selected_adt.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjBWt2wOWcifRuZcCatz-VyYRlcviUeeSeykaPzCQdnhrhfsmFm0JcH3TrLJWZvGl1ZIlGbc-bjsohKvILX6lJJm3wLgJsZ35AfH4J0Njgt4MzMFUbWFL-GQinUhB4U-pFbpfrbC10kBTU/s1600/selected_adt.png" height="480" width="640" /></a></div>
等値線は、渦や蛇行のなどのため等値線が4つあるが、黒潮続流にあたる本来欲しい線(東西にのびる線)だけを取り出すのが課題である。<br />
<br />
<br />
<h3>
Step 2 Ferretからpythonにデータを読み込む</h3>
<div>
matlabの機能を使うために、ferretからpythonへとデータを読み込む。これもferret magicを用いる。データがちゃんと読めているかを確認もしておく。</div>
<div>
<br /></div>
<div>
<h3>
Step 3 等値線の位置を求める</h3>
</div>
<div>
matlabのcontour routineの機能を使って、等値線の位置を求める。等値線の位置は4つ求まっている。</div>
<div>
<br /></div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgFedLSxi7aonT75NIpx_6bgIGpSA1kJLGd5P5M13CEL0bUJeWbijpviGj3T9Rw9H8A-E5bejZsD8ZOwW7rrOyVF_-AxR9biqCGXp37jheUn4KivHOP-LrVAiU42sxTAQ0YvG1i5B2_kiY/s1600/paths.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgFedLSxi7aonT75NIpx_6bgIGpSA1kJLGd5P5M13CEL0bUJeWbijpviGj3T9Rw9H8A-E5bejZsD8ZOwW7rrOyVF_-AxR9biqCGXp37jheUn4KivHOP-LrVAiU42sxTAQ0YvG1i5B2_kiY/s1600/paths.png" height="425" width="640" /></a></div>
<div>
<h3>
Step 4 等値線の位置情報から、もっとも長い等値線を選ぶ。</h3>
</div>
<div>
等値線の位置が数値的に求まったので、この中から、欲しいpathだけを選ぶ。最も長いpathが黒潮続流にあたる欲しいpathだと考えられる(渦などを避けるために、始点と終点がもっとも離れているpathという条件でも良いかもしれない。)。</div>
<div>
<br /></div>
<div>
まず、pathの長さを求める関数を作成する。</div>
<div>
get_distance は2点間の距離を km で求める関数。</div>
<div>
get_pathlength はget_distanceを用いてpath位置の情報の配列からpathの長さを求める関数。</div>
<div>
<br /></div>
<div>
上の4つのpathの長さを求め、path #2が最も長いことがわかる(約1432km)。</div>
<div>
<blockquote class="tr_bq">
path #= 0 length= 415.785621835 km<br />path #= 1 length= 527.959783762 km<br />path #= 2 length= 1432.08536051 km<br />path #= 3 length= 563.215908518 km</blockquote>
</div>
<div>
これが、求めるpathである。</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjzvqRfB89m7wdEzIAODtJezL4a01y4tL3Mvc4Q9jjld91GLZow75uSZba0cL7tfEH8pSJYwrC8RSHEeTxjawzMeET9Jh22mvSSvxXQBmA4oLEEf2msJRX8imvUd9yBXWaw45cjz-j1E70/s1600/wantedpath.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjzvqRfB89m7wdEzIAODtJezL4a01y4tL3Mvc4Q9jjld91GLZow75uSZba0cL7tfEH8pSJYwrC8RSHEeTxjawzMeET9Jh22mvSSvxXQBmA4oLEEf2msJRX8imvUd9yBXWaw45cjz-j1E70/s1600/wantedpath.png" height="426" width="640" /></a></div>
<div>
<br /></div>
<div>
<br /></div>
<div>
<b>Reference:</b></div>
<div>
<div>
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:<a href="http://dx.doi.org/10.1175/JPO2807.1" target="_blank">10.1175/JPO2807.1</a></div>
<div style="font-weight: bold;">
<br /></div>
</div>
<div>
<br /></div>
<br />
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-33061608010861912132014-10-18T18:02:00.000+09:002014-10-18T18:02:12.705+09:00pyferret, ipython ferret magic<br />
<br />
<a href="http://ferret.pmel.noaa.gov/Ferret/documentation/pyferret" target="_blank">pyferret</a>は、大気海洋分野のお絵描きソフト<a href="http://ferret.pmel.noaa.gov/Ferret/home" target="_blank">ferret</a>の、pythonとの結合を深めたバージョンである。<br />
<br />
さらに<a href="https://github.com/PBrockmann/ipython-ferretmagic" target="_blank">ipython-ferretmagic</a>は<a href="http://ipython.org/" target="_blank">ipython</a>の中でferretを使うのを可能にしてくれる。<br />
<br />
サンプルのノートブックはこちら。<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/avefneoyxdxf7ov/ferret_magic_example.ipynb?dl=0">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/avefneoyxdxf7ov/ferret_magic_example.ipynb?dl=0</a><br />
<br />
この中で、Levitusのデータを読み込んで、9月の温度26度の等温線を求めてプロットしている。<br />
ferretからpythonにデータを読み込むことも可能である。<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdggoUBjX8bdiV8tVdxXMzbeDdbVfd5uiYPqCnWGLm0h79WL_3jasYKkNr3mw_boN8bGkENGNKR3HoR99BnPgiLzoFUJaOueWpXAWSzNWv7wlb4mFeYVvrvObQXDtiT1ZLhavI5OGyB4E/s1600/ithothermdepth26.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdggoUBjX8bdiV8tVdxXMzbeDdbVfd5uiYPqCnWGLm0h79WL_3jasYKkNr3mw_boN8bGkENGNKR3HoR99BnPgiLzoFUJaOueWpXAWSzNWv7wlb4mFeYVvrvObQXDtiT1ZLhavI5OGyB4E/s1600/ithothermdepth26.png" height="240" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-2430324988469119672014-09-06T17:55:00.000+09:002014-09-06T17:55:53.481+09:00surface drifter<br />
<br />
<br />
NOAA, Atlantic Oceanographic and Meteorological Laboratory (AOML)のGlobal Drifter Program (GDP), <a href="http://www.aoml.noaa.gov/phod/dac/dacdata.php" target="_blank">Drifter Data Assembly Center (DAC)</a>には表層漂流ブイのアーカイブがある。このデータを読み込んで図示してみる。<br />
<br />
<a href="http://www.aoml.noaa.gov/envids/gld/FtpInterpolatedInstructions.php" target="_blank">サイト</a>からftpでデータをあらかじめダウンロードしておく。このうち、buoydata_10001_mar14.dat (ブイ1万個目以降から2014年3月までのデータ)を例として扱う。<br />
<br />
IPython Notebookはこちら。<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/uzv1bktgazcf821/Test_buoys.ipynb?dl=0" target="_blank">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/uzv1bktgazcf821/Test_buoys.ipynb?dl=0</a><br />
<br />
まず、<a href="http://pandas.pydata.org/" target="_blank">pandas</a>でデータを読み込む。データが多量なので、読み込むまで時間がかかるかもしれない。<br />
<br />
データ(data)の内で、一部に緯度が北緯25から45度、東経120から160度が含まれるデータを選びだす(selecteddata)。そのうちで、ユニークなid番号を持つもののは851個ある。<br />
<br />
あらためて、もとのデータ(data)から、そのid番号を持つデータを抜き出し(wanteddata)、pickle形式で保存する。<br />
<br />
次に、そのデータの中から、例として最初の10個だけ<a href="http://scitools.org.uk/cartopy/" target="_blank">cartopy</a>を用いてプロットした。<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmVSFJ0xZUBHWtq-oAfjFJDBDtF9Jdt0xmf_Yii57OjZVClTPAGd5gNr5YpiDKQ1FCu02nAj3rag7FuJh7IbcBs1z14FDyxXHvs21UICB3yfMeCqos5kh4B4n_r_logX4f7RP6vIAfPiQ/s1600/buoy.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgmVSFJ0xZUBHWtq-oAfjFJDBDtF9Jdt0xmf_Yii57OjZVClTPAGd5gNr5YpiDKQ1FCu02nAj3rag7FuJh7IbcBs1z14FDyxXHvs21UICB3yfMeCqos5kh4B4n_r_logX4f7RP6vIAfPiQ/s1600/buoy.png" height="241" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
参考文献</div>
<div class="separator" style="clear: both; text-align: left;">
<a href="http://www.amazon.co.jp/gp/product/4873116554/ref=as_li_ss_tl?ie=UTF8&camp=247&creative=7399&creativeASIN=4873116554&linkCode=as2&tag=climatechange-22">Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理</a><img alt="" border="0" src="http://ir-jp.amazon-adsystem.com/e/ir?t=climatechange-22&l=as2&o=9&a=4873116554" height="1" style="border: none !important; margin: 0px !important;" width="1" /></div>
<a href="http://www.amazon.co.jp/gp/product/4873116554/ref=as_li_ss_il?ie=UTF8&camp=247&creative=7399&creativeASIN=4873116554&linkCode=as2&tag=climatechange-22"><img border="0" src="http://ws-fe.amazon-adsystem.com/widgets/q?_encoding=UTF8&ASIN=4873116554&Format=_SL110_&ID=AsinImage&MarketPlace=JP&ServiceVersion=20070822&WS=1&tag=climatechange-22" /></a><img alt="" border="0" src="http://ir-jp.amazon-adsystem.com/e/ir?t=climatechange-22&l=as2&o=9&a=4873116554" height="1" style="border: none !important; margin: 0px !important;" width="1" /><br />
<br />
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-5552007694656497142014-08-07T23:36:00.000+09:002014-08-15T20:37:09.105+09:00IrisとETOPO1<br />
<br />
<a href="http://scitools.org.uk/iris/" target="_blank">Iris</a>はイギリスで開発されているpythonベースの気象・海洋データ解析パッケージである。<br />
<br />
Irisを使って地形データETOPO1を読み込んで、プロットしたり、補間してみた。<br />
<br />
ETOPO1は1/60°の海底地形で、<a href="http://www.ngdc.noaa.gov/mgg/global/global.html" target="_blank">NOAA NGDC</a>のサイトからETOPO1 Bedrock(氷無しの地形)、cell-registered(グリッドセルの位置で定義される)の中から、netcdfファイル ETOPO1_Bed_c_gmt4.grdをダウンロードした。<br />
<br />
IPython Notebookは<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/67khdbgzc2bf0sp/Read_ETOPO1.ipynb" target="_blank">こちら</a>である。<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/67khdbgzc2bf0sp/Read_ETOPO1.ipynb">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/67khdbgzc2bf0sp/Read_ETOPO1.ipynb</a><br />
<br />
下の図はプロットの一例。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjLCVh0KJ1JtsujlXG1K6gs-yTaKM2nsNdc6nWdipaS_WH3CDyO3sHmzdY1mAg80U4FQlXKZY2WYDvNDz7uOMVB5LuCKTcfPfkr9XEp2RdVG7FzRq1LQCmrefAEoAPL790KiYYG2u3bbSk/s1600/ETOPO1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjLCVh0KJ1JtsujlXG1K6gs-yTaKM2nsNdc6nWdipaS_WH3CDyO3sHmzdY1mAg80U4FQlXKZY2WYDvNDz7uOMVB5LuCKTcfPfkr9XEp2RdVG7FzRq1LQCmrefAEoAPL790KiYYG2u3bbSk/s1600/ETOPO1.png" height="227" width="320" /></a></div>
<span class="Apple-tab-span" style="white-space: pre;"> </span>tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-66437792234165150122014-06-14T08:14:00.000+09:002014-06-14T08:26:07.335+09:00PDO index<br />
<br />
Pacific Decadal Oscillation (PDO) index を<a href="http://jisao.washington.edu/pdo/" target="_blank">JISAOのPDOページ</a>から<a href="http://jisao.washington.edu/pdo/PDO.latest" target="_blank">データ</a>を取得し、ある月から最新のデータまで、時系列をプロットする。<br />
<br />
ipython notobookはこちら<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/y1lol5bl7h4xqdi/PDO.ipynb" target="_blank">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/y1lol5bl7h4xqdi/PDO.ipynb</a><br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnqFltwvGPIsCfV7nJkt6sslNotmc1oOuGzirORUOXoR-77i7Ec9zsN3sDz9q1kw2-b66bgd4vkVClbZKSMg-Z5Pd7BSKfucnskRNgzUGZb6X31P5PzOJXi2bfa9J0QeYfWv_yOJsSsoA/s1600/PDOsince2000-1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgnqFltwvGPIsCfV7nJkt6sslNotmc1oOuGzirORUOXoR-77i7Ec9zsN3sDz9q1kw2-b66bgd4vkVClbZKSMg-Z5Pd7BSKfucnskRNgzUGZb6X31P5PzOJXi2bfa9J0QeYfWv_yOJsSsoA/s1600/PDOsince2000-1.png" height="213" width="320" /></a></div>
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-91742600506962058922014-02-11T18:05:00.000+09:002014-02-11T18:05:57.732+09:00bokehを使って、串本浦神の潮位差をプロット<br />
<br />
串本と浦神の潮位差のデータを<a href="http://bokeh.pydata.org/" target="_blank">bokeh</a>を使ってプロットしてみる。bokehはpythonのプロットライブラリの一つで、ズームなどが可能なインタラクティブなプロットをweb上につくることができる。(bokehは日本語のボケから来てるらしい。何故そんな名前にしたんだか。)<br />
<br />
ipython notebookはこちら。<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/s/rry6nt3xs49o6s6/KushimotoUragami.ipynb" target="_blank">http://nbviewer.ipython.org/urls/dl.dropbox.com/s/rry6nt3xs49o6s6/KushimotoUragami.ipynb</a><br />
<br />
ズームなどができるのを確認してみてください。<br />
<br />
データは<a href="http://oceansciencehack.blogspot.jp/2012/07/pandas-version-0.html" target="_blank">以前のように</a>pandasで読み込んでいる。<br />
daily dataの黒線で、30日の移動平均を赤線でプロットしている。ここでの30日平均は基準日を中心とした30日平均ではなく、基準日をさかのぼっての30日平均であることに注意。<br />
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-6914804869621223432013-04-17T00:31:00.000+09:002013-04-17T00:31:44.599+09:00Numba<br />
<br />
<br />
Pythonでforをまわすような数値計算を行うとmatlabなどと同じで非常に遅い。<br />
それを克服するために<a href="http://pypy.org/" target="_blank">pypy</a>や<a href="http://www.cython.org/" target="_blank">cython</a>などがある。しかし基本的にコードを書き換える必要があるので、もっと楽をしたい。<br />
<br />
最近<a href="http://numba.pydata.org/" target="_blank">numba</a>というものを知った。ここのサンプルにある例だと、<br />
<br />
from numba import autojit<br />
<br />
@autojit<br />
<br />
のおまじないだけで関数が劇的に早くなる。<br />
<br />
少し、自分で例で実験してみた。<br />
全球のSSTに対して、<a href="http://www.data.jma.go.jp/gmd/cpd/db/elnino/index/dattab.html" target="_blank">NINO3</a>の相関をとるということをやってみる(これはnumbaの実験のために作った例であり、相関をとるだけならもっとエレガントな方法があるだろう)。<br />
<br />
そのノートブックが<br />
<a href="http://nbviewer.ipython.org/5396344">http://nbviewer.ipython.org/5396344</a><br />
<br />
説明<br />
(1) 図をプロットするためのおまじない。<br />
<br />
(2) SSTアノマリとNINO3を読む関数を定義。<br />
元データは <a href="http://www.esrl.noaa.gov/psd/data/gridded/data.noaa.oisst.v2.html" target="_blank">NOAA OISST v2</a>でmonthly dataのsst.mnmean.ncを使う。<br />
データは<a href="http://opengrads.org/wiki/index.php?title=Python_Interface_to_GrADS" target="_blank">Python Interface to GrADS</a>を使って読む。<br />
1982から2010年のデータから気候値からのanomaly(Bin Guan's GrADS Script Libraryの<a href="http://www.atmos.umd.edu/~bguan/grads/GrADS_Scripts.html#deseason.gs" target="_blank">deseason</a>を使った。古いバージョンを使ったので、今は使い方が違うかも知れない)とNINO3 SSTアノマリを定義し、返す。<br />
ここらへんは<a href="http://oceansciencehack.blogspot.jp/2010/02/blog-post.html" target="_blank">以前の記事</a>も参照。<br />
<span style="font-size: 13px; line-height: 20px; white-space: pre-wrap;">/usr/local/lib/python2.7 なんちゃらというメッセージは私の環境でインストール上のゴミなので気にしないで。</span><br />
<br />
(3) 上で定義した関数を呼んでいる。<br />
<br />
(4) 時系列xと二次元時系列(水平2次元と時間の3次元)から2次元の相関係数を返す関数。<br />
<br />
from numba import autojit<br />
@autojit<br />
の部分がnumba.<br />
<br />
<b>(5) ベンチマーク</b><br />
この場合、5.7秒かかっている。<br />
<br />
<complete id="goog_537716481"> @autojitの代わりに@jit('f8[:,:](f8[:],f8[:,:,:])')としても速度はかわらなかった。</complete><br />
<br />
もし上で@autojitを付けない場合、7.34秒だった。<br />
<br />
つまりnumbaを使っていることで、確かに速くなっているが、劇的にではなかった。<br />
<br />
これは、forでまわしている中身がプリミティブな式ではなくて、numpyの関数corrcoefであるからだろうか。<br />
<br />
ちなみに、欠損値処理は何もやっていないので、warningが出ている。この場合はプロットするときに陸地のデータはマスクするので、これで構わない。<br />
<br />
欠損値処理もできるma.corrcoefを使うとかなり遅くなる。<br />
numbaを使わない場合で50.2秒、numbaを47.5秒だった。<br />
<br />
(6) 実際に相関係数を求める。<br />
<br />
(7) land sea maskを読む関数<br />
<br />
(8) land sea maskを読む<br />
<br />
(9) プロット <br />
<br />
エルニーニョに典型的な馬蹄形が見られる。インド洋ダイポールへの相関も見られる。<br />
<br />
<br />
<b>結論</b><br />
今回試したケースでは、劇的というほどではないが、ちょっとした追加で、確かに速くはなった。場合によっては劇的に速くなる場合もあるだろう。。今後の進展に期待したい技術である。<br />
<br />
<br />
<br />
tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-39701106419599209312013-04-12T22:40:00.000+09:002013-04-12T22:50:21.370+09:00Lanczos filter と Wakari<br />
<br />
<br />
Lanczos filterによるlow-pass filterを理解するためにIPythonノートブックを作成した。<br />
<a href="https://www.wakari.io/nb/tmiyama/Lanczos_test">https://www.wakari.io/nb/tmiyama/Lanczos_test</a><br />
<br />
結論は、少ない点でもadjustすれば低周波で応答関数が1に近づくようにできるが、基本は点数を大きくする必要があるということだ。<br />
<br />
前回はIPython Viewerを使ってノートブックを公開したが、今回は<a href="http://continuum.io/blog/ipython-notebook-sharing-in-wakari" target="_blank">Wakari</a>を使用してみた。<br />
<br />
Wakariは端末がなくとも、クラウドでいろいろ計算できるので、なかなかおもしろい。<br />
<br />
<b>参考文献</b><br />
Duchon, C. E., 1979: Lanczos Filtering in One and Two Dimensions. Journal of Applied Meteorology, 18, 1016–1022.<br />
<div>
<br /></div>
<div>
気象研究ノート 221号 <a href="http://www.teda1.com/tsumura-shoten/shop/index.cgi?mode=item_view&no=BK-22217" target="_blank">気象学と海洋物理学で用いられるデータ解析法</a> 伊藤 久徳・見延 庄士郎/著 日本気象学会/刊 A4判 253p<br />
<br />
GFD電脳Ruby小物置き場 (Library) Lanczosフィルター<br />
<a href="http://davis.gfd-dennou.org/rubygadgets/ja/?(Library)+Lanczos%A5%D5%A5%A3%A5%EB%A5%BF%A1%BC">http://davis.gfd-dennou.org/rubygadgets/ja/?(Library)+Lanczos%A5%D5%A5%A3%A5%EB%A5%BF%A1%BC</a></div>
tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-82474408414783636042013-04-05T01:19:00.002+09:002013-04-05T01:20:24.634+09:00<br />
Torrence and Compo (1998)に出てくる カイ2条分布を理解するためのiPython Notebookを作った。<br />
<a href="http://nbviewer.ipython.org/5311742">http://nbviewer.ipython.org/5311742</a><br />
<br />
(追記するかも)<br />
<br />
<br />
<br />
<b>参考文献</b><br />
Torrence, C., and G. P. Compo, 1998: A Practical Guide to Wavelet Analysis. Bulletin of the American Meteorological Society, 79, 61–78. <a href="http://journals.ametsoc.org/doi/abs/10.1175/1520-0477%25281998%2529079%253C0061%253AAPGTWA%253E2.0.CO%253B2">http://journals.ametsoc.org/doi/abs/10.1175/1520-0477%25281998%2529079%253C0061%253AAPGTWA%253E2.0.CO%253B2</a><br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-27133768084468457952013-04-02T23:16:00.000+09:002013-04-02T23:50:25.202+09:00IPython Notebook<br />
<br />
<br />
<a href="http://ipython.org/notebook.html" target="_blank">IPython Notebook</a>で解析して、それをパブリックに表示するデモンストレーション。<br />
<br />
作ったnotebookは<a href="http://nbviewer.ipython.org/" target="_blank">IPytnon Notebook viewer</a>を使って、ここで見られる。<br />
<a href="http://nbviewer.ipython.org/urls/dl.dropbox.com/u/439394/ipyhon/seaice.ipynb">http://nbviewer.ipython.org/urls/dl.dropbox.com/u/439394/ipyhon/seaice.ipynb</a><br />
<br />
<ol>
<li>データとしては、<a href="http://www.ijis.iarc.uaf.edu/jp/seaice/extent.htm" target="_blank">北極海海氷面積データ</a>のcsvでダウンロード</li>
<li><a href="http://pandas.pydata.org/" target="_blank">pandas</a>でデータを読み、9月の平均データを作る。</li>
<li>IPythonの<a href="http://www.randalolson.com/2013/01/14/filling-in-pythons-gaps-in-statistics-packages-with-rmagic/?utm_source=rss&utm_medium=rss&utm_campaign=filling-in-pythons-gaps-in-statistics-packages-with-rmagic" target="_blank">Rmagic</a>を使い、線形回帰を行いトレンドラインを引く</li>
</ol>
<div>
<br /></div>
<br />
<br />
<div>
<b>参考文献</b></div>
<div>
<a href="http://www.amazon.co.jp/gp/product/1449319793/ref=as_li_ss_tl?ie=UTF8&camp=247&creative=7399&creativeASIN=1449319793&linkCode=as2&tag=climatechange-22">Python for Data Analysis</a><img alt="" border="0" height="1" src="http://www.assoc-amazon.jp/e/ir?t=climatechange-22&l=as2&o=9&a=1449319793" style="border: none !important; margin: 0px !important;" width="1" /></div>
<a href="http://www.amazon.co.jp/gp/product/1449319793/ref=as_li_ss_il?ie=UTF8&camp=247&creative=7399&creativeASIN=1449319793&linkCode=as2&tag=climatechange-22"><img border="0" src="http://ws.assoc-amazon.jp/widgets/q?_encoding=UTF8&ASIN=1449319793&Format=_SL110_&ID=AsinImage&MarketPlace=JP&ServiceVersion=20070822&WS=1&tag=climatechange-22" /></a><img alt="" border="0" height="1" src="http://www.assoc-amazon.jp/e/ir?t=climatechange-22&l=as2&o=9&a=1449319793" style="border: none !important; margin: 0px !important;" width="1" />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-60732479442465981552012-09-28T22:26:00.003+09:002012-09-28T22:27:21.664+09:00Sageで方程式<a href="http://www.sagemath.org/" target="_blank">Sage</a>で200万が5年後に202万になる利息を計算する。<br />
<br />
<blockquote class="tr_bq">
x=var(x)<br />
q=solve(200*(1+x)**5==202,x)<br />
for n in range(len(q)):<br />
print N(q[n].rhs())</blockquote>
<br />
<br />
結果<br />
<br />
-0.690367429042489 + 0.952951066209182*I<br />
-1.81062859479078 + 0.588956148532726*I<br />
-1.81062859479078 - 0.588956148532726*I<br />
-0.690367429042489 - 0.952951066209181*I<br />
0.00199204766653335<br />
<br />
したがって0.00199204766653335=約0.199%<br />
<br />
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-62676017161030505072012-07-26T01:36:00.000+09:002012-07-26T01:37:50.326+09:00NCEP Climate Forecast System Reanalysis<br />
<br />
<br />
NCEP Climate Forecast System Reanalysisというデータがある。<br />
<a href="http://cfs.ncep.noaa.gov/cfsr/">http://cfs.ncep.noaa.gov/cfsr/</a><br />
<br />
大気海洋モデル結合モデルに対してデータ同化が行われているのが特色だ。<br />
分解能は従来のNCEPの再解析よりも高い。<br />
<br />
大気: 水平分解能<span class="Apple-style-span" style="font-family: helvetica, arial, verdana, sans-serif;">~38 km (T382)、 鉛直64 levels</span><br />
<span class="Apple-style-span" style="font-family: helvetica, arial, verdana, sans-serif;">海洋: 熱帯で0.25度、それ以外で0.5度の水平分解能、鉛直</span><span class="Apple-style-span" style="font-family: helvetica, arial, verdana, sans-serif;">40 levels</span><br />
<br />
また変数によっては1時間間隔のデータがある。<br />
もっとも、データ同化が行われるのは6時間間隔だから、その間は予報データでうめてある。<br />
<br />
このデータを見てみた。<br />
<br />
UCARのサイトの<br />
<a href="http://rda.ucar.edu/pub/cfsr.html">http://rda.ucar.edu/pub/cfsr.html</a><br />
NCEP Climate Forecast System Reanalysis (CFSR) Selected Hourly Time-Series Products, January 1979 to December 2010<a href="http://rda.ucar.edu/datasets/ds093.1/">http://rda.ucar.edu/datasets/ds093.1/</a><br />
から、2001年の12月のSSTデータ(中身は海面下5mの温度)<br />
ocnsst.gdas.200112.grb2<br />
をダウンロードした。<br />
<br />
これは<a href="http://www.nco.ncep.noaa.gov/pmb/docs/grib2/" target="_blank">grib2形式</a>なので、これを<a href="http://blog.livedoor.jp/rootan2007/archives/51161844.html" target="_blank">wgrib2</a>でgrads形式に変換する。<br />
<br />
<span class="Apple-style-span" style="background-color: white;">g2ctl ocnsst.gdas.200112.grb2 > sst1hr.ctl</span><br />
<span class="Apple-style-span" style="background-color: white;">gribmap -i sst1hr.ctl</span><br />
<br />
で、grads ctlファイルとindexファイル<br />
<br />
sst1hr.ctl<br />
ocnsst.gdas.200112.grb2.idx<br />
<br />
ができる。水平分解能は0.5x0.5になっていた。<br />
<br />
gradsで最初のステップをプロットすると、こんなかんじ<br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg3-Ar6t_d_fHTlzSpfoasNBKBhV7gFF8dAptRKnDOXaX8QZPqrUutfwY1RDfT6jW9Jqngfxl20Ar3DttkzACaJ3IYfrnlcUJ3ke99Vdjtl3d_vzlNnua9TRDPGHpVFPP8AbRFxtQJPsww/s1600/sst1.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="308" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg3-Ar6t_d_fHTlzSpfoasNBKBhV7gFF8dAptRKnDOXaX8QZPqrUutfwY1RDfT6jW9Jqngfxl20Ar3DttkzACaJ3IYfrnlcUJ3ke99Vdjtl3d_vzlNnua9TRDPGHpVFPP8AbRFxtQJPsww/s400/sst1.png" width="400" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
0N, 180Eの最初の48ステップをプロットすると</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjenLHk3BDqpOnIfdfz3G12TsR9rxCdgRhYB-p2XcKcHhzM06szT36XjX2eeB6ZPxeoQ_U5iHFxn0LqfEZxJasmGdpHGGuUeb7iCe4XSgYR7f7l7sJFcyLZOGmvVAxRwUl97bH-JM8CJLs/s1600/sst2.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="247" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjenLHk3BDqpOnIfdfz3G12TsR9rxCdgRhYB-p2XcKcHhzM06szT36XjX2eeB6ZPxeoQ_U5iHFxn0LqfEZxJasmGdpHGGuUeb7iCe4XSgYR7f7l7sJFcyLZOGmvVAxRwUl97bH-JM8CJLs/s320/sst2.png" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
6時間ごとに値が飛んでいるように見えるのは、6時間ごとに同化(3D-VAR)で解析値に修正していることにともなうものであろう。</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
grads形式に変換する時に、解析値だけ取り出すには"-0"(ハイフンゼロ)をつけて</div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
</div>
<div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">
<span class="Apple-style-span" style="background-color: white;">g2ctl -0 ocnsst.gdas.200112.grb2 > sst1hr.ctl</span></div>
<div style="margin-bottom: 0px; margin-left: 0px; margin-right: 0px; margin-top: 0px;">
<span class="Apple-style-span" style="background-color: white;">gribmap -0 -i sst1hr.ctl</span></div>
<br />
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
この場合、0N, 180Eの最初の8ステップ(6時間ごとなので、48時間分)をプロットすると</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhbW3EMWaGRCThBU0E1LPLFJWol3LKn3VB0gtFpPpjKyljNiU6TizJkdssOg3A6KwAslFSR_9tN7AIRNYO9DGT1UFsu-36geqxiOfOsd3lFM6QBFE2_DldqLdrKksHVVsWYnEGJD0qqFIw/s1600/sst3.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="247" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhbW3EMWaGRCThBU0E1LPLFJWol3LKn3VB0gtFpPpjKyljNiU6TizJkdssOg3A6KwAslFSR_9tN7AIRNYO9DGT1UFsu-36geqxiOfOsd3lFM6QBFE2_DldqLdrKksHVVsWYnEGJD0qqFIw/s320/sst3.png" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-2896753967377970432012-07-24T23:37:00.001+09:002012-07-24T23:37:46.165+09:00Wvisで気象情報を3次元可視化<br />
<br />
<br />
Wvisという、気象庁の<a href="http://www.jma.go.jp/jma/kishou/know/whitep/1-3-6.html" target="_blank">メソモデル(MSM)</a>を3次元可視化してくれるソフトがある。<br />
今のところMSMしか可視化できないが、MSMなら手軽に3次元可視化できる。<br />
<br />
Wvisは以下からダウンロードできる。<br />
<a href="http://www.cybernet.co.jp/avs/download/wvis.html">http://www.cybernet.co.jp/avs/download/wvis.html</a>
<br />
<br />
Windowsのみのソフトだが、私の試した限り、Linux上でも、<a href="http://ja.wikipedia.org/wiki/Wine" target="_blank">wine</a>を用いることで問題なく使えた<br />
<br />
Wvisについては、気象学会の「天気」にも解説がある。<br />
<a href="http://www.metsoc.jp/tenki/pdf/2011/2011_09_0073.pdf" style="background-color: white;">http://www.metsoc.jp/tenki/pdf/2011/2011_09_0073.pdf</a><br />
<br />
使い方については、Wvisのマニュアルまたは、「天気」の記事を参照されたい。<br />
<br />
MSMのデータは<br />
<a href="http://database.rish.kyoto-u.ac.jp/arch/jmadata/gpv-original.html">http://database.rish.kyoto-u.ac.jp/arch/jmadata/gpv-original.html</a>
<br />
または<br />
<a href="http://gpvjma.ccs.hpcc.jp/~gpvjma/">http://gpvjma.ccs.hpcc.jp/~gpvjma/</a><br />
でgrib2形式のデータが入手可能である。<br />
<br />
例として、下の動画は2012年7月11-14日にかけて九州に豪雨をもたらした暖かく湿った空気(相当温位355K)をプロットしたものである。説明についてはYoutubeのキャプションの欄を参照されたい。<a href="http://www.youtube.com/watch?v=GO0ZrA3wx5s">http://www.youtube.com/watch?v=GO0ZrA3wx5s</a><br />
<br />
<div class="separator" style="clear: both; text-align: center;">
<iframe allowfullscreen='allowfullscreen' webkitallowfullscreen='webkitallowfullscreen' mozallowfullscreen='mozallowfullscreen' width='320' height='266' src='https://www.youtube.com/embed/GO0ZrA3wx5s?feature=player_embedded' frameborder='0'></iframe></div>
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-33327047036952454172012-07-24T01:03:00.000+09:002012-07-24T01:06:05.394+09:00浦神と串本の日平均潮位差をpandasで<br />
<br />
<a href="http://oceansciencehack.blogspot.jp/2012/03/pandas.html" target="_blank">以前書いた</a>、<a href="http://pandas.pydata.org/" target="_blank">pandas</a>での時系列処理の続き。バージョンアップ(これを書いてる時点でversion 0.8.1)により、さらに便利になった。<br />
<br />
<a href="http://www.data.kishou.go.jp/db/kobe/kuroshio/chouisa/chouisa.php" target="_blank">浦神と串本の日平均潮位差</a>をプロットしてみる。前回のようにwebから直接読むこともできるが、今回は<a href="http://www.data.kishou.go.jp/db/kobe/kuroshio/chouisa/chouisa.dat">http://www.data.kishou.go.jp/db/kobe/kuroshio/chouisa/chouisa.dat</a>からchouisa.datをダウンロードしてそれを読むことにする。<br />
<div>
<br /></div>
<div>
下がプログラム。</div>
<div>
read_psvでデータを読んでいる。<br />
前回は日付とデータの間にあるスペースをスマートに処理できなかったが、sep='\s*' (任意の数のスペース)を使えばよかったらしい。<br />
column=0(右端)をindexとして用い(index_col=0)、日付データとして解釈する(parse_dates=Trule)。<br />
”-”を欠損値と解釈する。<br />
<br />
version 0.8以降、月平均を計算するのが簡単になり、<span class="Apple-style-span" style="font-family: monospace; white-space: pre;">data.resample("M",how='mean')で良い。</span></div>
<div>
<br /></div>
<div>
<div>
<br />
<pre class="prettyprint linenums">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()
</pre>
</div>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjxKDQDMbIH18NdCeiuDfYv-3Z1gg3HOSd8tEI5HDUIfIf1QWTUwhmKQykBESRJIUQejTvum3ROEeFa6yxM38PUmq_XymLImks8PxU2g0JMWCI3u8egA__1sv6UgQaZHsLRxy8SSGxwH8I/s1600/choisa.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjxKDQDMbIH18NdCeiuDfYv-3Z1gg3HOSd8tEI5HDUIfIf1QWTUwhmKQykBESRJIUQejTvum3ROEeFa6yxM38PUmq_XymLImks8PxU2g0JMWCI3u8egA__1sv6UgQaZHsLRxy8SSGxwH8I/s320/choisa.png" width="320" /></a></div>
<br /></div>tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-21190595395254723572012-04-28T16:57:00.003+09:002012-04-28T16:58:05.410+09:00pythonでssh<br />
<br />
pythonで外部のサーバーにsshでアクセスして、そこでコマンドを実行する。<br />
そのために<a href="http://www.lag.net/paramiko/" target="_blank">paramiko</a>というパッケージを用いる。<br />
<br />
以下はその例。<br />
外部サーバにアクセス(ホスト名***, user名xxx, password +++)にアクセスし、コマンド"du . | sort -n -r"(どのディレクトリが容量を使っているかを調べる)を実行して、標準出力に表示する。<br />
<br />
<pre class="prettyprint linenums">#!/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()
</pre>
<br />
参考文献<br />
<a href="http://www.amazon.co.jp/gp/product/4774138053/ref=as_li_ss_tl?ie=UTF8&tag=climatechange-22&linkCode=as2&camp=247&creative=7399&creativeASIN=4774138053">Python ポケットリファレンス (Pocket Reference)</a><img alt="" border="0" height="1" src="http://www.assoc-amazon.jp/e/ir?t=climatechange-22&l=as2&o=9&a=4774138053" style="border: none !important; margin: 0px !important;" width="1" /> 317ページ<br />
<br />
(追記)<br />
私がparamikoを使っている時、時々timeoutの問題が起こることがあった。<br />
その場合<br />
<a href="http://www.stillhq.com/commentform.cgi?post=python/paramiko/000004">http://www.stillhq.com/commentform.cgi?post=python/paramiko/000004</a><br />
などを参考にして、回避する。<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-67486102393945495482012-04-28T11:41:00.000+09:002012-04-28T11:43:18.410+09:00ENDNOTE: MS Wordの文献リストのhanging indent<br />
<br />
MS WORDでENDNOTEを使ってHanging indentをするやり方がずっとわからなかったのだが、やっとわかった。<br />
<br />
Hanging indentというのはこういうの<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiE3UY-_Baop-uvsFfOgkyXUUmp2NpC-3eFn2oc3NgxBK_8iWq9ZUNvXBNua4y8Wa_vFL_bbt6WOm8rDq8ae_WQzI9eXLLMnkChW1mkAPzkP7rrLy_AVVSQiWKRWUJ0BFQRCTSH_RY9y5Y/s1600/c1.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="137" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiE3UY-_Baop-uvsFfOgkyXUUmp2NpC-3eFn2oc3NgxBK_8iWq9ZUNvXBNua4y8Wa_vFL_bbt6WOm8rDq8ae_WQzI9eXLLMnkChW1mkAPzkP7rrLy_AVVSQiWKRWUJ0BFQRCTSH_RY9y5Y/s640/c1.PNG" width="640" /></a></div>
<br />
<br />
設定はWordのEndnoteタブのBibliographyのLayoutでHanging indentの大きさを指定するとともに<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhU9Gt84pA7hpKZ_C6ieL6V67sasiD0NlbWsj59eES_z8zPKmv_bbHe4p1SmXibkX3085hIFd_h936bd54BcuLv7g59bOtb8PloyWHjvgdU36EUcOOwAk0x51HBxECCT1oJzR9XhkOBj2M/s1600/c2.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="241" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhU9Gt84pA7hpKZ_C6ieL6V67sasiD0NlbWsj59eES_z8zPKmv_bbHe4p1SmXibkX3085hIFd_h936bd54BcuLv7g59bOtb8PloyWHjvgdU36EUcOOwAk0x51HBxECCT1oJzR9XhkOBj2M/s320/c2.PNG" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<div class="separator" style="clear: both; text-align: left;">
Endnote本体のStyle ManagerのEditorで、BibliographyのLayoutの右隅(これに気が付かなかった)のHanging indentで"All paragraphs"など選んでやれば良い</div>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjGDfM8uT6QcOLMdAExElv3GuwHbhqCQ2sl4irvr-6e0roXdS2PAQn9R-PoddEFBjdPRyb0iSrPpa20IfrpwPzbPSDtEE5c41wbJ-c9EXNA2fSWNCCuZDwxZewD2ZHmtVkA95QiL-mzwoc/s1600/c3.PNG" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="231" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjGDfM8uT6QcOLMdAExElv3GuwHbhqCQ2sl4irvr-6e0roXdS2PAQn9R-PoddEFBjdPRyb0iSrPpa20IfrpwPzbPSDtEE5c41wbJ-c9EXNA2fSWNCCuZDwxZewD2ZHmtVkA95QiL-mzwoc/s320/c3.PNG" width="320" /></a></div>
<div class="separator" style="clear: both; text-align: left;">
<br /></div>
<br />
<br />
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com1tag:blogger.com,1999:blog-4918390655818864374.post-3572150467889272612012-03-16T23:46:00.000+09:002012-07-24T01:08:56.636+09:00黒潮の流路をpandasで<br />
<br />
<a href="http://oceansciencehack.blogspot.com/2010/02/blog-post_1073.html" target="_blank">以前</a>、気象庁黒潮の流軸の緯度のデータをネットから読み込んでプロットしてみた。<br />
<br />
同じことを<a href="http://pandas.pydata.org/" target="_blank">pandas (Python Data Analysis Library)</a>で行ってみる。<br />
<br />
<pre class="prettyprint linenums">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")
</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjDh9_3Rskm8bhztl1Cf31qsyxhuei7hQ54-YJBZr2EVpN2avpGKRZekeutoJbYIRt6iJ0w26U9J0KXLeH8eE4fPxhkDBvMI6u4fjmvoQzuecvwWbUNvxb-ZwUJXnuaUUovYdKLtB7vQfg/s1600/pandas_example.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="300" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjDh9_3Rskm8bhztl1Cf31qsyxhuei7hQ54-YJBZr2EVpN2avpGKRZekeutoJbYIRt6iJ0w26U9J0KXLeH8eE4fPxhkDBvMI6u4fjmvoQzuecvwWbUNvxb-ZwUJXnuaUUovYdKLtB7vQfg/s400/pandas_example.png" width="400" /></a></div>
<br />
(注記) 上のデータ読み込みでは、元データがスペースを使っている関係で、"space"という余計なデータ列(column)を作ってしまう。もっとスマートな方法があるかもしれない。<br />
(追記:2012/7/23) <a href="http://oceansciencehack.blogspot.jp/2012/07/pandas-version-0.html" target="_blank">後日</a>、改良版を書いた。<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-89186221483022179392012-02-06T16:48:00.001+09:002012-02-06T16:49:44.463+09:00GrADSで作成するpostscriptの不要な横線GrADSでpostscript形式の図を作ると不要な横線が入ってしまうことがある。<br />
<br />
従来の対策は、例えば<br />
http://www.sci.hokudai.ac.jp/~sasakiyo/tips_grads.html<br />
<br />
GrADS 2.0.0から従来のgxout shaded に加えて、gxout shade2が加わっている。<br />
http://www.iges.org/grads/gadoc/gradcomdsetgxout.html<br />
これで図を作ったところ、従来横線が入っていた場合に、横線が入らなくなった。常にそうであるかは要確認。tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-750664125543185202011-10-17T08:51:00.002+09:002013-04-13T07:34:01.377+09:00フィルタの応答関数<br />
<br />
フィルタの周波数応答関数をプロットしてみる。<br />
<br />
以下の例は、<a href="http://www.kumst.kyoto-u.ac.jp/kougi/time_series/">京都大学の講義「MATLABで時系列解析」</a>の<a href="http://www.kumst.kyoto-u.ac.jp/kougi/time_series/ex1030.html">線形システムの伝達関数を求める</a>を翻案したものである。<br />
<br />
想定としては、1ヶ月ずつのデータがあり、シグナルとしては、48ヶ月周期と6ヶ月周期とホワイトノイズが重なったシグナルに対して13ヶ月の移動平均をかけている。<br />
<br />
理論的な周波数応答関数を求める関数は<a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.freqz.html#scipy.signal.freqz">scipyのfreqz</a>関数である(<a href="http://www.mathworks.co.jp/help/toolbox/signal/ref/freqz.html">matlabのfreqz</a>関数に対応する)。応答関数は入力 x のパワースペクトル密度関数 Pxx と出力のクロススペクトル密度関数 Pyx の比としても推定できる。スペクトル密度関数を求める関数psdに関しては<a href="http://www.google.co.jp/search?rlz=1C1GGGE_jaJP404JP405&gcx=c&sourceid=chrome&ie=UTF-8&q=mlab.psd+site%3Aoceansciencehack.blogspot.com">何度か紹介した</a>。<br />
<br />
<pre class="prettyprint linenums">#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as plab
import matplotlib.mlab as mlab
from scipy.signal import freqz,lfilter
#--- Signal
Fs = 1 # frequency
T = 1000. # final time
t = np.arange(0.,T+1./Fs,1./Fs) # time
x = np.sin(2.*np.pi*t/6.) + np.sin(2.*np.pi*t/48.)+np.random.normal(0,0.1,len(t)) # input signal
#---- power spectral by Welch's method
plt.figure()
nfft = 256
window = np.hanning(256)
noverlap = 128
[Pxx,f] = mlab.psd(x,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # Power spectral
plt.subplot(2,2,1)
plt.plot(f,10*np.log10(Pxx)) # plot in decibel
plt.legend(('Pxx',))
plt.plot([1./12.,1./12],[-70.,20.],'k--')
plt.plot([1./6.,1./6],[-70.,20.],'k--')
plt.plot([1./48.,1./48],[-70.,20.],'k--')
plt.ylim((-70,20))
#--- Moving average
h = np.ones(13)/13.
#--- Calculate theoritical response
fs,Res = freqz(h)
plt.subplot(2,2,2)
plt.plot(fs/2./np.pi,np.abs(Res))
plt.legend(('theoretical response',))
plt.plot([1./12.,1./12],[0.,1.],'k--')
plt.plot([1./6.,1./6],[0.,1.],'k--')
plt.plot([1./48.,1./48],[0.,1.],'k--')
plt.ylim((0,1.))
#--- Apply Filter
y = lfilter(h,1,x) # filter
[Pxy,f] = mlab.csd(x,y,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # cross scectral density
[Pyy,f] = mlab.psd(y,NFFT=nfft,Fs=Fs,window=window,noverlap=noverlap,detrend=plab.detrend_none) # power spectral density
plt.subplot(2,2,3)
plt.plot(f,10.*np.log10(abs(Pyy)))
plt.legend(('Pyy',))
plt.plot([1./12.,1./12],[-70.,20.],'k--')
plt.plot([1./6.,1./6],[-70.,20.],'k--')
plt.plot([1./48.,1./48],[-70.,20.],'k--')
plt.ylim((-70,20))
#--- Calculate estimated response
Resc = Pxy/Pxx
plt.subplot(2,2,4)
plt.plot(f,abs(Resc));
plt.legend(('estimated response',))
plt.plot([1./12.,1./12],[0.,1.],'k--')
plt.plot([1./6.,1./6],[0.,1.],'k--')
plt.plot([1./48.,1./48],[0.,1.],'k--')
plt.ylim([0,1])
#--- Save fig
plt.savefig("test_response.png")
plt.show()
</pre>
<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdWA_lyes8cR-fHN-5XaSbz1qmRjb_sHGutDYnrbeQBa_ikmZRgxhAwu6tHU_HhZcHOZ1TO8XNO0u9CPvwJj6qthpMqjGFNUd5B7NlkB3yAeR0vufJL1rdFEaWpHFsiwxJuv1RQCjgJ_8/s1600/test_response.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjdWA_lyes8cR-fHN-5XaSbz1qmRjb_sHGutDYnrbeQBa_ikmZRgxhAwu6tHU_HhZcHOZ1TO8XNO0u9CPvwJj6qthpMqjGFNUd5B7NlkB3yAeR0vufJL1rdFEaWpHFsiwxJuv1RQCjgJ_8/s320/test_response.png" width="320" /></a></div>
<br />
以下がフィルタを12ヶ月移動平均にした場合。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRr1r0vfWHQ3bDJwq0yP9oXKBPm2uCjRvfXpQVkFG4LUB-fvWmUzdhzNyb5RNoGe3oQlPg7yCoe10KNvr_3bLjYswB4xZz9Tii7CnuvCF9bTwJ8HdCyfg68IsdQOhISSLisEJWI7TPO8g/s1600/test_response12.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgRr1r0vfWHQ3bDJwq0yP9oXKBPm2uCjRvfXpQVkFG4LUB-fvWmUzdhzNyb5RNoGe3oQlPg7yCoe10KNvr_3bLjYswB4xZz9Tii7CnuvCF9bTwJ8HdCyfg68IsdQOhISSLisEJWI7TPO8g/s320/test_response12.png" width="320" /></a></div>
以下はシグナルからホワイトノイズを除いた場合。この時はスペクトルからの応答関数の推定はうまくいかない(移動平均はこの場合は13ヶ月)。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi2ZcYKcEHrIS3XXUFcSl-cptoiW4wYmZWfPhGhMrGm4qLaMq9lHWIj1cLnmyN2VSX9b09a_hHr8Eo4y8vzARl8nOy6uNAxs3nr8SNBmG-bg9tP7Dj2Z7MwDP6tKNqh4zmtY_aQM40IL_A/s1600/test_response_now.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="240" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi2ZcYKcEHrIS3XXUFcSl-cptoiW4wYmZWfPhGhMrGm4qLaMq9lHWIj1cLnmyN2VSX9b09a_hHr8Eo4y8vzARl8nOy6uNAxs3nr8SNBmG-bg9tP7Dj2Z7MwDP6tKNqh4zmtY_aQM40IL_A/s320/test_response_now.png" width="320" /></a></div>
<br />
<br />
<br />
参考) <a href="http://www.scipy.org/Cookbook/FIRFilter">Scipy Cookbook FIRFilter</a>tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-2759057713930583542011-05-06T18:27:00.000+09:002011-05-06T18:27:10.659+09:00気象庁全球数値予報モデルGPV (GSM)<br />
<br />
<a href="http://www.jmbsc.or.jp/hp/online/f-online0a.html">気象庁全球数値予報モデルGPV (GSM) </a>の初期値データを<a href="http://database.rish.kyoto-u.ac.jp/arch/jmadata/gpv-original.html">京都大学のサイト</a>からダウンロードして、<a href="http://shimizus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%E2%A5%C7%A5%EB%B4%D8%CF%A2%A5%E1%A5%E2">清水慎吾さん作成のgsm2bin</a>でgrib2形式からgrads形式に変換するモジュール。<br />
<br />
以下がそのモジュール。サンプルのメインプログラムでは、2011年1月1日6時のデータを変換し、<a href="http://opengrads.org/wiki/index.php?title=Python_Interface_to_GrADS">python interface to grads</a>で海面気圧をプロットしている。<br />
<br />
get_gsm.py<br />
<pre class="prettyprint">def get_gsm(date=[2011,1,1,0],dir="."):
"""
download jma gsm initial data from http://database.rish.kyoto-u.ac.jp/arch/jmadata/da
ta/gpv/original/
and
convert to grads file using gsm2bin by Dr. Shingo Shimizu
(http://shimizus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%E2%A5%C7%A5%EB%B4%D8%CF%A2%A5%E1%A
5%E2)
usage:
get_gsm(date,dir)
input
date: [year,month,day,hour] hour must be multiple of 6
default [2011,1,1,0]
dir : directory to save data
"""
import os,urllib
#check
assert date[3] % 6 ==0, "Error, hour must be multiple of 6!!!"
#names from date
urldir,original,gradsbinname,gradsdate=name_maker(date)
# get original data
if not os.path.exists(dir+"/"+original):
if not os.path.exists(dir): os.mkdir(dir)
# Download the data
print 'Downloading data, please wait'
opener = urllib.urlopen(urldir+original)
open(dir+"/"+original, 'w').write(opener.read())
print 'Downloaded'
# convert to grads file
os.system("gsm2bin %(dir)s/%(original)s %(dir)s/%(gradsbinname)s" %locals())
# make ctlfile
make_ctl(dir,gradsbinname,gradsdate)
############################################
def name_maker(date):
"""
make filenames from date
"""
import datetime
time=datetime.datetime(*date)
cyear=time.strftime("%Y")
cmonth=time.strftime("%m")
cmonth2=(time.strftime("%b")).lower()
cday=time.strftime("%d")
chour=time.strftime("%H")
# URL of original data
urldir="http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/%(cyear)s/%(
cmonth)s/%(cday)s/" %locals()
# File name of original data
original="Z__C_RJTD_%(cyear)s%(cmonth)s%(cday)s%(chour)s0000_GSM_GPV_Rgl_FD0000_grib2.
bin" %locals()
# grads bin name after conversion
gradsbinname="gsm%(cyear)s%(cmonth)s%(cday)s%(chour)s.bin" %locals()
# date in grads format
gradsdate="%(chour)s:00Z%(cday)s%(cmonth2)s%(cyear)s" %locals()
return urldir,original,gradsbinname,gradsdate
###########################################
def make_ctl(dir,gradsbinname,gradsdate):
"""
make grads ctl file for GSM data
"""
#ctl file text
ctltxt="""
dset ^%(gradsbinname)s
options template
undef -999
xdef 720 LINEAR 0.0 0.5
ydef 361 LINEAR -90.0 0.5
zdef 17 LEVELS 1000 925 850 700 600 500 400 300 250 200 150 100 70 50
30 20 10
tdef 1 LINEAR %(gradsdate)s 6hr
vars 16
slp 0 0 sea level pressure
sp 0 0 surface pressure
su 0 0 surface westerly wind
sv 0 0 surface southerly wind
st 0 0 surface temperature
srh 0 0 surface relative humidity
lca 0 0 LowCloudAmount
mca 0 0 MidCloudAmount
hca 0 0 HighCloudAmount
tca 0 0 TotalCloudAmount
z 17 0 height
u 17 0 westerly wind
v 17 0 southerly wind
t 17 0 temperature
rh 17 0 relative humidity
w 17 0 p-velocity
endvars
"""
# write to file
gradsctlname=dir+"/"+gradsbinname.replace(".bin",".ctl")
f=open(gradsctlname,"w")
f.write(ctltxt %locals())
f.close()
#-----------------------------------------------------------------------
if __name__ == '__main__':
# get gsm data
get_gsm(date=[2011,1,1,6],dir="DATA")
# plot
from grads.galab import GaLab # for python interface for grads
ga = GaLab(Bin='grads',Window=False,Echo=True,Verb=0,Port=False)
ga.open("DATA/gsm2011010106.ctl")
script="""
set grads off
set gxout shaded
d slp/100
cbarn
draw title Sea level pressure (hpa)
printim gsm.png white
"""
ga(script)
</pre>
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhQWRZyymJEigUzFQYk7k0Try5tL2avqiTeF1L9i3JukTPjrD2Fl9cwMNAZOXw7gITYcgnV-boB7cNVT2ulY64sWMwbYcSXURN5Ep5zh_GfR9HFcOQkI5z5ptxl4fTUdxteq5TesCo6zUQ/s1600/gsm.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="247" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhQWRZyymJEigUzFQYk7k0Try5tL2avqiTeF1L9i3JukTPjrD2Fl9cwMNAZOXw7gITYcgnV-boB7cNVT2ulY64sWMwbYcSXURN5Ep5zh_GfR9HFcOQkI5z5ptxl4fTUdxteq5TesCo6zUQ/s320/gsm.png" width="320" /></a></div>
<br />
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0tag:blogger.com,1999:blog-4918390655818864374.post-43088170526134482442011-03-05T09:52:00.000+09:002011-03-05T09:52:13.647+09:00気象庁MGDSST<br />
<br />
気象庁の海面水温データ(0.25x0.25度、daily)であるMGDSSTをウェブからダンロードしてから、データを読み込むまでのモジュールをつくる。MGDSSTは<a href="http://goos.kishou.go.jp/">NEAR-GOOSのサイト</a>から入手が可能である。<br />
idとpasswordが必要なので、あらかじめ登録しておく。<br />
<br />
以下がスクリプト<br />
usrname="xxxxx"<br />
password="*****"<br />
のところを自分のidとpasswordで置き換える。<br />
<br />
<pre class="prettyprint">import numpy as np
import os
def downloadfile(infile="remotefile",outfile="localfile"):
"""
download file from NEAR-GOOS data site with id and password
input
infile: remote file on the web
outfile: local file
set password and id below
"""
import urllib2
host="goos.kishou.go.jp:80"
realm="NEAR-GOOS"
usrname="xxxxx"
password="*****"
auth_handler = urllib2.HTTPBasicAuthHandler()
auth_handler.add_password(realm, host, usrname, password)
opener = urllib2.build_opener(auth_handler)
urllib2.install_opener(opener)
f=urllib2.urlopen(infile)
open(outfile, 'w').write(f.read())
#----------------------------------------------------------
def readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)):
"""
load MGDSST data from NEAR-GOOSE site
read data with fortran (f2py)
select region specified by lonrange and latrange
Usage:
lon,lat,sst= readMGDSST(dir=".",year=2010,month=1,day=1,lonrange=(0,360),latrange=(-90,90)
Input
dir: directory for data download
year,month,day: date of the data
lonrange,latrange = longitude and latitude range for data
output
lon, lat: longitude, latitude
sst: sst data
"""
# date
import datetime
t=datetime.datetime(year,month,day)
syear=t.strftime("%Y")
smon=t.strftime("%b")
sday=t.strftime("%d")
# filename
dataname="mgdsst."+smon+sday
filename=dir+"/"+dataname+"."+syear
filenamegz=filename+".gz"
print filename
urlname="http://goos.kishou.go.jp/rrtdb/usr/pub/JMA/mgdsst/"+syear+"/"+dataname+".gz"
# file exist?
if not os.path.exists(filename):
if not os.path.exists(filenamegz):
downloadfile(infile=urlname,outfile=filenamegz)
os.system("gunzip %(filenamegz)s" %locals())
# temp.bin is used for temporary file
if os.path.lexists("temp.bin"):
raise IOError("temporary file temp.bin somehow exist! Remove it!")
else:
os.system("ln -s %(filename)s temp.bin" %locals())
# data read through fortran (f2py)
from read_sst import read_sst
im=1440
jm=720
sst=np.zeros([jm,im])
sst,ierr=read_sst()
assert ierr==0, "data error in fortran"
lon=0.125+0.25*np.arange(im)
lat=-89.875+0.25*np.arange(jm)
os.system("rm temp.bin")
# range selection
xrange=np.logical_and(lon>=lonrange[0],lon<=lonrange[1])
yrange=np.logical_and(lat>=latrange[0],lat<=latrange[1])
sst=(sst[yrange,:])[:,xrange]
return lon[xrange],lat[yrange], sst
#-----------------------------------------------------------------------
if __name__ == '__main__':
import cdms2,MV2,vcs #cdat
# read data
lon,lat,sst=readMGDSST(dir="DATA",year=2011,month=1,day=1,lonrange=(0,360),latrange=(-90,90))
# hearafter plot with cdat
# First let's create a variable
var = MV2.masked_array(sst,mask=(sst<-990.),id="SST")
# Now create the axis
lons=cdms2.createAxis(lon)
lats=cdms2.createAxis(lat)
# And decorate it
lons.id = 'longitude'
lats.id = 'latitude'
lons.units = 'degrees_east'
lats.units = 'degrees_north'
# Just to make sure we will describe this axis as a longitude one
lons.designateLongitude()
lats.designateLatitude()
# Now let's apply these to the variable
var.setAxisList((lats,lons))
#var.info() to see info
#plot
v = vcs.init()
v.plot(var)
v.png("MGDSST.png")
del v
</pre>
<br />
上のサンプルメインプログラムでは、2011年1月1日のデータをディレクトリ"DATA"にダウンロードし、経度0から360度、緯度-90度から90度の範囲で読み込んでいる。<br />
その上で、<a href="http://www2-pcmdi.llnl.gov/cdat">cdat</a>を用いて可視化した。<br />
結果は以下の図。<br />
<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgqM072hiUbRG8aQazBgGLjGMoWK8X4T6A3XbZ6QcYp5QBBX0Hn_Ve6Gvliz_WjplfmsgloV1yzzVjJFrGL8iS3bffWR1B_ZYwV3S0XKrhehJzonr7683kxcsmXpiin9re4NZoZeIAT0jY/s1600/MGDSST.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" height="308" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgqM072hiUbRG8aQazBgGLjGMoWK8X4T6A3XbZ6QcYp5QBBX0Hn_Ve6Gvliz_WjplfmsgloV1yzzVjJFrGL8iS3bffWR1B_ZYwV3S0XKrhehJzonr7683kxcsmXpiin9re4NZoZeIAT0jY/s400/MGDSST.png" width="400" /></a></div>
<br />
ダウンロードしたバイナリファイルを読むのにfortranのプログラムをモジュールとしてつかう。<a href="http://shimizus.hustle.ne.jp/wiki/wiki.cgi?page=%A5%E2%A5%C7%A5%EB%B4%D8%CF%A2%A5%E1%A5%E2">清水慎吾さんのデータを読むサンプルプログラム</a>を参考に、以下のようなfortranサブルーチンを作った(配列はpython風に緯度と経度を入れ替えている。欠損値は-999.)。<br />
このプログラムを<a href="http://cens.ioc.ee/projects/f2py2e/">f2py</a>を使って以下のようにコンパイルすればpythonのモジュールとして使用が可能になる(上のプログラムのread_sst)。<br />
<span class="Apple-style-span" style="background-color: #f3f3f3;">f2py --f90exec=gfortran -c -m read_sst read_sst.f90</span><br />
(ここではgfortranを使っている。)<br />
<br />
read_sst.f90<br />
<pre class="prettyprint">subroutine READ_SST(sst,ierr)
implicit none
integer :: year,mon,day
integer :: i,j,jj
integer,parameter :: im=1440
integer,parameter :: jm=720
!
real,intent(out):: sst(jm,im)
integer, intent(out):: ierr
!
integer :: isst(im,jm),ist
real :: undef = -999.
open(55,file="temp.bin",status='old',form='formatted',iostat=ierr,err=999)
read(55,*) year,mon,day
write(6,*) year,mon,day
do j=1,jm
read(55,10,iostat=ierr,err=999) (isst(i,j),i=1,im)
enddo
close(55)
10 format(1440i3)
do j=1,jm
do i=1,im
jj = jm -j +1
ist = isst(i,jj)
if(ist .eq. 888 .or. ist.eq.999) then
sst(j,i) = undef
else
sst(j,i) = real(ist)*0.1
endif
enddo
enddo
return
999 print *, "Read Error!"
end subroutine READ_SST
</pre>
<br />tmiyamahttp://www.blogger.com/profile/11413252495514096849noreply@blogger.com0