DIAMラジオ

2012 年 9 月 18 – 19 日 (千葉大・CEReS)
講師 : 眞子 直弘，カトリ プラディープ，齋藤 尚子, 樋口 篤志
第 6 回 VL 講習会@千葉大・CEReS
∼ エアロゾルと大気補正 ∼
演習 B テキスト
衛星画像の大気補正
担当： 眞子 直弘
i
目次
第 1 章 はじめに
1.1 テキスト本文中の表記 . . .
1.2 演習の準備 . . . . . . . . .
1.2.1 PC の立ち上げ . . .
1.2.2 演習ファイルの準備
1.3 演習の進め方 . . . . . . . .
1.4 演習ファイルの構成 . . . .
.
.
.
.
.
.
1
1
2
2
2
2
2
.
.
.
.
5
5
6
6
8
第 3 章 スカイラジオメーターのデータ処理
3.1 エアロゾル光学特性データの読み込み . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 エアロゾル物理特性データの読み込み . . . . . . . . . . . . . . . . . . . . . . . . . . . .
9
9
14
第 4 章 MODIS 衛星画像のデータ処理
4.1 MODIS とは . . . . . . . . . . .
4.2 MODIS L1B データ処理 . . . . .
4.2.1 HDF ファイルの閲覧 . . .
4.2.2 HDF ファイルの読み込み
4.2.3 地理位置情報の読み込み .
4.2.4 地図投影法の変換 . . . . .
第2章
2.1
2.2
2.3
2.4
第5章
5.1
5.2
5.3
5.4
Python の紹介
Python とは . . . . . . . .
はじめに . . . . . . . . . .
簡単な使用例 . . . . . . .
Python スクリプトの作成
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
15
15
18
18
20
24
25
6S を使った MODIS 衛星画像の大気補正
大気補正とは . . . . . . . . . . . . . . . .
6S とは . . . . . . . . . . . . . . . . . . .
6S の実行 . . . . . . . . . . . . . . . . . .
MODIS 衛星画像の大気補正 . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
31
31
32
33
39
.
.
.
.
.
.
.
.
43
43
43
44
45
45
48
48
49
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
付 録 A Python プログラミング簡易リファレンス
A.1 配列 . . . . . . . . . . . . . . . . . . . . .
A.1.1 リスト . . . . . . . . . . . . . . . .
A.1.2 タプル . . . . . . . . . . . . . . . .
A.1.3 文字列 . . . . . . . . . . . . . . . .
A.1.4 array (numpy.ndarray) . . . . . . .
A.1.5 辞書 . . . . . . . . . . . . . . . . .
A.2 制御構造 . . . . . . . . . . . . . . . . . . .
A.3 関数 . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
ii
例外処理 . . . . . . . . . . . . .
正規表現 (Regular Expression)
日付・時刻 . . . . . . . . . . .
数値データ処理 . . . . . . . . .
A.7.1 1 次元補間 . . . . . . . .
A.7.2 2 次元補間 . . . . . . . .
A.7.3 多項式最適化 . . . . . .
A.7.4 一般式最適化 . . . . . .
A.8 ファイル入力 . . . . . . . . . .
A.8.1 テキストファイル . . . .
A.8.2 バイナリファイル . . . .
A.8.3 画像ファイル . . . . . .
A.8.4 HDF ファイル . . . . .
A.8.5 NetCDF ファイル . . .
A.9 描画ライブラリ . . . . . . . . .
A.9.1 XY グラフ . . . . . . . .
A.9.2 散布図 . . . . . . . . . .
A.9.3 棒グラフ . . . . . . . . .
A.9.4 ヒストグラム . . . . . .
A.9.5 2 次元プロット . . . . .
A.9.6 地図投影 . . . . . . . . .
A.4
A.5
A.6
A.7
付 録 B 本演習に必要なソフトウェア
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
50
50
51
52
52
53
54
54
55
55
56
56
57
57
57
57
58
59
59
59
60
61
1
第 1 章 はじめに
本演習では、演習 A で得られたエアロゾル特性を使って MODIS 衛星画像の大気補正を行うことによ
り、エアロゾルの影響を実感しながら衛星データ解析の基礎を習得してもらう。
1.1
テキスト本文中の表記
演習の実行環境は Linux である。多くの作業はターミナルソフトにコマンドを入力して行うが、本テ
キストではシェル (bash) に対する入力を以下のように記述する。
bash に対する入力例
directory-name> some-command # some-comment
上の例で、directory-name、some-command にはそれぞれ作業ディレクトリ名および入力するコマンド
が入る。特に断らない限り、#以降の部分は説明用のコメントであり、実際に入力する必要はない。
本演習では、テキスト・バイナリデータ処理、数値演算、グラフ描画等を行うために IPython という
インタラクティブシェルを利用する。IPython への入力は以下のように記述する。
IPython に対する入力例
In : some-command # some-comment
Out: some-text
実際の IPython の表示は
In [1]: 1+1
Out[1]: 2
のようになるが、括弧内の番号は操作の手順によって変わり得るので省略する。ここでも some-command
は入力するコマンドであり、特に断らない限り、#以降の部分は説明用のコメントである。Out の部分は
IPython が出力するメッセージである。コマンドによっては何もメッセージを出力しない場合もある。
このように、In:というプロンプトを付けるのは IPython への入力であることを明示するためである
が、本テキストから直接コピー＆ペーストする場合には邪魔になるので、第 3 章以降の囲み記事では省
くことにする。その場合、囲み記事のタイトルに「(IPython 入力)」と明示するとともに、同等の動作
をするスクリプト名を書いておく。例えば、囲み記事のタイトルが
タイトル (IPython 入力 or some-script.py)
となっている場合、その内容を IPython に打ち込んでも良いし、some-script.py を実行しても良い。
(注：囲み記事ではそれ以前の内容を省いている。そのため、途中で IPython を再起動したりして、必要
なモジュールが読み込まれていなかったり、変数等が定義されていない場合は IPython に入力しても正
常に動作しないことがある。その場合は対応するスクリプトを実行すると良い。) 本演習は、後述する
「演習ファイルの準備」によって作られる~/work/2012VL_PartB/data ディレクトリで行うことを想定
している。一方、あらかじめ用意されたスクリプトは~/work/2012VL_PartB/script ディレクトリに置
いてある。このように、カレントディレクトリにないスクリプトを IPython で実行する場合、下の例の
ようにスクリプトのパスを指定する必要がある。
カレントディレクトリにないスクリプトの実行例
In : %run ../script/some-script.py
第1章
2
1.2
1.2.1
はじめに
演習の準備
PC の立ち上げ
計算機演習室の PC は Windows と Linux (CentOS) のデュアルブート方式になっている。PC の電源
スイッチを入れて暫くすると「DHCP..\」と表示され、その後「Press F8 to view menu.」と表示され
る。この時に F8 キーを押して OS 選択画面を表示させ、Linux を選択して立ち上げる。(注：Press F8
· · · のメッセージが表示されているのは 10 秒ほどで、それを見逃すと Windows が立ち上がってしまう。
Linux の選択に失敗した場合は電源ボタンを素早く (1 秒以下) 押して PC をシャットダウンし、しばら
く待ってから再起動する。) CentOS が立ち上がったら、各自に割り当てられたユーザー名とパスワー
ドを使ってログインする。
1.2.2
演習ファイルの準備
本演習に必要なファイルは、配布した USB メモリ中の 2012VL_PartB というディレクトリに入って
いる。まず、メニューバーのアプリケーション > アクセサリから GNOME 端末を立ち上げる。デフォル
トで tcsh というシェル環境になっているので、環境変数 PATH、LD_LIBRARY_PATH、PYTHONHOME、PS1
の設定をしてから bash というシェル環境に切り替える。2012VL_PartB.tgz をホームディレクトリに
展開し、data ディレクトリに HDF ファイルへのリンクを作成する。(注：計算機演習室の PC では、
ユーザーのホームディレクトリに 2GB の容量制限 (quota) が設けられている。HDF ファイルをホーム
ディレクトリにコピーすると容量の上限を越えてしまい、最悪の場合ログイン出来なくなる。ここでは
コピーは行わず、USB メモリ上にある HDF データへのシンボリックリンクを作成する。)
演習ファイルの準備
some-directory> cd # ホームディレクトリに移動 (別の場所に居た場合)
~> echo export PATH=.:${HOME}/local/bin:/bin:/usr/bin:/usr/local/bin >>.bashrc
~> echo export LD_LIBRARY_PATH=${HOME}/local/lib >>.bashrc
~> echo export PYTHONHOME=${HOME}/local >>.bashrc
~> echo export PS1=\"\\w\> \" >>.bashrc
~> bash
~> tar xzf /media/disk/2012VL_PartB/2012VL_PartB.tgz # 数分待つ
~> cd work/2012VL_PartB/data
~/work/2012VL_PartB/data> ln -s /media/disk/2012VL_PartB/*.hdf .
1.3
演習の進め方
本演習は主に学部 4 年生、大学院生、および若手研究者を対象としているが、演習内容に関する前提
知識は人それぞれであると考えられる。基本的には Linux 初心者に合わせて進めていく予定であるが、
余裕がある人は各自のペースで演習を進め、各章にある演習問題に挑んでほしい。
1.4
演習ファイルの構成
演習ファイルの構成は以下のようになっている。(前述のように、data ディレクトリにある HDF ファ
イルは USB メモリ上にある HDF ファイルへのシンボリックリンクである。)
1.4. 演習ファイルの構成
演習ファイルの構成
.
|-- local
‘-- work
3
‘-- 2012VL_PartB
|-- data
|
|
|
|-- 12051900opticalL2R01.TXT
|-- 12051900physicalL2R01.TXT
|-- MOD021KM.J20120519011621.20120519012656.hdf
|
|
|-- MOD02HKM.J20120519011621.20120519012656.hdf
|-- MOD02QKM.J20120519011621.20120519012656.hdf
|
|
|
|-- MOD03.J20120519011621.20120519012656.hdf
|-- MOD09GA.A2012140.h29v05.005.2012142054114.hdf
|-- input.dat
|
‘-- input2.dat
|-- doc
|
|
|
|-- 6S
|-- MODIS
‘-- VL2012_ChibaUniversity_PartB.pdf
‘-- script
|-- example.py
|-- modis_basemap.py
|-- modis_geolocation.py
|-- modis_plot_1km.py
|-- modis_plot_hkm.py
|-- modis_plot_qkm.py
|-- modis_read_1km.py
|-- modis_reprojection.py
|-- modis_sinusoidal.py
|-- sixs_atmos_correct.py
|-- sixs_calc_coeff.py
|-- sixs_lobs_vs_rho.py
|-- sixs_make_lut.py
|-- sixs_plot_coeff.py
|-- sixs_run_custom.py
|-- sixs_run_simple.py
|-- sixs_tau_vs_rho.py
|-- skyrad_aot_vs_time.py
|-- skyrad_aot_vs_wlen.py
|-- skyrad_read_opt.py
|-- skyrad_read_radius.py
‘-- skyrad_read_wlen.py
5
第 2 章 Python の紹介
本演習ではデータ処理を行うためのプログラミング言語として、主に Python (パイソン) を使用する。
実際にスカイラジオメーターや衛星データの解析を行う前に、ここで Python の基本的な使い方を学ぶ。
2.1
Python とは
Python とはオランダ人の Guido van Rossum 氏が作ったプログラミング言語の一種であり、オブジェ
クト指向を取り入れたインタプリタ型スクリプト言語である。C 言語や Fortran のようなコンパイラ言
語に比べ、少ない記述で同等の処理が可能であるために新規ソースの作成および既存ソースの解読が容
易 、コンパイルが不要であるために手軽に実行可能といった長所がある。ただし、実行速度が問題にな
る場合は C 言語や Fortran のようなコンパイラ言語を使った方が良く、場合によって使い分けるのが良
い。実際、実行速度が問題になる Python ライブラリの多くは C 言語や Fortran で実装されており、イ
ンタプリタ型スクリプト言語の中では高速に動作する。これらライブラリの多くは「その道のプロ」に
よって書かれており、下手に自分でコードを書くよりはライブラリを使った方が効率の良いプログラム
が出来る。
Python はとにかくライブラリが充実しており、テキスト・バイナリデータ処理、数値演算、グラフ
描画等、必要なデータ解析プログラムが Python だけで簡単に作れることが利点である。 プログラムが
Python だけで書けるのでプログラムに一貫性が保たれる。また、オーソドックスなシェルスクリプト
で多用されるファイルやパイプを介したプログラム間のデータ入出力が不要であるために処理が直接的
で見通しが良く、高速に動作するプログラムを作成することができる。さらに、Python を使った解析ス
キームが他のオーソドックスな解析スキームより有利な点として、プログラム開発やデバッグの容易さ
が挙げられる。IPython というインタラクティブシェルを利用すると、プログラム実行後に変数の中身
を確認したり、値を変更したりといったことが自由に行える。 そのため、ソースを書き換えてファイル
や画面に書き出して見ないと変数の中身が分からない C 言語や Fortran 等に比べてプログラム開発やデ
バッグが効率良く行える。IPython のオンラインヘルプ機能や入力補完機能もプログラム開発の効率化
に役立つ。
現在使われている Python にはバージョン 2 とバージョン 3 の 2 系統あるが、バージョン 3 は理想を
追求するために後方互換性を犠牲にして策定された仕様を含んでいる。今のところバージョン 3 に対応
していないライブラリも少なくないので、本演習ではバージョン 2 を使用する。現時点でバージョン 2
の最新版は 2.7.3 である。付録 A に Python プログラミングの簡易リファレンスがあるので、必要に応
じて参照してもらいたい。Python に関する豊富な資料がオンラインで容易に入手できるが、特に以下
の Web ページはお勧めである。
• Python チュートリアル
http://www.python.jp/doc/release/tutorial/index.html
• Python ライブラリリファレンス
http://www.python.jp/doc/release/library/index.html
• Numpy (数値計算ライブラリ) チュートリアル
http://www.scipy.org/Tentative NumPy Tutorial
• Matplotlib (描画ライブラリ) ギャラリー
http://matplotlib.sourceforge.net/gallery.html
• Python Scientific Lecture Notes
http://scipy-lectures.github.com/
第 2 章 Python の紹介
6
2.2
はじめに
ここでは Python プログラムを IPython というインタラクティブシェル上で実行する。IPython は普
通の Python インタラクティブシェルに比べてヘルプ機能や補完機能がある等、使いやすさが強化され
ている。
IPython の起動
~/work/2012VL_PartB/data> ipython --pylab
--pylab は Matplotlib ライブラリを IPython 上で簡単に使えるようにするためのオプションである。
IPython が立ち上がったら、まず最初にお決まりの”Hello World!”の出力をやってみよう。
Hello World!
In : print "Hello World!"
Hello World!
Python では文字列リテラルを作るのにダブルクォーテーション (”) とシングルクォーテーション (’) の
両方が使える (Ruby のように文字列展開に関する違いはない)。普段はシングルクォーテーションを使
うことも多いが、そうすると本テキストからコマンドラインにコピー＆ペーストした時に正常に動作し
ない場合があるので、ここではダブルクォーテーションを使う。
IPython は Ctrl+d を 2 回押すか、quit() と入力して終了できる。
2.3
簡単な使用例
さて、いきなりではあるが、簡単なので、Python の使用例として三角関数のグラフを作成してみよう。
三角関数グラフの作成
In : x = arange(0.0,pi,0.01) # 0<=x<πの範囲で 0.01 刻みの array を作成する
In : x # x の中身を確認
Out:
array([ 0.
,
0.01,
0.02,
0.03,
0.04,
0.05,
0.06,
0.07,
0.08,
0.09,
0.1 ,
0.11,
0.12,
0.13,
0.14,
0.15,
0.16,
0.17,
3.06,
3.07,
3.08,
3.09,
3.1 ,
3.11,
3.12,
3.13,
3.14])
(中略)
In : x.shape # x のサイズを確認
Out: (315,)
In : y = sin(x)
In : y # y の中身を確認
Out:
array([ 0.
,
0.04997917,
0.00999983,
0.05996401,
0.01999867,
0.06994285,
0.0299955 ,
0.07991469,
0.03998933,
0.08987855,
0.0315874 ,
0.02159098,
0.01159239,
0.00159265])
(中略)
0.04158066,
In : plot(x,y)
Out: [<matplotlib.lines.Line2D at 0xa931d0c>] # at の後のアドレスは実行ごとに変わる
In : plot(x,cos(x)) # このように直接 cos(x) をプロットすることも出来る
Out: [<matplotlib.lines.Line2D at 0xace8dcc>]
In : savefig("example.png") # 拡張子から画像形式を自動的に判断して保存する
In : !display example.png & # ビックリマーク (!) を使うとシェルコマンドを呼び出せる
2.3. 簡単な使用例
7
1.0
0.5
0.0
0.5
1.00.0
0.5
1.5
1.0
2.0
2.5
図 2.1: 三角関数のグラフ
3.0
3.5
上のようなグラフが得られただろうか？このように、Python の使用法は簡単で、直感的に分かりやすい。
それでは、もう少し IPython の便利な機能を見てみよう。
IPython の便利な機能
In : x.shape # 上矢印 (↑) を 7 回押して表示させる (履歴機能)
Out: (315,)
In : x.shape # x.s まで入力した後、上矢印 (↑) を 1 回押して表示させる (履歴機能)
Out: (315,)
In : x.shape # x.sh まで入力した後、Tab を入力して表示させる (入力補完機能)
Out: (315,)
In : x. # x. まで入力した後、Tab を入力すると入力候補リストが表示される (入力補完機能)
x.T
x.all
x.copy
x.ctypes
x.getfield
x.imag
x.put
x.ravel
(中略)
x.conjugate
x.flatten
x.ptp
x.sort
x.squeeze
x.std
In : x.std? # オンラインヘルプ機能
Type:
builtin_function_or_method
(中略)
a.std(axis=None, dtype=None, out=None, ddof=0)
Returns the standard deviation of the array elements along given axis.
Refer to ‘numpy.std‘ for full documentation.
(後略)
このように、入力補完機能を使えば使用中のオブジェクトにどういったメソッド・属性があるかが分か
る。(Python では ID、型、値を持つデータをオブジェクトと呼び、それに関連付けて定義された関数、
変数をそれぞれメソッド、属性と言う。上の例では、x は array オブジェクト、std、shape はそれぞれ
x のメソッドおよび属性である。) 多くのメソッド・属性は self explaining になっており、名前を見れば
何に使うかが大体分かる。さらにメソッドの引数等の詳細が知りたい場合にはオンラインヘルプ機能が
役立つ。(上の例では IPython の出力が分かりやすいように青文字にしているが、実際にそのように表
示されるわけではない。)
第 2 章 Python の紹介
8
2.4
Python スクリプトの作成
これまでは Python コマンドを IPython に逐次入力しながら実行してきたが、これだと再度同じこと
をしたい場合に面倒であるし、時間がかかる。そこで、これらのコマンドをファイルに書いたもの (ス
クリプト) を用意しよう。スクリプトを使えばたくさんのコマンドを一度に実行することもできる。基
本的に IPython に入力した内容と全く同じで良いが、以下の点に注意する。
• --pylab オプションと同じ効果を得るには matplotlib.pylab モジュールの import が必要。
• IPython が自動で有効にしているインタラクティブモードの設定が必要。
• IPython 独自のコマンド (!、%、?等で呼び出すもの) はスクリプトの中では使えない。
vi や gedit 等、好みのエディタを使って下記内容のスクリプトを作成し、example.py というファイル
名で保存しよう。
三角関数グラフ作成スクリプト (example.py)
from matplotlib.pylab import *
interactive(True)
x = arange(0.0,pi,0.01)
x
y = sin(x)
y
plot(x,y)
plot(x,cos(x))
savefig("example.png")
ファイルに保存されたスクリプトを IPython 上で実行するには、%run というマジック関数を利用する。
(%で始まる IPython 独自の関数をマジック関数と呼ぶ。他に%cpaste 等がある。)
三角関数グラフ作成スクリプトの実行
In : %run example.py
In : x # x の中身を確認
Out:
array([ 0.
,
0.01,
0.02,
0.03,
0.04,
0.05,
0.06,
0.07,
0.08,
0.09,
0.1 ,
0.11,
0.12,
0.13,
0.14,
0.15,
0.16,
0.17,
3.06,
3.07,
3.08,
3.09,
3.1 ,
3.11,
3.12,
3.13,
3.14])
(中略)
In : y # y の中身を確認
Out:
array([ 0.
,
0.04997917,
(中略)
0.00999983,
0.05996401,
0.01999867,
0.06994285,
0.0299955 ,
0.07991469,
0.03998933,
0.08987855,
0.04158066, 0.0315874 , 0.02159098, 0.01159239, 0.00159265])
上の例のように、IPython を使えばスクリプト実行後に変数の中身を確認することが出来る。この機能
はスクリプトのデバッグ等に便利である。
9
第 3 章 スカイラジオメーターのデータ処理
ここでは、演習 A でスカイラジオメーターのデータを解析して得られたエアロゾル光学特性・物理特
性のデータを扱う。
3.1
エアロゾル光学特性データの読み込み
スカイラジオメーターで得られたエアロゾル光学特性のデータを Python で読み込んでみよう。この
データは以下のような形式になっている。(データを gedit で開き、編集 > 設定 画面で「行番号を表示
する」にチェックを入れ、「テキストの折り返しを有効にする」のチェックを外すと見やすくなる。)
エアロゾル光学特性データ (12051900opticalL2R01.TXT)
































ヘッダー 

39 行




































データ 

104 行 



========================================= · · ·
Aerosol optical thickness, single scatter · · ·
========================================= · · ·
Observation area · · ·
---------------- · · ·
Serial number of · · ·
---------------- · · ·
NOTATIONS:
(中略)
---------------------------------------------------------------- · · ·
Level
Year
Month
Day
Hour(LT)
Jday(UT)
SZA
···
···
---------------------------------------------------------------- · · ·
L2.20
L2.20
2012
2012
5
5
19
19
(中略)
L2.20
|
2012
5
19
5.740
5.830
15.500
{z
データ 44 列
139.86417
139.86792
76.790 · · ·
75.728 · · ·
140.27083
53.277 · · ·
まず、細かいことにはこだわらずに最短のコードを示す。
エアロゾル光学特性データの読み込み+描画 (IPython 入力 or skyrad_read_opt.py)
}
data = loadtxt("12051900opticalL2R01.TXT",skiprows=39,usecols=range(1,44))
plot(data[:,3],data[:,10])
これで図 3.1 のような図が得られる。このように、数値テーブル形式のテキストデータは loadtxt
というメソッドを使えば一発で読み込める。この例では、skiprows で先頭 39 行を読み飛ばしている。
また、0 列目には"L2.20"という数値に変換出来ない文字列が入っているので、usecols で 1∼43 列目
(1 ≤ 列番号 < 44) だけ読み込んでいる。読み込まれたデータは data という 104×43 の 2 次元 array に
格納される。各列の内容は表 3.1 のようになっている。
第 3 章 スカイラジオメーターのデータ処理
10
0.40
0.35
0.30
0.25
0.20
0.154
14
6
8
10
12
図 3.1: エアロゾル光学的厚さの時系列変化
16
表 3.1: 読み込んだエアロゾル光学特性データ (data) の内容
0
Year
1
Month
2
Day
3
Hour(LT)
4
Jday(UT)
5
SZA
6
AOT(Ch. 1)
7
AOT(Ch. 2)
8
AOT(Ch. 3)
9
AOT(Ch. 4)
10
AOT(Ch. 5)
11
AOT(Ch. 6)
12
AOT(Ch. 7)
13
AOT(Ch. 8)
14
AOT(Ch. 9)
15
AOT(Ch.10)
16
AOT(Ch.11)
17
Ang. Exp.
18
Beta(0.5 µm)
19
Irrad.(Ave.)
20
Irrad.(S.D.)
21
SSA(Ch. 1)
22
SSA(Ch. 2)
23
SSA(Ch. 3)
24
SSA(Ch. 4)
25
SSA(Ch. 5)
26
SSA(Ch. 6)
27
SSA(Ch. 7)
28
SSA(Ch. 8)
29
SSA(Ch. 9)
30
SSA(Ch.10)
31
SSA(Ch.11)
32
ASY(Ch. 1)
33
ASY(Ch. 2)
34
ASY(Ch. 3)
35
ASY(Ch. 4)
36
ASY(Ch. 5)
37
ASY(Ch. 6)
38
ASY(Ch. 7)
39
ASY(Ch. 8)
40
ASY(Ch. 9)
41
ASY(Ch.10)
42
ASY(Ch.11)
それでは、少し図を改善してみよう。
エアロゾル光学特性データの描画の改善 (IPython 入力 or skyrad_aot_vs_time.py)
NCH = 11
aot = data[:,arange(NCH)+6] # select column 6, 7, ..., 16
figure() # open new canvas
grid() # draw grids
minorticks_on() # draw small tick marks
for i in range(NCH):
if all(aot[:,i] < 0.0): # -9.99999 -> Data not available
continue
plot(data[:,3],aot[:,i],"o",label="Ch.%d"%(i+1),color=cm.jet(float(i)/NCH))
legend(loc=2,numpoints=1,frameon=False)
xlabel("Hour")
ylabel("Aerosol Optical Thickness")
3.1. エアロゾル光学特性データの読み込み
11
ここでは data の 6∼16 列目を取り出して aot という array を作っている。このスカイラジオメー
ターは全部で NCH = 11 のチャンネルを持っているが、描画の際、光学的厚さが全て負の値になってい
るチャンネルは使えないのでスキップしている。各データ点を線でつなぐと見栄えが悪いので、"o"オ
プションを指定して丸いマーカーでプロットしている。各チャンネルに対応するマーカーは cm という
モジュールの jet というカラーマップを使って色分けしている。引数は 0 から 1 までの実数値であり、
0 が青系、1 が赤系の色に対応する。
ここで for 文や if 文といった制御文が使われていることに注意しよう。Python は C 言語や他の多
くのプログラミング言語と異なり、for 文や if 文などのブロック区切りに括弧を使わず、インデントを
使う。(詳細は付録 A.2 参照) このようなインデントを含むコードを IPython にコピー＆ペーストで入力
する場合は、%cpaste というマジック関数を使うと良い。ペーストが終わったら最後に-- だけの 1 行を
入力してペーストを終了する。
cpaste を使った IPython への入力
In: %cpaste
Pasting code; enter ’--’ alone on the line to stop or use Ctrl-D.
:NCH = 11
:aot = data[:,arange(NCH)+6] # select column 6, 7, ..., 16
(中略)
:xlabel("Hour")
:ylabel("Aerosol Optical Thickness")
:-Out: <matplotlib.text.Text at 0x9d5b4ac>
これで図 3.2 のような図が得られる。ラベルや凡例の位置など、まだ細かい改善の余地があるが、と
りあえずこれで良しとしよう。
0.6
Aerosol Optical Thickness
0.5
0.4
Ch.2
Ch.3
Ch.4
Ch.5
Ch.6
Ch.7
Ch.9
0.3
0.2
0.1
0.04
14
10
12
Hour
図 3.2: エアロゾル光学的厚さの時系列変化
6
8
16
演習問題 1 スカイラジオメーターのエアロゾル光学特性データには、光学的厚さ (AOT) の他
にも単一散乱アルベド (SSA) や非対称性パラメータ (ASY) といったパラメータが入っている。上
の例を参考に、単一散乱アルベド、非対称性パラメータの時系列変化図を作成してみよう。
ヒント: SSA は data の 21∼31 列目、ASY は 32∼42 列目である。
第 3 章 スカイラジオメーターのデータ処理
12
先ほどはエアロゾル光学特性データの先頭 39 行を読み飛ばしたが、ここにはスカイラジオメーター
の各チャンネルに対応する波長が書かれている。この情報は後で MODIS 画像の大気補正をするために
必要なので、これを読み出してみよう。
エアロゾル光学特性データのヘッダー (一部省略)
(前略)
(i)
(ii)
Data analysis software (SKYRAD.pack)
Cloud screening software (CSSR)
: version4.2
: version1.0
Gobal irradiance data for cloud screening
: Applied · · ·
Calibration constant (F0) for direct radiation : Improved · · ·
-9.99999
: Data not available
(iii)
(iv)
(v)
Channel no.
----------1
Wavelength(um)
-------------0.315
Calibration constant(F0)
------------------------9.99999
2
3
0.340
0.380
5.303E-06 (I)
2.024E-05 (I)
4
5
6
0.400
0.500
0.675
9.028E-05 (I)
2.413E-04 (I)
2.709E-04 (I)
7
8
0.870
0.940
2.235E-04 (I)
-9.99999
9
10
11
1.020
1.627
2.200
1.762E-04 (I)
-9.99999
-9.99999
(後略)
波長データの読み込み (IPython 入力 or skyrad_read_wlen.py)
import re
with open("12051900opticalL2R01.TXT","r") as fp:
lines = fp.readlines()
indx = None
for i,line in enumerate(lines): # (i,line) = (0,lines[0]),(1,lines[1]),...
if re.search("Channel no.",line):
indx = i+2
break
wlen = []
for line in lines[indx:indx+NCH]:
items = line.split()
wlen.append(float(items[1]))
wlen = array(wlen) # list -> array
まず、ファイル全体を読み込んで各行を lines というリストに保存する。次に、re という正規表現
(Regular Expression) モジュールを使って Channel no. と書かれた行を探し出す。その 2 行後より (0
列目から数えて)1 列目の要素をチャンネル個数分だけ取り出して wlen というリストに保存する。最後
に、リストを array に変換している。これは array の方が全要素に対する演算や要素の選択等が簡単に
出来て使いやすいからである。
3.1. エアロゾル光学特性データの読み込み
13
エアロゾル光学的厚さの波長依存性プロット (IPython 入力 or skyrad_aot_vs_wlen.py)
figure() # open new canvas
grid() # draw grids
minorticks_on() # draw small tick marks
cnd = aot[0]>0.0 # cnd = array([False, True, ...], dtype=bool)
aot550 = interp(0.55,wlen[cnd],aot[0][cnd])
plot(wlen[cnd],aot[0][cnd],"o-")
plot([0.55],[aot550],"ro")
xlabel("Wavelength ($\mu$m)") # TeX-like special symbol
ylabel("Aerosol Optical Thickness")
これで図 3.3 のような図が得られる。ここでは後の解析で必要になる波長 550 nm の光学的厚さ (aot550)
を直線補間 (interp) によって求めている。光学的厚さの値はしばしば波長 550 nm の値で規格化される
ことがある。
0.30
Aerosol Optical Thickness
0.25
0.20
0.15
0.10
0.050.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
Wavelength (µm)
図 3.3: エアロゾル光学的厚さの波長依存性
1.1
第 3 章 スカイラジオメーターのデータ処理
14
3.2
エアロゾル物理特性データの読み込み
エアロゾル物理特性データは、エアロゾル光学特性データとほぼ同じ形式になっている。ヘッダーに
は波長とともに体積分布に対応する粒半径の値 (20 個) も含まれている。
エアロゾル物理特性データのヘッダー (粒半径部分)
Radius:
1.209E-06
|
1.768E-06
2.586E-06
3.782E-06
{z
粒半径 20 個
5.530E-06
···
}
この中にある粒半径の値は以下のように読み出すことが出来る。
粒半径データの読み込み (IPython 入力 or skyrad_read_radius.py)
with open("12051900physicalL2R01.TXT","r") as fp:
for line in fp:
if re.search("Radius",line):
radius = float64(line.split()[1:])
break
radius *= 1.0e4 # cm -> um
演習問題 2 スカイラジオメーターの物理特性データ (12051900physicalL2R01.TXT) には
波長 (Wavelength) ごとの複素屈折率実部 (mr)、虚部 (mi)、および粒半径 (Radius) ごとの体積分
布 (Volume) のデータが入っている。このデータを読んで、複素屈折率実部 vs. 波長、複素屈折率
虚部 vs. 波長、体積分布 vs. 粒半径のグラフを作ってみよう。
ヒント: ヘッダー部分は光学特性データと同じく 39 行であり、この中に複素屈折率に対応する波
長および体積分布に対応する粒半径が書かれている。データ部分は 49 列あるが、0 列目を除く 1
∼48 列目までを data に読み込む。mr は data の 6∼16 列目、mi は 17∼27 列目、Volume は 28
∼47 列目である (表 3.2 参照)。
表 3.2: スカイラジオメーターのエアロゾル物理特性データ (Level を除く)
0
Year
1
Month
2
Day
3
Hour(LT)
4
Jday(UT)
5
SZA
6
mr(Ch. 1)
7
mr(Ch. 2)
8
mr(Ch. 3)
9
mr(Ch. 4)
10
mr(Ch. 5)
11
mr(Ch. 6)
12
mr(Ch. 7)
13
mr(Ch. 8)
14
mr(Ch. 9)
15
mr(Ch.10)
16
mr(Ch.11)
17
mi(Ch. 1)
18
mi(Ch. 2)
19
mi(Ch. 3)
20
mi(Ch. 4)
21
mi(Ch. 5)
22
mi(Ch. 6)
23
mi(Ch. 7)
24
mi(Ch. 8)
25
mi(Ch. 9)
26
mi(Ch.10)
27
mi(Ch.11)
28
Volume( 0)
29
Volume( 1)
30
Volume ( 2)
31
Volume( 3)
32
Volume( 4)
33
Volume( 5)
34
Volume( 6)
35
Volume( 7)
36
Volume( 8)
37
Volume( 9)
38
Volume(10)
39
Volume(11)
40
Volume(12)
41
Volume(13)
42
Volume(14)
43
Volume(15)
44
Volume(16)
45
Volume(17)
46
Volume(18)
47
Volume(19)
15
第 4 章 MODIS 衛星画像のデータ処理
4.1
MODIS とは
MODIS (Moderate Resolution Imaging Spectroradiometer) は Terra (EOS AM)、Aqua (EOS PM)
衛星に搭載されたイメージセンサである (図 4.1、表 4.1、表 4.2 参照)。Terra は現地時間の午前に赤道を
北から南に、Aqua は午後に南から北に通過する。Terra/Aqua MODIS は 1∼2 日ごとに全球をカバー
し、36 波長帯のデータを取得する。MODIS のデータは全球の陸域、海域、低層大気における力学や様々
な過程の理解に役立てられている。
図 4.1: MODIS センサー (http://mcst.gsfc.nasa.gov/calibration/information)
表 4.1: MODIS 主要諸元
Orbit
705 km, 10:30 a.m. descending node (Terra) or 1:30 p.m. ascending node (Aqua), sun-synchronous, near-polar, circular
Scan Rate
Swath Dimensions
Telescope
20.3 rpm, cross track
2330 km (cross track) by 10 km (along track at nadir)
17.78 cm diam. off-axis, afocal (collimated), with intermediate
Size
Weight
Power
Data Rate
Quantization
Spatial Resolution
Design Life
field stop
1.0 x 1.6 x 1.0 m
228.7 kg
162.5 W (single orbit average)
10.6 Mbps (peak daytime); 6.1 Mbps (orbital average)
12 bits
250 m (bands 1-2)
500 m (bands 3-7)
1000 m (bands 8-36)
6 years
http://modis.gsfc.nasa.gov/about/specifications.php
第 4 章 MODIS 衛星画像のデータ処理
16
表 4.2: MODIS 観測波長帯
Bandwidth
Spectral
Radiance1
Required
Sensitivity2
1
620 - 670 nm
21.8
128 σ
2
841 - 876 nm
24.7
201 σ
3
4
459 - 479 nm
545 - 565 nm
35.3
29.0
243 σ
228 σ
5
6
7
1230 - 1250 nm
1628 - 1652 nm
2105 - 2155 nm
5.4
7.3
1.0
74 σ
275 σ
110 σ
8
9
10
405 - 420 nm
438 - 448 nm
483 - 493 nm
44.9
41.9
32.1
11
12
526 - 536 nm
546 - 556 nm
27.9
21.0
880 σ
838 σ
802 σ
754 σ
750 σ
13
14
15
662 - 672 nm
673 - 683 nm
743 - 753 nm
9.5
8.7
10.2
910 σ
1087 σ
586 σ
16
17
862 - 877 nm
890 - 920 nm
6.2
10.0
516 σ
167 σ
18
19
20
931 - 941 nm
915 - 965 nm
3.660 - 3.840 µm
3.6
15.0
0.45 (300K)
57 σ
250 σ
0.05 K
Surface/Cloud
Temperature
21
22
3.929 - 3.989 µm
3.929 - 3.989 µm
2.38 (335K)
0.67 (300K)
2.00 K
0.07 K
Atmospheric
Temperature
23
24
25
4.020 - 4.08 µm0
4.433 - 4.498 µm
4.482 - 4.549 µm
0.79 (300K)
0.17 (250K)
0.59 (275K)
0.07 K
0.25 K
0.25 K
Cirrus Clouds
Water Vapor
26
27
1.360 - 1.390 µm
6.535 - 6.895 µm
6.00
1.16 (240K)
150 σ
0.25 K
28
29
30
7.175 - 7.475 µm
8.400 - 8.700 µm
9.580 - 9.880 µm
2.18 (250K)
9.58(300K)
3.69 (250K)
0.25 K
0.05 K
0.25 K
31
32
10.780 - 11.280 µm
11.770 - 12.270 µm
9.55 (300K)
8.94 (300K)
0.05 K
0.05 K
33
34
35
13.185 - 13.485 µm
13.485 - 13.785 µm
13.785 - 14.085 µm
4.52 (260K)
3.76 (250K)
3.11 (240K)
0.25 K
0.25 K
0.25 K
36
14.085 - 14.385 µm
2.08 (220K)
0.35 K
Primary Use
Land/Cloud/Aerosols
Boundaries
Land/Cloud/Aerosols
Properties
Ocean Color/
Phytoplankton/
Biogeochemistry
Atmospheric
Water Vapor
Cloud Properties
Ozone
Surface/Cloud
Temperature
Cloud Top
Altitude
Band
1· · · Spectral Radiance in W/m2 /µm/sr
2· · · Signal-to-noise ratio in σ or Noise-equivalent temperature difference in K
Spatial
Resolution
250 m
500 m
1 km
4.1. MODIS とは
17
MODIS の各観測バンドの応答関数を図 4.2 に示す。この図の作成に使ったデータは
~/work/2012VL_PartB/doc/MODIS/MODIS_Band_Pass に置いてある。
Relative Spectral Response
8
9
3
10
11 4
12
1.0
0.5
0.0
0.40
0.45
0.50
0.55
0.60
15
2 16
17
1.0
0.5
0.0
0.70
0.75
0.80
0.85
0.90
26
6
1.0 5
0.5
0.0
1.4
1.2
1.6
1.8
20
21 23
24
1.0
22
0.5
0.0
4.4
3.6
3.8
4.0
4.2
29
30
31
32
1.0 27 28
0.5
0.0
7
8
9
11
10
Wavelength (µm)
12
1
13
0.65
18
19
0.95
14
0.70
1.00
7
2.0
25
2.2
4.6
4.8
33 34 36
35
13
14
図 4.2: MODIS 観測波長帯
MODIS のデータ取得方法を知っておくと L1B データの構造が理解しやすい。MODIS を搭載している
Terr/Aqua 衛星が Along-Track 方向に進む間、MODIS の回転鏡を掃引して Cross-Track 方向にスキャ
ンが行われる。1 スキャン当り 1.4771 秒の時間を要し、1354 フレームの観測が行われる。1km、500m、
250mMCST
の観測バンドでそれぞれ
10MODIS
、20、40
個の検出器が
L1B
Page 9
Level
1B ProductAlong-Track
User’s Guide方向に並べられており、各検出器
Software Support Team 1、MOD_PR02
February 27, 2009 L1B
は 1 フレーム当りそれぞれ
2、4 個のサンプルを取得する。例えば約
5 分間のデータをまとめた
V6.1.0 (TERRA)/V6.1.1 (AQUA)
データ (Granule) には 203 または 204 スキャンが含まれている (図 4.3 参照)。
Cross Track (Along Scan)
1 measurement per frame
Resolution
1km
500m
250m
Measurements
1354
2708
5416
Along Track (203 scans)
1 measurement per scan
Resolution
1km
500m
250m
Measurements
2030
4060
8120
図 4.3: MODIS データ取得のイメージ (MODIS Level 1B Product User ’s Guide より)
Figure 2-3: Schematic of a MODIS Observation Swath.
This illustrates that scan lines are approximately perpendicular to the direction of motion of the
satellite ground track. (obtained from [4]).
Section 2.2.1
Overview of the HDF-EOS Swath Structure
第 4 章 MODIS 衛星画像のデータ処理
18
4.2
4.2.1
MODIS L1B データ処理
HDF ファイルの閲覧
MODIS L1B データは HDF と呼ばれるデータフォーマットになっている。HDF とは Hierarchical
Data Format(階層的データフォーマット) のことであり、現在は主にバージョン 4 と 5 が使われている。
衛星データには HDF や NetCDF(Network Common Data Form) といったデータフォーマットが良く
使われる。HDF ファイルを閲覧するには、HDFView、Panoply といった GUI ベースのツールが便利で
ある。(Panoply は NetCDF ファイルの閲覧も出来る。) 他にも hdp や ncdump といった CUI ベースの
ツールがある。本来 ncdump は NetCDF ファイルの閲覧ツールであるが、NetCDF のバージョン 4 以降
では HDF ファイルにも対応している。ここでは HDFView を使って L1B データの中身を見てみよう。
HDFView の起動
~/work/2012VL_PartB/data> hdfview MOD021KM.J20120519011621.20120519012656.hdf &
図 4.4: HDFView を使った MODIS L1B データの閲覧
4.2. MODIS L1B データ処理
19
HDFView で MODIS L1B データを開くと、左のパネルにデータの階層構造が表示される。ファイル
名の横にある「4」と書かれたアイコンは HDF バージョンが 4 であることを示している。その下にあ
る「MODIS SWATH Type L1B」のフォルダをダブルクリックすると、「Geolocation Fields」、「Data
Fields」、
「Swath Attributes」という 3 つのフォルダが表示される。
「Data Fiels」をダブルクリックする
と、たくさんのデータセット名が表示される。一番上にある「EV 1KM RefSB」を左クリックすると、
下のパネルにいろいろな情報 (attribute) が表示される。右クリックして表示されるメニューから「Open
As」を選択すると「Dataset Selection - · · · 」という画面 (図 4.4) が表示されるので、「Image」を選択
し、
「Height」、
「Width」、
「Depth」を設定して「Ok」ボタンを押す。雲のようなモノクロ画像が表示さ
れるが、メニューバーにある矢印ボタンを使って「Depth」を 0 (MODIS band 8) から 4 (MODIS band
12) に変えると色が明るくなる。縦のスライドバーを真ん中辺りに持ってくると日本の本州らしい地形
が見えるはずである。メニューバーにあるズームボタンとスライドバーを使って千葉周辺を拡大表示し
てみよう。ここで「 Image」メニューの中から「Set Value Range」を選ぶと図 4.6 左側のような画面が
表示されるので、適当に調節する。1 つの観測波長帯のデータは本来モノクロ画像であるが、メニュー
バーにある「Palette」ボタンもしくは「 Image」メニューの「Change Palette」を選択して表示される
画面で図 4.6 右側のようにカラーパレットを「Rainbow」にすると、図 4.5 右側のように値に応じて色
付けすることもできる。
図 4.5: HDFView による MODIS L1B データの画像化
図 4.6: HDFView のレンジ調整 (左)、カラーパレットの選択 (右)
第 4 章 MODIS 衛星画像のデータ処理
20
4.2.2
HDF ファイルの読み込み
前述のように、MODIS L1B のデータフォーマットは HDF4 形式になっており、HDFView というツー
ルを使うとデータの閲覧や簡単な描画が出来る。しかしながら、HDFView では複雑な解析は出来ない
ため、Python で読み込んでみよう。これは PyHDF というモジュールを使って簡単に出来る。いかに簡
単かを実感してもらうため、ここでは詳細にはこだわらずに最短のコードを示す。前提として、
~/work/2012VL_PartB/data> ipython --pylab
のように data ディレクトリで IPython を立ち上げているものとする。(2012VL_PartB ディレクトリで
IPython を立ち上げている場合は、IPython 上で
In : cd data
のようにしてディレクトリを移動する。)
MOD02 1km 分解能データの読み込み+描画 (IPython 入力 or modis_read_1km.py)
from pyhdf.SD import SD,SDC
f = SD("MOD021KM.J20120519011621.20120519012656.hdf",SDC.READ)
v = f.select("EV_1KM_RefSB")
imshow(v[0]) # 0 => MODIS band 8
このように、たった 4 行で図 4.7 のような図が得られる。ちなみに、同じ図を得るだけであれば、2 行
目から 4 行目を以下のような 1 行のコマンドにまとめることも出来る。
In : imshow(SD("HDF ファイル名",SDC.READ).select("EV_1KM_RefSB")[0])
0
1000
2000
3000
4000
0200400600800
1000
1200
図 4.7: MOD02 1km 分解能データ (バンド 8、ピクセル座標)
最初なので、ここで上記コードを 1 行ずつ詳しく説明する。
• 1 行目 (from pyhdf.SD import SD,SDC)
pyhdf.SD モジュールから SD クラス、SDC クラスを import する。(SD は Scientific Dataset の略。
SDC の C は Constants の略。) これは IPython を立ち上げてから SD、SDC の使用前に 1 度だけ行え
ば良いが、2 度以上行っても害はない。(ただし、IPython を立ち上げた後に import したモジュー
ルのソースを書き換えたとしても、更新は反映されない点に注意。)
• 2 行目 (f = SD("HDF ファイル名",SDC.READ))
"HDF ファイル名"を読み込み、そのハンドル (SD インスタンス) を f という変数に格納する。
4.2. MODIS L1B データ処理
21
【補足事項 】
In : f.datasets()
とタイプすると、f からアクセス可能なデータセットの辞書型配列 (表 4.3 参照) が得られる。
In : f.datasets().keys()
でデータセットの名前だけを取得することも出来る。このままだと順番がバラバラで少し見にく
いので、
In : sorted(f.datasets().keys())
のようにした方が分かりやすい。これらデータセットのサイズを調べたければ、
In : for k in sorted(f.datasets().keys()): print k,f.select(k)[:].shape
とすれば良いだろう。Latitude、Longitude 等のデータセットは along-track、along-scan 方向で
それぞれ 1/5 に間引かれている。
表 4.3: MOD02 の主要データセット
データセット名
サイズ (例)
Band_1KM_RefSB
(15,)
EV_1KM_RefSB
Latitude
(15, 4300, 1354)
(860, 271)
Longitude
SensorAzimuth
SensorZenith
(860, 271)
(860, 271)
(860, 271)
SolarAzimuth
SolarZenith
(860, 271)
(860, 271)
In : sorted(f.attributes().keys())
で f が持つ属性のリストが得られる。(f.attributes メソッドは辞書型配列を返す。) 例えば、
In : f.attributes()["Earth-Sun Distance"]
を見ると地球–太陽間距離 (A.U.) が得られる。
• 3 行目 (v = f.select("EV_1KM_RefSB"))
"EV_1KM_RefSB"というデータセットのハンドル (SDS インスタンス) を v という変数に格納する。
【補足事項 】
In : sorted(v.attributes().keys())
で v が持つ属性のリスト (表 4.4 参照) が得られる。(v.attributes メソッドは辞書型配列を返
す。) この中にはバンド名や校正係数等、後の解析で必要な情報が含まれている。band_names は
カンマ区切りの文字列になっているが、
In : v.attribute()["band_names"].split(",")
とすればバンド名のリストが得られる。
例えば IPython に
In : v[:]
と入力して中身を表示させて見ると分かるが、変数 v のデータは整数値 (Scaled Integer: SI) に
なっている。SI は以下の式を使って物理量 (放射輝度 L または地表面反射率 ρ cos θ、θ は太陽天頂
角) に変換できる。
第 4 章 MODIS 衛星画像のデータ処理
22
L
ρ cos θ
= radiance_scales × (SI − radiance_offsets)
(4.1)
= reflectance_scales × (SI − reflectance_offsets)
(4.2)
0
SI
H
YH
I
@
HH
@ scaling
H
HH @
H @
Lmin L
32767
*
Lmax
図 4.8: 放射輝度のスケーリング
例えば、放射輝度の測定値 L は、放射輝度の最小値 Lmin 、最大値 Lmax を使って図 4.8 のように
スケーリングされる。ここで、radiance_scales、radiance_offsets は Lmin 、Lmax と以下の
ような関係がある。
radiance_scales = (Lmax − Lmin )/32767
radiance_offsets = −(32767 × Lmin )/(Lmax − Lmin )
(4.3)
(4.4)
式 4.1 を使うよりも、式 4.2 と f の属性にある"Solar Irradiance on RSB Detectors over pi"、
"Earth-Sun Distance" を使うとより正確な放射輝度を求めることができる。
L
=
(Solar Irradiance on RSB Detectors over pi)
× ρ cos θ
π(Earth-Sun Distance)2
(4.5)
表 4.4: EV_1KM_RefSB データセットの主要属性
属性
値 (例)
_FillValue
band_names
65535
’8,9,10,11,12,13lo,13hi,14lo,14hi,15,16,17,18,19,26’
radiance_offsets
radiance_scales
[316.9721984863281,316.9721984863281,...]
[0.01298768911510706,0.008816360495984554,...]
radiance_units
reflectance_offsets
reflectance_scales
’Watts/m^2/micrometer/steradian’
[316.9721984863281,316.9721984863281,...]
[2.3924452762003057e-05,1.489247824792983e-05,...]
reflectance_units
valid_range
’none’
[0,32767]
• 4 行目 (imshow(v[0]))
変数 v のデータ (EV_1KM_RefSB) から 2 次元プロットを作成する。変数 v のデータはもともと
[band, along-track, along-scan] の 3 次元 array になっている。このうち、band = 0 (MODIS
band 8) のデータを使用している。imshow は 2 次元 array データを直交座標系の等間隔グリッド
上に割り当てて画像化する Matplotlib ライブラリの関数である。IPython を起動する時に--pylab
のオプションを付けたためにこのように直接呼び出すことが出来る。スクリプトで同様の呼び出
しをするには、
from matplotlib.pylab import *
という 1 文が必要である。
4.2. MODIS L1B データ処理
23
これまでは詳細にこだわらず、とにかくデータを表示してみたが、以下の点で改善が必要である。
1. 整数値 (SI) を物理量に変換する。
2. 表示する値のレンジを限定する。
3. 注目している領域 (千葉周辺) を拡大する。
ここでは整数値をレンジが分かりやすい地表面反射率 (reflectance) に変換する。波長帯は、データが比
較的綺麗なバンド 12 (波長 546∼556 nm、0 から数えて 4 番目のバンド) を選ぶ。
MOD02 1km 分解能データの描画の改善 (IPython 入力 or modis_plot_1km.py)
v_attr = v.attributes()
reflectance_offsets = v_attr["reflectance_offsets"]
reflectance_scales = v_attr["reflectance_scales"]
band = 4 # MODIS band 12
reflectance = reflectance_scales[band]*(v[band]-reflectance_offsets[band])
figure() # open new canvas
imshow(reflectance,vmin=0.0,vmax=0.2) # reflectance 0 -- 0.2
colorbar()
x1,x2 = 372,513
y1,y2 = 1829,1983
xlim(x1,x2)
ylim(y2,y1) # (y1,y2) -> flip vertical
これらのコマンドを追加すると、以下のような図が得られる。大分マシになった。しかしながら、まだ
不満が残る。縦軸、横軸の値がピクセル番号になっていることだ。これについては、次節で改善するこ
とにしよう。
0.20
1840
0.18
0.16
1860
0.14
1880
0.12
1900
0.10
1920
0.08
1940
0.06
0.04
1960
1980
0.02
380
400
420
440
460
480
500
0.00
図 4.9: MOD02 1km 分解能反射率 (バンド 12、ピクセル座標)
第 4 章 MODIS 衛星画像のデータ処理
24
4.2.3
地理位置情報の読み込み
まず、先ほどの図の縦軸、横軸がそれぞれ緯度、経度になるように修正してみよう。そのためには各
ピクセルの緯度、経度の値が必要である。1km 分解能の MOD02 データにも緯度、経度の情報が含まれ
ているが、これらは前述のように 1/5 に間引かれている。実はこのデータに対応する Geolocation デー
タ (MOD03) というものがあり、1km 分解能の緯度、経度の値はこのデータに含まれている。その他、
太陽方向、センサー方向等、MODIS 衛星画像の大気補正に必要な情報もこのデータに含まれている。
MOD03 緯度、経度データの読み込み+描画 (IPython 入力 or modis_geolocation.py)
g = SD("MOD03.J20120519011621.20120519012656.hdf",SDC.READ)
lat = g.select("Latitude")[y1:y2,x1:x2]
lon = g.select("Longitude")[y1:y2,x1:x2]
dat = reflectance[y1:y2,x1:x2]
figure() # open new canvas
pcolor(lon,lat,dat,vmin=0.0,vmax=0.2)
colorbar()
これで下のような図が得られる。見てすぐ分かることは、描画された領域が軸に対して斜めになってい
ることである。これは、along-track、along-scan 方向が緯度、経度方向と異なることが原因である。こ
のまま解析することもできるが、使いやすいように地図投影法を変換しておこう。
その前に、上のコードでいくつか注意する点を挙げておく。まず、
In : lat = g.select("Latitude")[y1:y2,x1:x2]
の部分であるが、ここでは緯度データのうち、[y1:y2,x1:x2] の部分を切り出している。(これをスラ
イスと言う。) このように、2 次元 array のインデックスは [y,x] の順序で書く。(後方のパラメータの
ループが先に回る。) これは C 言語と同じである。次に、
In : pcolor(lon,lat,dat,vmin=0.0,vmax=0.2)
の部分であるが、pcolor メソッドの引数は (x,y,z) の順番になっている。(x、y、z はそれぞれ横軸、
縦軸、カラーバーの変数である。) このように、大抵のメソッドの引数は (x,y,z) の順番になっている
ので、インデックスの順番と混同しないように注意しよう。
36.2
0.20
36.0
0.18
0.16
35.8
0.14
35.6
0.12
35.4
0.10
35.2
0.08
0.06
35.0
0.04
34.8
34.6
139.0
0.02
139.5
140.0
140.5
141.0
141.5
0.00
図 4.10: MOD02 1km 分解能反射率 (バンド 12、緯度経度座標)
4.2. MODIS L1B データ処理
4.2.4
25
地図投影法の変換
2 次元補間を利用すれば、2 次元 array のグリッドを様々な地図投影法のグリッドに変換することが出
来る。このような地図投影法の変換は MODIS 用に開発された MODIS Reprojection Tool (MRTSwath)
というツールを使っても可能であるが、ここでは Python を使ってやってみよう。(MRTSwath は Land
Processes Distributed Active Archive Center (LP DAAC) に保存されているデータ (MOD11 L2、
MYD11 L2、MOD14、MYD14) のためにデザインされており、L1B データに適用した結果は保証さ
れていない。)
2 次元補間 (IPython 入力 or modis_reprojection.py)
lat_min,lat_max = 34.8,36.0
lon_min,lon_max = 139.5,141.0
pixel_size = 0.009 # 0.009 deg ~ 1 km
lat_map = arange(lat_min,lat_max,pixel_size)
lon_map = arange(lon_min,lon_max,pixel_size)
dat_map = griddata(lon.flatten(),lat.flatten(),dat.flatten(),lon_map,lat_map)
figure() # open new canvas
pcolor(lon_map,lat_map,dat_map,vmin=0.0,vmax=0.2)
colorbar()
xlim(lon_min,lon_max)
ylim(lat_min,lat_max)
これで図 4.11 のような図が得られる。2 次元補間には griddata という Matplotlib のメソッドを使って
いる。本演習で使われる補間アルゴリズムは NCAR natgrid library というライブラリに実装された強
固なものである。(付録 A.7.2 参照)
36.0
0.20
0.18
35.8
35.6
0.16
0.14
0.12
35.4
0.10
0.08
35.2
35.0
0.06
0.04
0.02
34.8 139.6 139.8 140.0 140.2 140.4 140.6 140.8 141.0
0.00
図 4.11: MOD02 1km 分解能反射率 (バンド 12、2 次元補間データ)
なお、図 4.11 の領域をぎりぎり内包するピクセル番号は以下のように調べられる。
In : indx = np.where((lat>=lat_min)&(lat<=lat_max)&(lon>=lon_min)&(lon<=lon_max))
In : x1 = max(int(indx[1].min())-1,0)
In : x2 = min(int(indx[1].max())+2,lat.shape[1])
In : y1 = max(int(indx[0].min())-1,0)
In : y2 = min(int(indx[0].max())+2,lat.shape[0])
第 4 章 MODIS 衛星画像のデータ処理
26
【補足事項 1 】
本演習の目的には上記の変換で十分であるが、Matplotlib には地図を扱う上で便利な Basemap という
モジュールが提供されているので、ここでその使い方を見てみよう。
Basemap モジュールの使用例 (IPython 入力 or modis_basemap.py)
from mpl_toolkits.basemap import Basemap
m = Basemap(projection="cyl",
lat_0=lat.mean(),lon_0=lon.mean(),
llcrnrlat=lat.min(),llcrnrlon=lon.min(), # lower left corner
urcrnrlat=lat.max(),urcrnrlon=lon.max(), # upper right corner
resolution="h", # c(crude),l(low),i(intermediate),h(high),f(full)
area_thresh=1000.0, # threshold in km^2
fix_aspect=False)
figure() # open new canvas
m.pcolor(lon,lat,dat,vmin=0.0,vmax=0.2)
m.drawcoastlines()
m.drawparallels(arange(30.0,41.0,0.5),labels=[1,0,0,0])
m.drawmeridians(arange(130.0,150.0,0.5),labels=[0,0,0,1])
colorbar()
# labels ... left,
# right,top,bottom
0.20
0.18
36°N
0.16
0.14
0.12
35.5°N
0.10
0.08
0.06
35°N
0.04
0.02
139.5°E
140°E
140.5°E
141°E
0.00
図 4.12: MOD02 1km 分解能反射率 (バンド 12、正距円筒図法)
上の例は Equidistant Cylindrical Projection の例であるが、その他にも以下のように様々な地図投影
法に対応している。
Azimuthal Equidistant, Gnomonic, Orthographic, Geostationary, Near-Sided Perspective, Mollweide,
Hammer, Robinson, Eckert IV, Kavrayskiy VII, McBryde-Thomas Flat Polar Quartic Sinusoidal,
Equidistant Cylindrical, Cassini, Mercator, Transverse Mercator, Oblique Mercator, Polyconic, Miller
Cylindrical, Gall Stereographic, Lambert Conformal, Lambert Azimuthal Equal Area, Stereographic,
Equidistant Conic, Albers Equal Area, Polar Stereographic, Polar Lambert Azimuthal, Polar Azimuthal Equidistant, van der Grinten
centers of the first 1000, 500, and 250 meter pixels across track on each earth scan are coregistered, and that each earth scan must be interpolated separately because of the overlap
between successive earth scans.
General purpose interpolation routines such as CONGRID and REBIN in IDL should not
be used. Figure 1 shows how the MODIS 1000 meter and 500 meter pixel locations are
co-registered at the start of an earth scan. The first three of the ten MODIS 1000 meter
4.2. MODIS L1B データ処理
27
detectors along-track are depicted. The blue outlines represent the 1000 meter pixels,
while the green outlines represent the 500 meter pixels. The solid blue squares are the
centers
【補足事項
2 of
】the 1000 meter pixels, and the slid green squares are the centers of the 500
pixels. Note
the centers of the
leftmost
pixels 500m
are co-registered.
Also note how
MODIS L1Bmeter
データには
1km how
グリッドの他に、
250m
グリッド、
グリッドのデータがあるが、
Gethe topmost row of 500 meter pixel centers lie outside the bounds of the 1000 meter pixel
olocation データには
centers. 1km グリッドの緯度経度しか含まれない。したがって、250m、500m グリッドの
値は補間によって求める必要がある。1km グリッドと 500m グリッドの配置は図 4.13 のようになって
Figure 1: MODIS 1000 and 500 meter pixel registration
いる。
1000 meter pixels
500 meter pixels
図 4.13:
MODIS the
1000m
500m グリッド
True Color
Tutorial より)
To interpolate
1000グリッドと
meter geolocation
data to 500(Modis
meter resolution,
a bilinear
interpolation algorithm is applied to each earth scan. For the first and last row of 500
meter pixels, bilinear extrapolation is used. A similar method is used to interpolate the
geolocation
data to 250 metermeshgrid
resolution.という Matplotlib のメソッドで作成出来る。緯度経度
このような等間隔
2 次元グリッドは
の補間には、地図投影法の変換で使った
が利用出来る。以下のコードで図
4.14 のような図が
IDL routines: get_geodata.pro,griddata
modis_geo_interp_500.pro,
modis_geo_interp_250.pro
得られる。
MOD02 500m 分解能データの読み込み+描画 (IPython 入力 or modis_plot_hkm.py)
11
xk,yk = meshgrid(arange(x1,x2),arange(y1,y2)) # Make 1km grid
xh,yh = meshgrid(arange(x1,x2,0.5),arange(y1,y2,0.5)-0.25) # Make 500m grid
lat_hkm = griddata(xk.flatten(),yk.flatten(),lat.flatten(),xh,yh)
lon_hkm = griddata(xk.flatten(),yk.flatten(),lon.flatten(),xh,yh)
f = SD("MOD02HKM.J20120519011621.20120519012656.hdf",SDC.READ)
v = f.select("EV_500_RefSB")
v_attr = v.attributes()
reflectance_offsets = v_attr["reflectance_offsets"]
reflectance_scales = v_attr["reflectance_scales"]
band = 1 # MODIS band 4
xh_min,xh_max = x1*2,x2*2
yh_min,yh_max = y1*2,y2*2
v_data = v[band,yh_min:yh_max,xh_min:xh_max]
dat = reflectance_scales[band]*(v_data-reflectance_offsets[band])
figure() # open new canvas
pcolor(lon_hkm,lat_hkm,dat,vmin=0.0,vmax=0.2)
colorbar()
28
第 4 章 MODIS 衛星画像のデータ処理
図 4.14: MOD02 500m 分解能反射率 (バンド 4)
以上は 500m グリッドの場合であるが、250m グリッドもほぼ同様である (modis_plot_qkm.py 参照)。
図 4.15: MOD02 250m 分解能反射率 (バンド 1)
4.2. MODIS L1B データ処理
29
【補足事項 3 】
7. Data product grids
MODIS Level 2G、Level 3、Level 4 プロダクトでは Sinusoidal Projection という投影法が用いられて
7.1. MODIS sinusoidal grid
おり、全球が下図のように 1 区画のサイズが約 10◦ × 10◦ のタイル状に区切られている。
The Earth's surface is divided into tiles (10° x 10°).
Figure
1. The MODIS
sinusoidal
grid with an
example tile
shown asReflectance
a RGB-image produced
the より)
図 4.16:
MODIS
Sinusoidal
Projection
(MODIS
Surface
User ’sfrom
Guide
MODIS data acquired on December 3, 2006 over the US East cost. This product Granule ID is
MOD09A1.A2000337.h11v05.005.2006342055602.hdf.
v = V 、h = H のタイル内にあるグリッド (i, j) の中心緯度経度 (ϕ,λ) は以下の様に求められる。
w
= T /N
(4.6)
x
= (j + 0.5)w + HT + xmin
(4.7)
y
= ymax − (i + 0.5)w − V T
y
=
R
x
=
R cos ϕ
(4.8)
ϕ
λ
(4.9)
(4.10)
ここで、N はタイルの分割数 (i, j = 0, 1, ..., N )、T はタイルサイズ (1111950 m)、xmin 、ymax はマッ
プの西端 (-20015109 m) および北端 (10007555 m)、R は地球半径 (6371007.181 m) である。
7
第 4 章 MODIS 衛星画像のデータ処理
30
MOD09 500m 分解能データの読み込み+描画 (IPython 入力 or modis_sinusoidal.py)
H,V = 29.0,5.0
R,T = 6371007.181,1111950.0
xmin,ymax = -20015109.0,10007555.0
N = 2400
w = T/N
i = arange(N).astype("float")
j = arange(N).astype("float")
x = (j+0.5)*w+H*T+xmin
y = ymax-(i+0.5)*w-V*T
lat = hstack([y.reshape(N,1)/R]*N)
lon = x/(R*cos(lat))
f = SD("MOD09GA.A2012140.h29v05.005.2012142054114.hdf",SDC.READ)
v = f.select("sur_refl_b04_1") # Band 4
v_attr = v.attributes()
x1,x2 = 686,1388
y1,y2 = 959,1249
lat_map = degrees(lat)[y1:y2,x1:x2]
lon_map = degrees(lon)[y1:y2,x1:x2]
dat_map = 1.0/v_attr["scale_factor"]*(v[y1:y2,x1:x2]-v_attr["add_offset"])
figure() # open new canvas
pcolor(lon_map,lat_map,dat_map,vmin=0.0,vmax=0.2)
xlim(lon_min,lon_max)
ylim(lat_min,lat_max)
colorbar()
上のコードを実行すると、図 4.17 のような図が得られる。
図 4.17: MOD09GA 500m 分解能反射率 (バンド 4)
31
第 5 章 6S を使った MODIS 衛星画像の大気補正
5.1
大気補正とは
Sensor
Sun
Lobs
Lps Lgd Lpm Lgi L0gi L0pm
Ltar
Lenv Latm
ρ
ρ
Target Environment
図 5.1: 衛星センサーに入射する太陽光
太陽光によって照らされている地表面を衛星から見ると、主に以下の 6 成分からなる放射輝度の総和
(Lobs ) が観測される。
1. Lps : 大気中で単散乱された成分
2. Lpm : 大気中で多重散乱された成分
3. L′pm : 対象周辺 (Environment) で反射された後に大気中で単散乱された成分
1
4. Lgd : 直接的に対象領域 (Target) に到達し、反射された成分
5. Lgi : 大気中で散乱されて対象領域に到達し、反射された成分
6. L′gi : 対象周辺で反射された後に大気中で単散乱されて対象領域に到達し、反射された成分
これらはさらに以下のように分類することができる。
1. Latm : 大気中だけで散乱された成分 (Lps ,Lpm )
2. Lenv : 対象周辺で反射された成分 (L′pm )
3. Ltar : 対象領域で反射された成分 (Lgd 、Lgi 、L′gi )
Lgd 以外の成分は大気科学の研究にとって重要な情報を含んでいるが、地表面の研究にとっては余分な
情報なので、これを取り除く必要があり、そのことを大気補正と呼んでいる。
第 5 章 6S を使った MODIS 衛星画像の大気補正
32
5.2
6S とは
6S とは放射伝達シミュレーションコードの 1 つであり、人工衛星等で測定される物体の放射輝度を計
算したり、放射輝度値に含まれる大気の影響を見積り、それを取り除くこと (大気補正) が出来る。6S と
いう名前は”Second Simulation of a Satellite Signal in the Solar Spectrum” に由来している。元々は
Second Simulation
a Satellite Signal
偏光を考慮しないスカラーコードであったが、
2005 of
年に偏光を考慮したベクターコードがリリースされ
in the Solar
- Vector (6SV) ) からダウンロード出来る。
た。6S 最新版のコードやマニュアルは
WebSpectrum
ページ ( http://6s.ltdri.org/
図 5.2: 6S でシミュレートする地表面観測のイメージ (6S マニュアルより)
Vermote E.(1,*), D. Tanré(2), J. L. Deuzé(2), M. Herman(2),
J. J. Morcrette(3), and S. Y. Kotchenova(1)
6S 内部では前述の Latm 、Lenv 、Ltar は以下のように計算される。
∫
λ
λ
S(λ)EUniversity
θv )ρaλ dλ
(1) Department of1Geography,
ofs ,Maryland
s µs Tg (θ
∫
Latm
=
address for correspondence:
π
S(λ) dλ
4321 Hartwick Road, Suite 209
∫
T λ (θs )tλd (θv )⟨ρλ ⟩
College Park, MD 20740 λ
S(λ)Es µs Tgλ (θs , θv )
dλ
1
(2) Laboratoire d'Optique
Atmosphérique, URA CNRS 1
713
− ⟨ρλ ⟩sλ
∫
LUniversité
env = des Sciences et Technologies de Lille
π
S(λ) dλ
59655 Villeneuve d'Ascq Cédex, France
∫
T λ (θs ) exp(−τ /µv )ρλ
λ Range
λ Weather Forecast
(3) European Centre for
Medium
S(λ)E
dλ
s µs Tg (θs , θv )
Shinfield Park,1Reading
1 − ⟨ρλ ⟩sλ
LBerkshire
= RG2 9AX, United Kingdom ∫
tar
π
S(λ) dλ
(5.1)
(5.2)
(5.3)
ここで、λ:波長、θs :太陽天頂角、
θv :視線方向天頂角、
µs :cos
θs 、µv :cos θv 、S:フィルター関数、Esλ :大気
* Formerly affiliated
with the Laboratoire d'Optique
Atmosphérique
上端太陽放射照度、τ :大気の光学的厚さ、tλd :拡散透過率、T :透過率 (e−τ /µ + td )、Tg :気体分子の吸収に
よる透過率、sλ :球面アルベド、ρλ :地表面反射率、ρaλ :大気の反射率である。記号 ⟨⟩ は空間平均を表す。
6S の大気補正モードを有効にすると、入力された (大気補正されていない) 見かけの反射率 ρ∗i または見
かけの放射輝度 Lobs を再現するように大気補正された反射率 ρac を計算する。(大気補正モードでは、
地表面は一様反射する Lambertian として扱われる。)
ρ∗i
=
ρac
=
ρ′ac
=
πLobs
Es µs
ρ′ac
1 + ρ′ac s
ρ∗i − ρatm
Tg (θs , θv )T (θs )T (θv )
(5.4)
(5.5)
(5.6)
(5.7)
ここで、Es は S(λ)、T 、Tg 、s、ρa は S(λ)Eλ の重み付き平均値、ρatm = ρa Tg である。
5.3. 6S の実行
5.3
33
6S の実行
6S 入力ファイルの主な設定項目を表 5.1 に挙げる。(場合によっては表 5.1 にない項目の設定が必要。
組込みモデル等の詳細は 6S Manual Part 1 参照。)
表 5.1: 6S 入力ファイルの主な設定項目
項目
内容
IGEOM
幾何学条件
IDATM
大気分子モデル
IAER
エアロゾルモデル
V
エアロゾル濃度
−1: エアロゾルなし、 0: 光学的厚さ † を指定、 V>5: 視程 † V km
XPS
ターゲット高度
XPS≥0: 海面高度、 XPS<0: −XPS km
XPP
センサー高度
−1000: 衛星高度、 0: 地上高度、 −100<XPS<0: −XPP km
IWAVE
センサー波長特性
設定値
0: ユーザー指定 (太陽方向、視線方向、観測月日)
1∼7: 組込み衛星
0: 吸収なし、 1∼7: 組込みモデル、 8∼12: ユーザー指定
−1: ユーザー指定 (プロファイル)、 0: エアロゾルなし
1∼7: 組込みモデル、 8∼12: ユーザー指定 (粒径分布)
−2: 波長の上下限を指定。フィルタ関数は一定。出力は 2.5nm
ごと、 −1: 単一波長、 0: 波長の上下限を指定。フィルタ関数は
一定。出力はバンド中心、1: 波長の上下限、2.5nm ごとのフィルタ
関数を指定、2∼165: 組込み衛星
INHOMO
地表面反射モデル
0: 一様、 1: 非一様
IDIREC
反射率角度依存性
0: 角度依存なし、 1: 角度依存あり
IGROUN
地表面反射率
-1: 2.5nm ごとに指定、 0: 一定値、 1∼4: 組込みモデル
IRAPP
大気補正モード
−1: 大気補正なし、 0,1: 大気補正あり
RAPP
大気上端条件
−1 <RAPP<0: ρTOA = −RAPP、
RAPP>0: LTOA = RAPP W/m2 /st/µm
† · · · 光学的厚さ、視程はともに波長 550 nm における値
簡単な設定例を以下に挙げる。
6S 入力ファイルの例 (input.dat)
0
(User-defined geometric conditions)
22.84 128.87 22.90 99.73 5 19 (SZA, SAZ, VZA, VAZ, month, day)
3
(Atmosphere model: Midlatitude Winter)
2
0
0.2
(Aerosol model: Maritime)
(Constant value for tau)
(Aerosol optical thickness at 550 nm tau)
0
-1000
(Target level)
(Satellite level)
45
0
0
(Chosen band: MODIS Band4)
(Homogeneous surface)
(No directional effects)
0
0.1
(Constant value for rho)
(Surface reflectance rho)
0
50.0
(Atm. correction Lambertian)
(Radiance)
第 5 章 6S を使った MODIS 衛星画像の大気補正
34
例えば上の内容を input.dat というファイル名で保存して、以下のように 6S を実行出来る。
6S の実行
~/work/2012VL_PartB/data> sixsV1.1 <input.dat >output.dat
実行結果は標準出力に出力されるが、上のようにファイルにリダイレクトして保存することが出来る。
6S 出力ファイルの例 (output.dat)
******************************* 6SV version 1.1 *******************************
*
*
*
*
geometrical conditions identity
-------------------------------
*
*
*
*
*
user defined conditions
month:
*
*
*
*
*
*
solar zenith angle:
view zenith angle:
scattering angle:
*
*
5 day :
19
22.84 deg
22.90 deg
168.78 deg
solar azimuthal angle:
128.87 deg
view azimuthal angle:
99.73 deg
azimuthal angle difference: 29.14 deg
*
*
*
*
*
atmospheric model description
*
*
*
----------------------------atmospheric model identity :
midlatitude winter (uh2o=.853g/cm2,uo3=.395cm-atm)
*
*
*
*
*
aerosols type identity :
Maritime aerosol model
*
*
*
*
optical condition identity :
visibility : 26.74 km opt. thick. 550 nm :
0.2000
*
*
(中略)
*******************************************************************************
*
atmospheric correction result
*
*
----------------------------*
*
*
*
input apparent reflectance
:
measured radiance [w/m2/sr/mic]
:
atmospherically corrected reflectance
*
*
Lambertian case :
BRDF
case :
0.094
50.000
0.04685
0.04685
*
*
*
*
*
*
coefficients xa xb xc
: 0.00234 0.06997 0.12119
*
*
y=xa*(measured radiance)-xb; acr=y/(1.+xc*y)
*
*******************************************************************************
出力ファイル中の xa、xb、xc、y はそれぞれ π/ (Es µs Tg (θs , θv )T (θs )T (θv ))、ρatm / (Tg (θs , θv )T (θs )T (θv ))、
s、ρ′ac に対応しており、大気補正された地表面反射率が次式によって求められる。
y = xa × Lobs − xb
y
ρac =
1 + xc × y
(5.8)
(5.9)
5.3. 6S の実行
35
今後 6S を実行するのに便利なように、太陽方向、視線方向、観測月日、および光学的厚さを入力し
たら 6S を実行して結果を取り出してくれる Python 関数を用意しよう。
6S 実行関数 (IPython 入力 or sixs_run_simple.py)
import re
from subprocess import Popen,PIPE
def run_6s(sza,saz,vza,vaz,mon,day,tau):
lines = ""
lines += "0\n" # User defined
lines += "%.10f %.10f %.10f %.10f %d %d\n"%(sza,saz,vza,vaz,mon,day)
lines += "3\n" # Midlatitude Winter
lines += "2\n" # Maritime Model
lines += "0\n" # Constant value for tau
lines += "%.10f\n"%(tau) # tau
lines += "0\n" # Target level
lines += "-1000\n" # Satellite level
lines += "42\n" # Chosen band
lines += "0\n" # Homogeneous surface
lines += "0\n" # No directional effects
lines += "0\n" # Constant value for rho
lines += "0.1\n" # rho
lines += "0\n" # Atm. correction Lambertian
lines += "50.0\n" # Radiance (positive value)
p = Popen("sixsV1.1",stdin=PIPE,stdout=PIPE,stderr=PIPE,shell=True)
out,err = p.communicate(lines)
xa = nan
xb = nan
xc = nan
for line in out.splitlines():
m = re.search("coefficients\s+xa\s+xb\s+xc\s*:\s*(\S+)\s+(\S+)\s+(\S+)",
line)
if m:
xa = float(m.group(1))
xb = float(m.group(2))
xc = float(m.group(3))
return xa,xb,xc
まず、6S に入力したい内容を lines という変数に格納している。それから POpen というモジュールを
使って 6S に lines を入力し、標準出力、標準エラー出力をそれぞれ out、err という変数に取り込ん
でいる。out の内容を読み、その中から見つけた xa、xb、xc の値を return するようになっている。
簡単のため、とりあえず大気モデルやエアロゾルモデルには 6S に組み込まれているものを使おう。大
気補正機能を有効にすると 6S に入力する地表面反射率や放射輝度の値は結果に影響しなくなるので、こ
こでは適当な値 (0.1 および 50) を入れておく。
第 5 章 6S を使った MODIS 衛星画像の大気補正
36
下の例では run_6s 関数に与える光学的厚さ τ550 を 0 から 1 まで 0.1 刻みで変化させて大気補正係数
テーブルを作成している。各大気補正係数は τ550 に応じて図 5.3 のように変化する。
大気補正係数テーブル作成スクリプト (IPython 入力 or sixs_calc_coeff.py)
sza,saz = 22.84,128.87
vza,vaz = 22.90,99.73
mon,day = 5,19
tau = arange(0.0,1.01,0.1)
with open("coeff.dat","w") as fp:
for t in tau:
print "Tau: %.3f"%(t)
xa,xb,xc = run_6s(sza,saz,vza,vaz,mon,day,t)
fp.write("%4.1f %8.5f %8.5f %8.5f\n"%(t,xa,xb,xc))
0.16
0.0031
0.14
0.0030
0.0029
xa
0.0028
0.0027
0.0026
0.0
0.2
0.4
0.6
0.8
0.12
0.10
0.08
0.06
0.0025
0.04
1.00.0024
0.02
0.22
0.20
0.18
0.16
0.14
xc
xa
xb
xc
0.0032
xb
0.12
0.10
0.08
0.06
0.04
τ
図 5.3: 大気補正係数 vs. 光学的厚さ (sixs_plot_coeff.py にて作成)
次の例では放射輝度、光学的厚さ、地表面反射率を関係付ける Lookup テーブル (LUT) を作成する。
LUT 作成スクリプト (IPython 入力 or sixs_make_lut.py)
lobs = arange(10.0,100.1,10.0)
tau = arange(0.0,1.01,0.1)
rho = []
for t in tau:
print "Tau: %.3f"%(t)
xa,xb,xc = run_6s(sza,saz,vza,vaz,mon,day,t)
y = xa*lobs-xb
rho.append(y/(1.0+xc*y))
with open("lut.dat","w") as fp:
for i in range(tau.size):
for j in range(lobs.size):
fp.write(" %10.3e"%(rho[i][j]))
fp.write("\n")
5.3. 6S の実行
37
ここで作成した LUT を使って、Lobs vs. ρac および τ550 vs. ρac の関係をプロットしたものをそれぞれ
図 5.4、図 5.5 に示す。
0.0
0.1
100
τ550
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
90
Lobs (W/m2 /sr/µm)
80
70
60
50
40
30
20
10
0.00
0.05
0.10
ρac
0.15
0.20
0.25
図 5.4: Lobs vs. ρac (sixs_lobs_vs_rho.py にて作成)
Lobs (W/m /sr/µm)
2
10.0
20.0
1.0
30.0
40.0
50.0
60.0
70.0
80.0
90.0
100.0
0.8
τ550
0.6
0.4
0.2
0.0
0.00
0.05
0.10
ρac
0.15
0.20
0.25
図 5.5: τ550 vs. ρac (sixs_tau_vs_rho.py にて作成)
第 5 章 6S を使った MODIS 衛星画像の大気補正
38
以前示した 6S の入力ファイルは非常に単純であったが、これを修正すればいろいろなカスタマイズが
可能である。ここではエアロゾルモデルとしてスカイラジオメーターで得られた物理特性 (複素屈折率、
体積分布) を使ってみよう。また、センサー波長特性にはデフォルトにない MODIS バンド 12 のフィル
タ関数を使ってみよう。
6S 入力ファイルのカスタマイズ (input2.dat)
1:
2:
0
(User-defined geometric conditions)
22.84 128.87 22.90 99.73 5 19 (SZA, SAZ, VZA, VAZ, month, day)
3:
4:
3
11
(Atmosphere model: Midlatitude Winter)
(Sun Photometric Distribution)
5:
6:
7:
20
1.209e-02 3.652e-10
1.768e-02 2.202e-09
(Number of data records)
(Radius in um, dV/d(logr) in cm3/cm2)
(Radius in um, dV/d(logr) in cm3/cm2)
(中略)
24: 1.131e+01 3.809e-06
(Radius in um, dV/d(logr) in cm3/cm2)
25:
26:
27:
28:
29:
30:
1.654e+01 3.083e-06
1.50 1.50 1.50 1.50 1.50 1.50
7.225e-04 5.000e-05 8.408e-04
0
0
0.2
(Radius in um, dV/d(logr) in cm3/cm2)
1.50 1.50 1.50 1.50 1.50 1.50 1.50 1.50 · · ·
2.884e-03 4.663e-03 5.849e-03 6.265e-03 · · ·
(No results saved)
(Constant value for tau)
(Aerosol optical thickness at 550 nm tau)
31:
32:
33:
0
-1000
1
(Target level)
(Sensor level)
34:
0.530 0.560
35:
36:
0.0000 0.0018 0.0020 0.0130 0.1792 0.6463 0.9806 0.9950 0.8617 0.4262 · · ·
0
(Homogeneous surface)
(User-defined filtered function)
(Min. and Max. wavelength in um)
(後略)
以前の入力ファイルとの違いは以下の通り。
• 4 行目: エアロゾルモデルに 11(Sun Photometric Distribution) を指定
• 5 行目: 体積分布のデータ数を指定
• 6∼25 行目: 粒半径 (µm)、体積分布 (cm3 /cm2 /µm)
• 26∼27 行目: 以下の 20 波長 (µm) における複素屈折率実部、虚部
0.3500, 0.4000, 0.4120, 0.4430, 0.4700, 0.4880, 0.5150, 0.5500, 0.5900, 0.6330, 0.6700, 0.6940,
0.7600, 0.8600, 1.2400, 1.5360, 1.6500, 1.9500, 2.2500, 3.7500
• 33 行目: センサー波長特性に 1 (ユーザー関数) を指定
• 34 行目: 波長の最小値・最大値 (µm) を指定
• 35 行目: 波長 2.5 nm 間隔のフィルター関数
5.4. MODIS 衛星画像の大気補正
39
MODIS 衛星画像の大気補正
5.4
これまでに学んだことを統合して、MODIS 衛星画像の大気補正にチャレンジしよう。使用する MODIS
衛星画像は、千葉大周辺に雲がなく、かつ比較的光学的厚さが厚いと思われる日のものである。簡単の
ため、太陽方向、視線方向には千葉大にもっとも近いピクセルの値を使う。また、エアロゾルの光学的
厚さには千葉大にあるスカイラジオメーターで測定された値を使う。これら (太陽方向、視線方向、エ
アロゾル光学的厚さ) の値が一様であるという仮定は、実際には千葉大のごく周辺にしか当てはまらな
いことに注意する。
以下に大気補正の手順を示す。
1. MODIS バンド 12 放射輝度データ (MOD021KM)、Geolocation データ (MOD03) を読み込む
2. 千葉大 (北緯 35.622 度、東経 140.124 度) にもっとも近いピクセルを探す
3. 衛星観測時刻、緯度、経度、太陽方向、視線方向を求める
4. スカイラジオメーターの光学特性データから、衛星観測時刻に最も近いデータを選ぶ
5. スカイラジオメーターの物理特性データから、衛星観測時刻に最も近いデータを選ぶ
6. 6S を実行し、大気補正係数 (xa、xb、xc) を求める
7. 千葉大周辺領域について、Lobs 、τ から ρac を求める
これらは script ディレクトリにある atmos_correct.py というスクリプトにまとめてある。data ディ
レクトリで IPython を立ち上げ、%run ../script/sixs_atmos_correct.py を実行すると下のような
図が得られる。
36.0
Apparent Reflectance (MODIS Band 12)
Atmospherically Corrected Reflectance (MODIS Band 12) 0.20
36.0
0.20
0.18
35.8
35.6
0.18
35.8
0.16
Chiba-U
0.14
35.6
0.16
Chiba-U
0.10
0.08
35.2
0.06
0.04
35.0
35.4
0.00
0.10
0.08
35.2
35.0
0.02
34.8 139.6 139.8 140.0 140.2 140.4 140.6 140.8 141.0
Longitude
0.14
0.12
Latitude
Reflectance
Latitude
0.12
35.4
Reflectance
0.06
0.04
0.02
34.8 139.6 139.8 140.0 140.2 140.4 140.6 140.8 141.0
Longitude
0.00
図 5.6: 千葉周辺の地表面反射率 (MODIS バンド 12)
左：大気補正前 (ρ∗i )、右：大気補正後 (ρac )
第 5 章 6S を使った MODIS 衛星画像の大気補正
40
6S カスタマイズモデル実行関数 (IPython 入力 or sixs_run_custom.py)
def run_6s(sza,saz,vza,vaz,mon,day,wlen,rsr,radius,volume,refr,refi,tau):
lines = ""
lines += "0\n" # User defined
lines += "%.10f %.10f %.10f %.10f %d %d\n"%(sza,saz,vza,vaz,mon,day)
lines += "3\n" # Midlatitude Winter
lines += "11\n" # Sun Photometric Distribution
lines += "%d\n"%(radius.size) # Number of data records
for r,v in zip(radius,volume):
lines += "%.3e %.3e\n"%(r,v) # Radius in um, dV/d(logr) in cm3/cm2
for r in refr:
lines += " %4.2f"%(r)
lines += "\n"
for r in refi:
lines += " %4.2f"%(r)
lines += "\n"
lines += "0\n" # No results saved
lines += "0\n" # Constant value for tau
lines += "%.10f\n"%(tau) # tau
lines += "0\n" # Target level
lines += "-1000\n" # Sensor level
lines += "1\n" # User ’s defined filtered function
lines += "%.3f %.3f\n"%(wlen[0],wlen[-1]) # Min. and Max. wavelength in um
for r in rsr:
lines += " %.4f"%(r)
lines += "\n"
lines += "0\n" # Homogeneous surface
lines += "0\n" # No directional effects
lines += "0\n" # Constant value for rho
lines += "0.1\n" # rho
lines += "0\n" # Atm. correction Lambertian
lines += "50.0\n" # Radiance (positive value)
p = Popen("sixsV1.1",stdin=PIPE,stdout=PIPE,stderr=PIPE,shell=True)
out,err = p.communicate(lines)
xa = nan
xb = nan
xc = nan
for line in out.splitlines():
m = re.search("coefficients\s+xa\s+xb\s+xc\s*:\s*(\S+)\s+(\S+)\s+(\S+)",
line)
if m:
xa = float(m.group(1))
xb = float(m.group(2))
xc = float(m.group(3))
return xa,xb,xc
5.4. MODIS 衛星画像の大気補正
演習問題
• MOD03 データ (MOD03.J20120519011621.20120519012656.hdf) から千葉周辺 (緯度 34.87
∼35.74 度、経度 139.60∼140.60 度) における太陽方向 (SolarZenith、SolarAzimuth)、視
線方向 (SensorZenith、SensorAzimuth) のデータを読み出し、各々の最大値、最小値、平
均値、標準偏差を計算してみよう。
• 余裕がある人は太陽天頂角、視線天頂角等の緯度経度マップを作成してみよう。
• 千葉大から離れた地点 (例えば緯度 34.87 度、経度 139.60 度) における太陽方向、視線方向
の値を 6S に入力して大気補正をしてみよう。千葉大の値を使って得た値とどう違うか？
41
43
付 録A
Python プログラミング簡易リファレンス
ここで、テキスト本文であまり説明しなかった Python プログラミングの方法を簡単にまとめる。
本文では IPython に--pylab オプションを与えて Matplotlib を簡単に使えるようにしていたが、ここ
では Numpy、Matplotlib のモジュールを以下のように読み込む (import する)。
In : import numpy as np
In : import matplotlib.pyplot as plt
こうすると、Numpy、Matplotlib のメソッドはそれぞれ np. メソッド名、plt. メソッド名のように書
かれるため、各メソッドがどのモジュールに属するか明確になる。
A.1
配列
複数の要素を持つオブジェクトを配列と呼ぶこと
にすると、Python にはいろいろな種類の配列 (リ
スト、タプル、文字列といったシーケンス型に分類
されるオブジェクト、Numpy の array オブジェク
ト、辞書のようなマップ型のオブジェクト等) が用
意されている。これらには共通することも多いの
で、まず代表的なリストについて述べる。
A.1.1
リスト
Python のリストは”[]”の中に要素を”,” で区切
って並べることで作成出来る。任意のオブジェクト
を要素にすることが出来、異なる種類のオブジェク
トが混ざっていても構わない。リストの要素数は
len という組み込み関数で調べられる。リストの
中に特定の要素が含まれるかどうかは in という演
算子で、また、最初の要素番号 (インデックス) は
index というメソッドで調べられる。
使用例:
In : a = [0,1,2,3,"four",5]
In : a
Out: [0, 1, 2, 3, ’four’, 5]
In : len(a)
Out: 6
In : 2 in a
Out: True
In : 4 in a
Out: False
In : a.index("four")
Out: 4
• リスト要素の指定
リストの要素は、演算子”[]” を使って一部を切
り出すことが出来る。(これをスライスと言う。)
スライスに代入を行ってリストの要素を変更する
ことも出来る。この際、要素数も変えられる。
使い方:
l[start=0[:stop=len(l)][:step=1]]
ã リスト l の start から stop まで step 間隔の
要素を指定する。stop は含まれないことに注意。
start や stop に負の値を入れると要素を後ろから
指定出来る。start 等は省略することも出来、単
に”:”だけ入れると全ての要素を指定出来る。
使用例:
In :
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
In :
Out:
a = [0,1,2,3,"four",5]
a[4]
’four’
a[-1]
5
a[2:5]
[2, 3, ’four’]
a[-2:1:-2]
[’four’, 2]
a[:2]
[0, 1]
a[::2]
[0, 2, ’four’]
a[:]
[0, 1, 2, 3, ’four’, 5]
a[4:] = [4]
a
[0, 1, 2, 3, 4]
• range 関数
使い方:
range([start=0,]stop[,step=1])
ã start から stop まで step 間隔の数列 (リスト)
を作る。stop は含まれないことに注意。
付 録A
44
使用例:
In :
Out:
In :
Out:
In :
Out:
range(10)
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
range(3,10)
[3, 4, 5, 6, 7, 8, 9]
range(3,10,2)
[3, 5, 7, 9]
• リストの結合、繰り返し
演算子”+”を使ってリストを結合したり、演算子”*”
を使ってリストの要素を繰り返した別のリストを作
ることができる。( 元のリストは変更されない。 )
使用例:
In :
In :
Out:
In :
Out:
In :
In :
Out:
In :
Out:
a = range(3)
a
[0, 1, 2]
a*2
[0, 1, 2, 0, 1, 2]
b = range(10,13)
b
[10, 11, 12]
a+b
[0, 1, 2, 10, 11, 12]
• リストの append メソッド
リストは変更可能なオブジェクト であり、append
メソッドを使って要素を追加することが出来る。
使用例:
In :
In :
Out:
In :
In :
Out:
a = range(3)
a
[0, 1, 2]
a.append("three")
a
[0, 1, 2, ’three’]
• リストの内包表記
式の後ろに for 節や if 節を付けることでリストを
作成することが出来る。後述の for 文を使うより
も高速に動作する。
使用例:
In :
Out:
In :
Out:
[i*3 for i in range(3)]
[0, 3, 6]
[i*4 for i in range(10) if i%2 == 0]
[0, 8, 16, 24, 32]
• リストのコピー
オブジェクトをコピーする際、やり方によっては
同じ (ID が等しい) オブジェクトを参照する別名を
Python プログラミング簡易リファレンス
作ってしまうことになるので注意が必要。 ID が等
しいかどうかは”is” という演算子で調べられる。
使用例:
In : a = range(3)
In : a
Out: [0, 1, 2]
In : b = a # b is an "alias" of a
In : b
Out: [0, 1, 2]
In : b is a
Out: True
In : b.append("three")
In : b
Out: [0, 1, 2, ’three’]
In : a
Out: [0, 1, 2, ’three’]
In : a = range(10,20,2) # redefine a
In : a
Out: [10, 12, 14, 16, 18]
In : b
Out: [0, 1, 2, ’three’]
In : b is a
Out: False
In : b = a[:] # b is a copy of a
In : b
Out: [10, 12, 14, 16, 18]
In : b is a
Out: False
In : a = [[]]*3 # 3 []s are identical
In : a
Out: [[], [], []]
In : a[0].append(1)
In : a
Out: [[1], [1], [1]]
In : a = [[] for i in range(3)] # different []s
In : a
Out: [[], [], []]
In : a[0].append(1)
In : a
Out: [[1], [], []]
A.1.2
タプル
タプルは”()”の中に要素を”,” で区切って並べ
ることで作成出来る。一度作成されたタプルは要素
の変更が出来ない。 (immutable と言う。) 要素の
指定、結合、繰り返し等はリストと同様に出来る。
A.1. 配列
使用例:
In : a = (0,1,2,3,"four",5)
In : a[4]
Out: ’four’
In : a[4] = 4 # -> error
TypeError
• タプルの pack、unpack
複数のオブジェクトを 1 つのタプルにまとめること
を pack、その逆を unpack と言う。unpack はタプ
ルの要素数と与えられた変数の個数が等しくない
と出来ない。 unpack はリスト等でも可能である。
使用例:
In : t = 1,2,3,"OK" # pack
In : t
Out: (1, 2, 3, ’OK’)
In : a,b,c,d = t # unpack
In : a
Out: 1
In : d
Out: ’OK’
45
Out:
In :
In :
Out:
In :
In :
Out:
[’hello’, ’world!’]
s = "abC"
s.upper(),s.lower(),s.swapcase()
(’ABC’, ’abc’, ’ABc’)
s = "ABC"
s.ljust(5),s.rjust(5),s.center(5)
(’ABC ’, ’ ABC’, ’ ABC ’)
A.1.4
array (numpy.ndarray)
Numpy モジュールの array オブジェクトは数値
計算に適した多次元配列である。ベクトルや行列の
ように、array は個々の要素を意識せず、あたかも
1 つの変数のように一括処理出来る。 リストと異
なり、
• 全要素は同じ型のオブジェクト
• 同じ次元の要素数は等しい
(Numpy では次元のことを axis と呼ぶ。)
• 全要素数の変更は効率的でない
(全要素数 · · · 全 axis の要素数の和)
という制約があるが、その代わりに全要素に渡る
高速処理が可能 である。上の条件を満たすリスト
A.1.3 文字列
は array に変換することが出来る。各 axis の要
文字列も配列の一種であり、ダブルクォーテーショ 素数は shape、全要素数は size という属性を調べ
ン ("") またはシングルクォーテーション (’’) の中 ると分かる。(shape は各 axis の要素数を成分に
に文字を並べて作成出来る。改行 (\n)、タブ (\t) 持つタプルである。) axis の順序は swapaxes メ
等、一部の特殊文字は”\”+文字列で表す。これをエ ソッドで入れ替えられる。全要素数が同じであれ
スケープシーケンスと言う。文字列の前に r を付 ば、reshape メソッドで shape を変更出来る。
けると raw 文字列になり、エスケープシーケンス 使用例:
が無効になる。raw 文字列は後述の正規表現で”\” In : a = np.array([0.0,1.0,2.0,3.0])
を入力する時に便利である。文字列の前に u を付 In : a
けると Unicode 文字列になり、\u0020 (空白文字) Out: array([ 0., 1., 2., 3.])
のような Unicode が扱える。一度作成された文字 In : a.shape
列は、タプル同様、要素の変更が出来ない。 要素 Out: (4,)
の指定、結合、繰り返し等はリストと同じように出 In : a.size
来る。文字列オブジェクトは文字列処理に便利な多 Out: 4
くのメソッドを持っている。
In : a = np.array([[0.0,1.0],[2.0,3.0]])
使用例:
In : a
Out:
In : "
hello world! \n".strip()
array([[ 0., 1.],
Out: ’hello world!’
[ 2., 3.]])
In : "a1\ta2 a3\na4".split()
In : a.shape
Out: [’a1’, ’a2’, ’a3’, ’a4’]
Out: (2, 2)
In : "a1,a2,a3".split(",")
In : a.size
Out: [’a1’, ’a2’, ’a3’]
In : "a1, a2, a3".replace(","," ").split() Out: 4
In : a.swapaxes(0,1)
Out: [’a1’, ’a2’, ’a3’]
Out:
In : "hello\nworld!\n".splitlines()
46
付 録A
array([[ 0., 2.],
[ 1., 3.]])
In : a.reshape(4)
Out: array([ 0., 1., 2., 3.])
In : a.reshape(1,4)
Out: array([[ 0., 1., 2., 3.]])
In : a.reshape(4,1)
Out:
array([[ 0.],
[ 1.],
[ 2.],
[ 3.]])
• empty、zeros、ones メソッド
これらを使うと指定した shape の array を作成
出来る。zeros、ones は内容をそれぞれ 0、1 に
初期化するのに対し、empty は初期化を行わない。
(empty で作成された array の内容は不定。)
使い方:
empty(shape, dtype=float, order=’C’)
ã zeros、ones も同様。shape はタプルでもリス
トでも良い。例えば dtype を int にすると整数型
の array になる。order を’F’ にすると Fortran
形式 (前方の添字が先に回る) の array になる。
使用例:
In : np.empty(3)
Out: array([ 6.93501252e-310,
6.93501252e-310,
1.58101007e-322])
In : np.zeros(4)
Out: array([ 0., 0., 0., 0.])
In : np.ones([2,2])
Out:
array([[ 1., 1.],
[ 1., 1.]])
• arange メソッド
Numpy の arange メソッドは整数に対する range
関数を実数に拡張したようなものである。
使い方:
arange([start=0.0,]stop[,step=1.0])
ã start から stop まで step 間隔の数列 (array)
を作る。stop は含まれないことに注意。
使用例:
In : np.arange(5.0)
Out: array([ 0., 1., 2., 3., 4.])
In : np.arange(3.5,4.6,0.5)
Out: array([ 3.5, 4. , 4.5])
Python プログラミング簡易リファレンス
• linspace、logspace メソッド
使い方:
linspace(start,stop,num=50)
logspace(start,stop,num=50,base=10.0)
ã linspace はリニアスケール、logspace は基底
base のログスケールで start から stop まで等間
隔の数列 (array) を作る。要素数は num で指定出
来る。デフォルトで stop は含まれる。
使用例:
In :
Out:
In :
Out:
np.linspace(1.0,4.0,4)
array([ 1., 2., 3., 4.])
np.logspace(1.0,4.0,4)
array([ 10., 100., 1000.,10000.])
• array 要素の指定
array の要素はリストのようにスライスを使って
指定出来るほか、インデックスの array (リストも
可) や、bool 型の array (リストは不可) を使って
指定出来る。 指定要素への代入で、内容の変更は
可能だが、要素数は変更出来ない。bool 型 (True
or False) の array は”>”、”<”、”==”、”|”、”&”
のような演算子で作られる。bool 型の array に
True が含まれるかどうかは、any メソッドや sum
メソッドを使って調べられる。where メソッドに
bool 型の array を与えると True がある位置のイ
ンデックスが得られる。(axis ごとのインデックス
の array がタプルとして得られる。)
使用例:
In :
In :
Out:
In :
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
In :
Out:
In :
Out:
a = np.array([10.0,11.0,12.0,13.0])
a[[0,1,3]]
array([ 10., 11., 13.])
cnd = (a < 11.5) | (a > 12.5)
cnd
array([ True, True, False, True],
dtype=bool)
np.any(cnd)
True
cnd.sum() # True = 1, False = 0
3
a[cnd]
array([ 10., 11., 13.])
indx = np.where(cnd)
indx
(array([0, 1, 3]),)
a[indx[0]]
array([ 10., 11., 13.])
A.1. 配列
47
• array の和、積
リストと異なり、array の”+”、”*” はそれぞれ要
素ごとの和、積を表す。(行列の積は dot メソッド
で計算出来る。)
使用例:
In : a = np.arange(4.0)
In : a
Out: array([ 0., 1., 2.,
In : a+a
Out: array([ 0., 2., 4.,
In : a*2
Out: array([ 0., 2., 4.,
In : a = a.reshape(2,2)
In : a
Out:
array([[ 0., 1.],
[ 2., 3.]])
In : a*a
Out:
array([[ 0., 1.],
[ 4., 9.]])
In : np.dot(a,a)
Out:
array([[ 2.,
3.],
[ 6., 11.]])
3.])
6.])
6.])
• array の結合
concatenate、append、vstack、hstack といった
メソッドを使って array を結合することが出来る。
(元の array は変更されない。)
使い方:
concatenate(tup,axis=0)
append(a1,a2,axis=None)
vstack(tup)
hstack(tup)
ã tup· · · 結合する array のタプル、a1,a2· · · 結合
する array、axis· · · 結合する方向、axis が None
の場合、flat な array となる。
使用例:
array([[ 20., 21.],
[ 22., 23.]])
In : np.concatenate((a,b))
Out:
array([[ 10., 11.],
[ 12., 13.],
[ 20., 21.],
[ 22., 23.]])
In : np.concatenate((a,b),axis=1)
Out:
array([[ 10., 11., 20., 21.],
[ 12., 13., 22., 23.]])
In : np.append(a,b)
Out: array([ 10., 11., 12., 13.,
20., 21., 22., 23.])
In : np.vstack((a,b))
Out:
array([[ 10., 11.],
[ 12., 13.],
[ 20., 21.],
[ 22., 23.]])
In : np.hstack((a,b))
Out:
array([[ 10., 11., 20., 21.],
[ 12., 13., 22., 23.]])
• その他
Numpy は他にも多くのメソッドを提供している。
(www.scipy.org/Numpy Functions by Category
参照。) Numpy のメソッドの中には、デフォル
トで array の全要素をまとめて扱うが、引数で
axis を指定すると axis ごとにまとめて扱うもの
もある。
使用例:
In : a = np.arange(4.0).reshape(2,2)
In : a
Out:
In : a = np.arange(10.0,14.0).reshape(2,2) array([[ 0., 1.],
[ 2., 3.]])
In : b = np.arange(20.0,24.0).reshape(2,2)
In : a.mean()
In : a
Out: 1.5
Out:
In : a.mean(axis=0)
array([[ 10., 11.],
Out: array([ 1., 2.])
[ 12., 13.]])
In : a.mean(axis=1)
In : b
Out: array([ 0.5, 2.5])
Out:
付 録A
48
A.1.5
Python プログラミング簡易リファレンス
辞書
In : w.update({"CO2":44})
In : for molec,mass in w.iteritems():
:
print molec,mass
例えば窒素 (N2 )、酸素 (O2 )、二酸化炭素 (CO2 )
:
の分子量をテーブルに入れてプログラムで参照す
る場合、w = [28,32,44] というリストを作って N2 28
w[0]、w[1]、w[2] のように参照しても良いが、 O2 32
w["N2"]、w["O2"]、w["CO2"] のように参照出来 CO2 44
た方が分かりやすいし、バグの抑制にもなる。これ
を可能にするのが辞書型オブジェクトである。
辞書は”{}” の中に要素を”,” で区切って並べる
ことで作成出来る。各要素は key と value の組で
あり、key:value のように書く。update メソッド
で新たな要素を追加出来る。(key が既存の場合は
value が更新される。) keys、values、items とい
うメソッドで、それぞれ key、value、(key,value)
のリストが得られる。iteritems メソッドと後述
の for 文を使うと key、value を順番に取り出せ
る。
使用例:
In : w = {"N2":28,"O2":32}
In : w
Out: {’N2’: 28, ’O2’: 32}
In : w["O2"]
Out: 32
In : w.keys()
Out: [’N2’, ’O2’]
In : w.values()
Out: [28, 32]
In : w.items()
Out: [(’N2’, 28), (’O2’, 32)]
In : w.update({"CO2":44})
In : w
Out: {’CO2’: 44, ’N2’: 28, ’O2’: 32}
In : for molec,mass in w.iteritems():
:
print molec,mass
:
N2 28
CO2 44
O2 32
上の例のように、辞書の順番は作成時の順番に無関
係である。 順番を保つには、OrderedDict を使う
と良い。
使用例:
from
In :
In :
In :
collections import OrderedDict
w = OrderedDict()
w.update({"N2":28})
w.update({"O2":32})
A.2
制御構造
• if 文
使い方:
if 条件:
␣␣␣␣処理内容
ã 上 の”␣”は 空 白 を 表 し て い る 。こ の よ う に 、
Python はブロックの区切りに括弧を使わず、イン
デントを使う。すなわち、同じブロックの行は同じ
インデントで揃える。 このようにブロック区切り
にインデントを使うことで制御構造が視覚的に明
確になる。インデントの量は空白 4 個分くらいが適
当である。例外的に、()、[]、{}のような括弧の
途中で改行した場合は任意のインデントで継続出
来る。if 文が終わった後に行を続ける場合、if の
行と同じインデントでなければならない。残念なが
ら、ブロックの終わりを機械的に判断することは不
可能なので、インデントを戻すのは手動で行う必要
がある。条件が偽の場合、else 節があればそれが
実行される。elif 節は else 節に if 文を置くのと
同じであるが、インデントを減らすことが出来る。
使用例:
In : i = 20
In : if i > 100:
:
print 100
: elif i > 10:
:
print 10
: else:
:
print 1
:
10
• for 文
使い方:
for 変数 in 配列:
␣␣␣␣処理内容
ã 途中でループを抜ける場合は break 文、途中
でループ先頭に戻る場合は continue 文を使う。
Python の for 文は、他の大抵のインタプリタ型
A.3. 関数
49
スクリプト言語同様、あまり効率的ではないので、 A.3
関数
ループ回数が多い場合は使用を避けた方が良い。
(リストの内包表記や Numpy 等のライブラリを使 定義方法:
def 関数名 (引数):
えば多くの場合 for 文の使用を避けられる。)
␣␣␣␣処理内容
使用例:
␣␣␣␣return [戻り値]
ã 引数には”=”に続けてデフォルト値を与えること
In : for i in range(10):
が出来る。(デフォルト値を持つ引数は持たない引
:
if i > 8:
数の後ろに置く。) 戻り値にはゼロ個以上の複数の
:
break
オブジェクトを指定出来る。
:
elif i < 4:
呼び出し方法:
:
continue
[変数名=] 関数名 (引数)
:
else:
ã 変数名を与えると戻り値が格納される。(戻り値
:
print i
が複数ある場合、変数名を 1 つ与えると戻り値のタ
:
プルが得られる。戻り値の個数に等しい変数名を与
4
えるとタプルが unpack される。付録 A.1.2 参照。)
5
デフォルト値がある引数は省略出来る。引数名=値
6
の形 (キーワード引数) でも与えられる。(キーワー
7
ド引数は非キーワード引数の後ろに置く。) ”*”演
8
算子に続けて引数を入れたリストやタプルを与え
たり、”**”演算子に続けてキーワード引数を入れ
た辞書を与えることも出来る。
• while 文
使用例:
使い方:
In : def get_limit(a):
while 条件:
:
return min(a),max(a)
␣␣␣␣処理内容
:
ã 途中でループを抜ける場合は break 文、途中
In : x = [3,5,8,2]
でループ先頭に戻る場合は continue 文を使う。
In : xmin,xmax = get_limit(x)
前述の for 文同様、ループ回数が多い場合は使用
In : xmin,xmax
を避けた方が良い。
Out: (2, 8)
使用例:
In : def print_ab(a,b="world"):
:
print a,b
In : i = 0
:
In : while i < 10:
In : print_ab("hello")
:
if i > 8:
hello world
:
break
In : print_ab(*["Hello","World"])
:
elif i < 4:
Hello World
:
i += 1
In : print_ab(**{"a":"Hello,","b":"World!"})
:
continue
Hello, World!
:
else:
:
print i
• 変数の有効範囲 (スコープ)
:
i += 1
関数内で定義した変数 (ローカル変数) は、関数外
:
で使用出来ない。関数外で定義した変数は、関数内
4
で参照出来るが、普通に代入しようとするとローカ
5
ル変数が定義されるため、関数外の変数は変更され
6
ない。(global 文を用いてグローバル変数を定義出
7
来るが、あまり推奨されない。)
8
付 録A
50
A.4
例外処理
使い方:
try:
␣␣␣␣処理内容
except 例外名 [, 変数名]:
␣␣␣␣例外処理内容
ã 例外名の後ろに変数名を置くと例外インスタン
スを取得出来る。例外インスタンスの中身を調べる
と例外の詳細が分かる。一番後ろに else 節を置い
て例外が発生しなかった時の処理を指定出来る。
使用例:
In : try:
:
0.0/0.0
: except ZeroDivisionError:
:
print "Divided by zero"
:
Divided by zero
In : try:
:
raise ValueError("Error example")
: except ValueError,inst:
:
print inst
:
Error example
A.5
正規表現
(Regular Expression)
正規表現とは、通常文字と特殊文字を組み合わせ
て出来る表現である。これは対象の文字列に特定の
パターンが含まれるかどうか調べるために用いら
れる。Python で使える主な特殊文字を表 A.1 に挙
げる。
• re モジュール
使い方:
re.search(pattern,string)
re.sub(pattern,repl,string[,count=0])
ã pattern· · · 検索パターン、string· · · 対象文字
列、repl· · · 変換文字列、最大置換数は count で
指定でき、これが 0 であれば全てが置換される。
使用例:
In : import re
In : for item in ["N2","O2","CO2"]:
In :
if re.search("O2",item):
In :
print item
:
O2
CO2
In :
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
Out:
Python プログラミング簡易リファレンス
m = re.search("(\d+)(\D+)(\d+)","18O2")
m.group(1)
’18’
m.group(2)
’O’
m.group(3)
’2’
re.sub("(\d+)(\D+)(\d+)",
r"$^{\1}$\2$_{\3}$", # r -> raw string
"18O2")
’$^{18}$O$_{2}$’
re.sub("(\d+)",r"$_{\1}$","Na2CO3")
’Na$_{2}$CO$_{3}$’
re.sub("(\d+)",r"$_{\1}$",
"O2 32 g/mol",count=1)
’O$_{2}$ 32 g/mol’
表 A.1: Python で使える正規表現の主な特殊文字
特殊文字
意味
"."
改行以外の任意の文字
^
$
*
文字列の先頭
文字列の末尾
0 回以上の繰り返し (最大一致)
+
?
1 回以上の繰り返し (最大一致)
0 回か 1 回の繰り返し (最大一致)
*?, +?, ??
{m}
{m, n}
対応する特殊文字の最小一致
m 回の繰り返し
m ∼ n 回の繰り返し (最大一致)
{m, n}?
|
m ∼ n 回の繰り返し (最小一致)
前後のどちらか
()
(?:)
[]
後方参照に使えるグループ
[^]
\d
補集合
\D
\w
\W
\s
\S
単純なグループ
集合
数字 ([0-9])
非数字
英数文字 ([a-zA-Z0-9_])
非英数文字
空白文字 ([\t\n\r\f\v])
非空白文字
A.6. 日付・時刻
A.6
日付・時刻
51
• datetime モジュール
使用例:
In :
:
• time モジュール
In
:
time モジュールは C 言語の time.h にある関数と
In :
同様の関数を提供する。
Out:
使用例:
In :
In :
In :
Out:
In :
In :
Out:
In :
Out:
In :
Out:
import time
t = time.time()
t
1345440277.0
tl = time.localtime(t)
tl
time.struct_time(
tm_year=2012, tm_mon=8, tm_mday=20,
tm_hour=14, tm_min=24, tm_sec=37,
tm_wday=0, tm_yday=233, tm_isdst=0)
time.mktime(tl)
1345440277.0
time.gmtime(t)
time.struct_time(
tm_year=2012, tm_mon=8, tm_mday=20,
tm_hour=5, tm_min=24, tm_sec=37,
tm_wday=0, tm_yday=233, tm_isdst=0)
以下のように環境変数 TZ を設定するとスクリプト
内部におけるローカルタイムを変更出来る。
In :
In :
In :
In :
In :
Out:
In :
In :
In :
Out:
import os
os.environ["TZ"] = "UTC"
time.tzset()
t = 1345440277.0
time.localtime(t)
time.struct_time(
tm_year=2012, tm_mon=8, tm_mday=20,
tm_hour=5, tm_min=24, tm_sec=37,
tm_wday=0, tm_yday=233, tm_isdst=0)
os.environ["TZ"] = "JST-9"
time.tzset()
time.localtime(t)
time.struct_time(
tm_year=2012, tm_mon=8, tm_mday=20,
tm_hour=14, tm_min=24, tm_sec=37,
tm_wday=0, tm_yday=233, tm_isdst=0)
In :
Out:
In :
Out:
In :
In :
In :
Out:
In :
Out:
In :
Out:
In :
Out:
In :
In :
Out:
In :
In :
Out:
In :
Out:
from datetime import datetime,\
timedelta
t1 = datetime.now()
t1
datetime.datetime(
2012, 8, 20, 14, 24, 37, 718808)
t1.year,t1.month,t1.day
(2012, 8, 20)
t1.hour,t1.minute,t1.second
(14, 24, 37)
t2 = datetime.strptime(
"2012-8-20T14:00:00",
"%Y-%m-%dT%H:%M:%S")
dt = t1-t2
dt
datetime.timedelta(0, 1477, 718808)
dt.days,dt.seconds,dt.microseconds
(0, 1477, 718808)
t1+timedelta(1)
datetime.datetime(
2012, 8, 21, 14, 24, 37, 718808)
t1+timedelta(seconds=3600)
datetime.datetime(
2012, 8, 20, 15, 24, 37, 718808)
t = time.mktime(t1.timetuple())
t
1345440277.0
t += t1.microsecond*1.0e-6
datetime.fromtimestamp(t)
datetime.datetime(
2012, 8, 20, 14, 24, 37, 718808)
datetime.utcfromtimestamp(t)
datetime.datetime(
2012, 8, 20, 5, 24, 37, 718808)
後 述 す る Matplotlib の plot メ ソッド 等 は
datetime オブジェクトを X 軸にとって時系列図
を作成出来る。
使用例 (IPython 入力):
t1 = datetime.fromtimestamp(1345440277.0)
x = [t1+timedelta(i) for i in range(12)]
y = np.random.rand(len(x))
しばらくプログラムの実行を停止させる場合、 fig = plt.figure()
sleep メソッドが使える。引数は浮動小数点で、1 秒 plt.plot(x,y)
以下のスリープにも対応している。
fig.autofmt_xdate()
In : time.sleep(0.1)
plt.grid()
付 録A
52
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
012
21 2
Aug
A.7
A.7.1
Aug
012
23 2
12
5 20
2
Aug
Aug
012
27 2
Aug
012
29 2
Aug
012
31 2
数値データ処理
1 次元補間
Python プログラミング簡易リファレンス
ã x1,y1· · · 入力グリッド、入力データ、x2· · · 出
力グリッド、入力グリッドより外側の出力グリッド
に対し、bounds_error が True ならばエラーを出
し、False ならば fill_value を返す。
使用例 (IPython 入力):
from scipy.interpolate import interp1d
x = np.arange(0.0,3.0,0.5)
y = np.sin(x)
x2 = np.arange(-1.0,4.01,0.1)
f = interp1d(x,y,bounds_error=False,
fill_value=0.0)
y2 = f(x2)
plt.figure()
plt.plot(x,y,"bo",ms=20)
plt.plot(x2,y2,"ro")
plt.ylim(-0.2,1.2001)
plt.grid()
1.2
• interp (直線補間)
使い方:
interp(x2,x1,y1,left=None,right=None)
ã x2· · · 出力グリッド、x1,y1· · · 入力グリッド、入
力データ、left,right· · · 入力グリッドより左側の
値、右側の値
使用例 (IPython 入力):
x = np.arange(0.0,3.0,0.5)
y = np.sin(x)
x2 = np.arange(-1.0,4.01,0.1)
y2 = np.interp(x2,x,y)
plt.figure()
plt.plot(x,y,"bo",ms=20)
plt.plot(x2,y2,"ro")
plt.ylim(-0.2,1.2001)
plt.grid()
1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.2 1
0
1
2
3
4
• interp1d (直線補間)
使い方:
interp1d(x1,y1,bounds_error=True,
fill_value=np.nan)(x2)
1.0
0.8
0.6
0.4
0.2
0.0
0.2 1
0
1
2
3
4
• interp_linear (直線補間)
上に挙げた interp、interp1d はいずれも外挿が
出来ない。 interp_linear は、interp1d を基に、
外挿出来るように修正したものである。
def interp_linear(x1,y1,x2,
dmax=None,last=False):
indx = np.searchsorted(
x1,x2).clip(1,len(x1)-1).astype(int)
lo = indx-1
hi = indx
x_lo = x1[lo]
x_hi = x1[hi]
if last:
y_lo = y1[...,lo]
y_hi = y1[...,hi]
else:
y_lo = y1[lo]
y_hi = y1[hi]
r_lo = (x_hi-x2)/(x_hi-x_lo)
r_hi = 1.0-r_lo
y2 = r_lo*y_lo+r_hi*y_hi
A.7. 数値データ処理
53
1.5
if dmax is not None:
y2[x_hi-x_lo>dmax] = np.nan
return y2
1.0
0.5
ã x1,y1· · · 入力グリッド、入力データ、x2· · · 出
力グリッド、dmax が None でない場合、補間間隔
が dmax より広い出力グリッドに対して np.nan を
返す。last が False ならば最初の axis、True な
らば最後の axis に対して補間を行う。
使用例 (IPython 入力):
x = np.arange(0.0,3.0,0.5)
y = np.sin(x)
x2 = np.arange(-1.0,4.01,0.1)
y2 = interp_linear(x,y,x2)
plt.figure()
plt.plot(x,y,"bo",ms=20)
plt.plot(x2,y2,"ro")
plt.ylim(-1.0,1.5)
plt.grid()
0.0
0.5
1.0 1
A.7.2
0
1
2
3
4
2 次元補間
• griddata
使い方:
griddata(x1,y1,z1,x2,y2)
ã x1,y1,z1· · · 入力 X グリッド、入力 Y グリッド、
入力データ、これらは flatten メソッドを使って
1 次元データにしておく、x2,y2· · · 出力 X グリッ
ド、出力 Y グリッド
このメソッドではデフォルトで Robert Kern 氏に
よる Delaunay triangulation に基づいた natural
1.5
neighbor interpolation というアルゴリズムが用い
1.0
られる。ただし、NCAR natgrid library というラ
0.5
イブラリがインストールされていれば、それに基づ
いたより強固なアルゴリズムが用いられる。
0.0
使用例 (IPython 入力):
0.5
from matplotlib.pylab import griddata
1.0 1
x = np.arange(-3.0,3.1,0.5)
1
4
0
2
3
y = x[:]
x1,y1 = np.meshgrid(x,y)
z1 = np.exp(-0.5*(x1*x1+y1*y1))
• splrep,splev (スプライン補間)
xe = np.arange(-3.25,3.26,0.5)
使い方:
ye = xe[:]
tck = splrep(x1,y1,w=None)
plt.figure()
splev(x2,tck)
plt.pcolor(xe,ye,z1)
ã x1,y1· · · 入力グリッド、入力データ、w· · · 重み、 plt.xlim(-3.0,3.0)
x2· · · 出力グリッド
plt.ylim(-3.0,3.0)
使用例 (IPython 入力):
plt.colorbar()
from scipy.interpolate import splrep,splev x = np.arange(-3.0,3.1,0.1)
x = np.arange(0.0,3.0,0.5)
y = x[:]
y = np.sin(x)
x2,y2 = np.meshgrid(x,y)
x2 = np.arange(-1.0,4.01,0.1)
z2 = griddata(x1.flatten(),
y2 = splev(x2,splrep(x,y))
y1.flatten(),
plt.figure()
z1.flatten(),x2,y2)
plt.plot(x,y,"bo",ms=20)
cnd = [z2.mask == True]
plt.plot(x2,y2,"ro")
z2.data[cnd] = 0.0
plt.ylim(-1.0,1.5)
xe = np.arange(-3.05,3.06,0.1)
plt.grid()
ye = xe[:]
付 録A
54
Python プログラミング簡易リファレンス
120
plt.figure()
plt.pcolor(xe,ye,z2.data)
plt.xlim(-3.0,3.0)
plt.ylim(-3.0,3.0)
plt.colorbar()
100
Data
Fit
80
60
3
1.0
40
0.9
2
0.8
20
0.7
1
0.6
0
00
2
4
6
8
10
0.5
0.4
1
0.3
0.2
2
0.1
33
2
1
0
1
2
3
0.0
A.7.4
元データ
3
一般式最適化
1.0
0.9
2
0.8
0.7
1
0.6
0
0.5
0.4
1
0.3
0.2
2
0.1
33
2
1
0
1
2
3
0.0
補間データ
A.7.3
多項式最適化
• polyfit
使い方:
polyfit(x,y,deg)
ã x,y· · · 入力データ、deg· · · 次数
使用例 (IPython 入力):
x = np.arange(10.0)
y = np.random.normal(0.0,3.0,10)+x*x
result = np.polyfit(x,y,2)
x2 = np.arange(0.0,10.0,0.1)
y2 = np.polyval(result,x2)
plt.figure()
plt.plot(x,y,"bo",label="Data")
plt.plot(x2,y2,"r-",label="Fit")
plt.legend(loc=2,numpoints=1)
plt.grid()
• leastsq
使い方:
leastsq(func,x0,args=(),maxfev=0)
ã func· · · フィッティング関数、x0· · · 初期値、
args· · · func のフィッティングパラメータ以外の引
数、maxfev· · · func の最大呼び出し回数
このメソッドは内部で MINPACK ライブラリの
lmdif、lmder アルゴリズムを使い、LevenbergMarquardt 法によって非線型最小二乗問題を解く。
使用例 (IPython 入力):
from scipy.optimize import leastsq
def func(p,x):
return p[0]*x*x+p[1]*x+p[2]
def residuals(p,x,y):
return y-func(p,x)
x = np.arange(10.0)
y = np.random.normal(0.0,3.0,10)+x*x
p = np.array([1.0,0.0,0.0])
result = leastsq(residuals,p,args=(x,y))
p2 = result[0]
x2 = np.arange(0.0,10.0,0.1)
y2 = func(p2,x2)
plt.figure()
plt.plot(x,y,"bo",label="Data")
plt.plot(x2,y2,"r-",label="Fit")
plt.legend(loc=2,numpoints=1)
plt.grid()
A.8. ファイル入力
100
55
100
Data
Fit
80
80
60
60
40
40
20
20
00
2
4
6
8
10
• minuit2 モジュール
このモジュールは欧州原子核研究機構 (CERN) で
開発された MINUIT という関数最適化ツールを使
用する。利用には ROOT というパッケージが必要。
使用例 (IPython 入力):
import minuit2
def func(p,x):
return p[0]*x*x+p[1]*x+p[2]
def residuals(p0,p1,p2):
e = yinp-func([p0,p1,p2],xinp)
return (e*e).sum()
x = np.arange(10.0)
y = np.random.normal(0.0,3.0,10)+x*x
p = np.array([1.0,0.0,0.0])
xinp,yinp = x,y
m = minuit2.Minuit2(residuals)
m.printMode = 1
m.values["p0"] = 1.0
m.values["p1"] = 0.0
m.values["p2"] = 0.0
m.errors["p0"] = 0.5
m.errors["p1"] = 0.5
m.errors["p2"] = 0.5
m.maxcalls = 10000
m.tol = 10.0
m.migrad()
p0 = m.values["p0"]
p1 = m.values["p1"]
p2 = m.values["p2"]
x2 = np.arange(0.0,10.0,0.1)
y2 = func([p0,p1,p2],x2)
plt.figure()
plt.plot(x,y,"bo",label="Data")
plt.plot(x2,y2,"r-",label="Fit")
plt.legend(loc=2,numpoints=1)
plt.grid()
00
A.8
A.8.1
Data
Fit
2
4
6
8
10
ファイル入力
テキストファイル
• loadtxt
使い方:
loadtxt(fname,comments="#",
delimiter=None,skiprows=0,
usecols=None,unpack=False)
ã fname· · · ファイ ル 名 、comments· · · コ メ ン
トを表す文字、delimiter· · · 列の区切り文字、
skiprows· · · 読み飛ばす先頭の行数、usecols· · ·
読み込む列番号 (0 から始まる)、unpack を True に
すると読み込んだデータが転置され、列の数と同じ
array に分けることが出来る。
使用例:
例えば以下のようなテキストファイル (test.dat)
があるとする。
############################
# This is a sample text file
############################
0.0 10.0
1.0 11.0
2.0 12.0
3.0 13.0
############################
これは以下のコマンドで読み込める。
In : data = np.loadtxt("test.dat")
In : data
Out:
array([[ 0., 10.],
[ 1., 11.],
[ 2., 12.],
[ 3., 13.]])
付 録A
56
loadtxt は、ファイルからだけでなく、文字列か
らデータを読み込むことも出来る。
from StringIO import StringIO
for open("test.dat","r") as fp:
lines = fp.readlines()
data = np.loadtxt(StringIO(lines))
A.8.2
バイナリファイル
Python プログラミング簡易リファレンス
In
In
In
In
: import struct
: fmt = "4s2f4d"
: n = struct.calcsize(fmt)
: with open("test.dat","rb") as fp:
:
line = fp.read(n)
In : data = struct.unpack(fmt,line)
ここで、data は (header,p1,p2,p3,p4,p5,p6)
のようなタプルになっており、
In : header,p1,p2,p3,p4,p5,p6 = data
の よ う に 個々の パ ラ メ ー タ を 取 り 出 し た り、
In : pars = np.array(data[3:])
のようにデータを array として取り出すことが出
来る。(付録 A.1.2 参照。)
• struct モジュール
使い方:
struct.unpack(fmt, バイト列)
ã fmt· · · フォーマットを表す文字列 (表 A.2 参照)
フォーマット文字列の先頭に<、>を付けることでバ
イトオーダーをそれぞれリトルエンディアン (下位 A.8.3 画像ファイル
バイトから並べる方式)、ビッグエンディアン (上位 • imread
バイトから並べる方式) に指定出来る。
使い方:
imread(ファイル名)
表 A.2: 主なフォーマット文字の意味
ã 画像形式は拡張子から自動的に判断される。
使用例:
fmt C 言語の型
In : im = plt.imread("lena.png")
c
char
In : im.shape
B
unsigned char
Out: (480, 480, 4) # ny,nx,RGBA
s
char[]
In : im(0,0)
h
short
Out: array([ 0.89411765, 0.52941179,
H
unsigned short
0.50980395, 1.
],
i
int
dtype=float32)
I
unsigned int
In : plt.imshow(im[::-1,::-1])
f
float
d
double
使用例:
例えば以下のようなバイナリファイル (test.dat)
があるとする。
先頭位置
バイト数
fmt
パラメータ名
0
4
4
4
s
f
header
p1
8
12
20
4
8
8
f
d
d
p2
p3
p4
28
36
8
8
d
d
p5
p6
これは以下のコマンドで読み込める。
• Python Imaging Library (PIL)
使用例 (IPython 入力):
from PIL import Image
im = Image.open("lena.png")
im.crop((80,80,400,400)).\
transpose(Image.FLIP_LEFT_RIGHT).\
save("lena2.png")
A.9. 描画ライブラリ
57
from netCDF4 import Dataset
f = Dataset("test.nc","r")
f_attr = f.ncattrs()
v = f.variables["TEST"]
v_attr = v.ncattrs()
test = v[:]
A.8.4
描画ライブラリ
A.9
HDF ファイル
XY グラフ
A.9.1
HDF4
• PyHDF モジュール
使用例 (IPython 入力):
例えば TEST というフィールドを持つ HDF4 ファイ
ル (test.hdf) があるとする。
from pyhdf.SD import SD,SDC
f = SD(fnam,SDC.READ)
f_attr = f.attributes()
v = f.select("TEST")
v_attr = v.attributes()
test = v[:]
1.0
0.8
HDF5
plot(x,y[, オプション文字列,**kwargs])
ã x,y· · · 入力データ、オプション文字列には表 A.3
にまとめた色、線種、マーカーを組み合わせたもの
を与えることが出来る。color=色のキーワード引
数で表 A.3 には無い色 (DarkBlue、#0000A0 のよ
うな HTML カラーコード、0.0 ∼ 1.0 までのグレー
スケールを表す文字列) も指定出来る。plot には
表 A.4 に挙げた以外にも非常に多くのキーワード
引数がある。詳細はオンラインマニュアル参照。
表 A.3: plot のオプション文字列の意味
’b’
’g’
1.0
• H5py モジュール
0.6
使用例 (IPython 入力):
例えば TEST というフィールドを持つ HDF5 ファイ
0.4
ル (test.hdf) があるとする。
import h5py
0.2
f = h5py.File("test.hdf","r")
f_attr = f.attrs
0.0
v = f["TEST"]
0.0
0.2
v_attr = v.attrs
test = v[:]
A.8.5
• plot
使い方:
’r’
0.8
’c’
’-’
’m’
0.6
’--’
’y’
’-.’
0.4
’k’
’:’
(a) 色
(b) 線種
0.2
0.4
0.0
0.0
0.6
0.2
0.8
0.4
NetCDF ファイル
• netCDF4 モジュール
使用例 (IPython 入力):
例えば TEST というフィールドを持つ NetCDF ファ
イル (test.nc) があるとする。
(c) マーカー
0.6
1.0
0.8
1.0
付 録A
58
Python プログラミング簡易リファレンス
1.5
表 A.4: plot の主なキーワード引数
1.0
color or c
色
0.5
alpha
linestyle or ls
透明度
0.0
線種
0.5
linewidth or lw
marker
markeredgecolor or mec
線幅
1.0
マーカー
1.51.0
markeredgewidth or mew
markerfacecolor or mfc
マーカー線幅
マーカー色
A.9.2
markersize or ms
markevery
label
マーカーサイズ
• scatter
使い方:
scatter(x,y,s=20,c="b",marker="o",
cmap=None,vmin=None,vmax=None,
alpha=None)
ã x,y· · · 入力データ、s· · · マーカーサイズ、c· · ·
色、marker· · · マーカー (表 A.5 参照)、cmap· · ·
カ ラ ー マップ、vmin,vmax· · · 値 の 下 限 、上 限 、
alpha· · · 透明度 (0∼1)
マーカーサイズ、色には固定値の他、x や y と size
が等しい array を与えることが出来る。
0.5
0.0
0.5
1.0
マーカー線色
マーカー間隔
凡例に使うラベル
使用例 (IPython 入力):
x = np.arange(-1.0,1.0,0.01)
y = x*x*x
plt.figure()
plt.plot(x,y)
plt.plot(x,y,"ro",markevery=10)
plt.grid()
散布図
表 A.5: scatter に使えるマーカー
1.0
0.5
0.0
0.5
1.01.0
0.5
0.0
0.5
1.0
• errorbar
使い方:
errorbar(x,y,yerr=None,xerr=None,
fmt="-",ecolor=None)
ã x,y· · · 入力データ、xerr,yerr· · · x、y の誤差、
ecolor· · · 誤差棒の色、fmt に None を指定すると
誤差棒のみをプロット出来る。
使用例 (IPython 入力):
x = np.arange(-1.0,1.0,0.01)
y = x*x*x
e = np.ones_like(x)*0.5
plt.figure()
plt.plot(x,y,"bo-",markevery=10)
plt.errorbar(x[::10],y[::10],yerr=e[::10],
fmt=None,ecolor="b")
plt.grid()
表 A.5 の verts、path、tuple はそれぞれ Path
vertice、Path instance、タプル (number of sides,
style, angle) を表す。
表 A.6: タプルの style
値
スタイル
0
regular polygon
1
2
star-like symbol
asterisk
3
circle
表 A.5 では以下の値が使われている。
pathdata = [
(Path.MOVETO,(1.58,-2.57)),
(Path.CURVE4,(0.35,-1.1)),
(Path.CURVE4,(-1.75,2.0)),
A.9. 描画ライブラリ
59
(Path.CURVE4,(0.375,2.0)),
(Path.LINETO,(0.85,1.15)),
(Path.CURVE4,(2.2,3.2)),
(Path.CURVE4,(3,0.05)),
(Path.CURVE4,(2.0,-0.5)),
(Path.CLOSEPOLY,(1.58,-2.57)),
]
codes,verts = zip(*pathdata)
path = mpath.Path(verts,codes)
tuple = (5, 1, 30)
使用例 (IPython 入力):
A.9.4
ヒストグラム
• hist
使い方:
hist(x,bins=10,range=None,normed=False)
ã x· · · 入力データ、bins· · · 分割数、range· · · 入
力レンジ、normed が True の場合、積分値が 1 に
なるようにスケールが調整される。
使用例 (IPython 入力):
x = np.random.normal(size=1000)
plt.hist(x)
250
x = np.random.normal(0.0,1.0,100)
y = np.random.normal(0.0,1.0,100)+x
plt.figure()
plt.scatter(x,y,c=x,s=np.fabs(y)*20.0)
plt.grid()
200
150
100
50
4
3
03
2
1
2
1
0
2
3
4
1
0
1
A.9.5
2
2 次元プロット
3
• imshow
使い方:
54
1
1
3
2
0
2
3
imshow(z,cmap=None,interpolation=None,
vmin=None,vmax=None,
origin=None,extent=None)
A.9.3 棒グラフ
ã z· · · 入 力 デ ー タ、cmap· · · カ ラ ー マップ、
interpolation· · · 補 間 モ ー ド、vmin,vmax· · ·
• bar
値 の 下 限 、上 限 、origin· · · 上 下 の 基 準 、
使い方:
extent· · · X、Y の範囲
bar(left,height,width=0.8,bottom=0)
使用例 (IPython 入力):
ã left· · · 棒 の 左 端 、height· · · 棒 の 高 さ 、 x = np.arange(-4.0,4.1,0.5)
width· · · 棒の幅、bottom· · · 棒の下端
y = np.arange(-3.0,3.1,0.5)
使用例 (IPython 入力):
xi,yi = np.meshgrid(x,y)
x = arange(10.0)
zi = np.exp(-0.5*(xi*xi+yi*yi))
y = np.random.rand(x.size)
plt.imshow(zi,interpolation="none",
plt.figure()
extent=(x.min(),x.max(),
plt.bar(x,y)
y.min(),y.max()))
4
0.9
3
0.8
2
0.7
0.6
1
0.5
0
0.4
0.3
1
0.2
2
0.1
0.00
2
4
6
8
10
34
3
2
1
0
1
2
3
4
60
付 録A
Python プログラミング簡易リファレンス
4
• pcolor
3
使い方:
2
pcolor(x,y,z,cmap=None,vmin=None,vmax=None)
1
ã x,y· · · 入力グリッド、z· · · 入力データ、cmap· · ·
0
1
カラーマップ、vmin,vmax· · · 値の下限、上限
2
入力グリッドはグリッド中心ではなく、グリッド境
3
界を与える点に注意。
44
1
1
4
3
2
0
2
3
使用例 (IPython 入力):
n = 10
xo = np.arange(-4.0,4.0001,0.01)
A.9.6 地図投影
yo = np.arange(-4.0,4.0001,0.01)
xc = np.random.normal(0.0,0.1,n)
• Basemap モジュール
yc = np.random.normal(0.0,0.1,n)
使用例 (IPython 入力):
a = np.random.normal(0.02,0.01,n)
from mpl_toolkits.basemap import Basemap
b = np.random.normal(0.02,0.01,n)
lat_min,lat_max = 34.871,35.747
c = np.linspace(-3.0,3.0,n)
lon_min,lon_max = 139.633,140.621
hlines = []
lat = 35.622
vlines = []
lon = 140.124
for i in range(n):
m = Basemap(projection="cyl",
y2 = a[i]*np.square(xo-xc[i])+c[i]
lat_0=0.5*(lat_min+lat_max),
x2 = b[i]*np.square(yo-yc[i])+c[i]
lon_0=0.5*(lon_min+lon_max),
plt.plot(xo,y2,"b-")
llcrnrlat=lat_min,
plt.plot(x2,yo,"b-")
llcrnrlon=lon_min,
hlines.append(y2[:])
urcrnrlat=lat_max,
vlines.append(x2[:])
urcrnrlon=lon_max,
xp = []
resolution="h",
yp = []
area_thresh=1000.0)
for i in range(n):
m.drawcoastlines()
y2 = hlines[i]
m.drawparallels(arange(30.0,41.0,0.5),
for j in range(n):
labels=[1,0,0,0])
x2 = vlines[j]
m.drawmeridians(arange(130.0,150.0,0.5),
l = []
labels=[0,0,0,1])
for k in range(xo.size):
x,y = m(lon,lat)
l.append((np.square(xo[k]-x2)+
m.plot(x,y,"ro")
np.square(y2[k]-yo)).min())
plt.text(x,y,"Chiba-U")
l = np.array(l)
indx = np.argmin(l)
Chiba-U
xp.append(xo[indx])
35.5°N
yp.append(y2[indx])
xp = np.array(xp).reshape(n,n)
yp = np.array(yp).reshape(n,n)
xq = xp[:,:-1]+0.5*np.diff(xp,axis=1)
yq = yp[:-1,:]+0.5*np.diff(yp,axis=0)
35°N
xr = xq[:-1,:]+0.5*np.diff(xq,axis=0)
140°E
140.5°E
yr = yq[:,:-1]+0.5*np.diff(yq,axis=1)
zr = np.exp(-0.5*(xr*xr+yr*yr))
plt.pcolor(xp,yp,zr)
61
付 録B
本演習に必要なソフトウェア
本演習に必要なソフトウェアのうち、標準ではインストールされない可能性があるものについて、ソー
スの入手先とホームディレクトリへのインストール方法を簡単に示す。 は必須、 は非必須。
•
• ツール
◦
• HDFView (2.8)
http://www.hdfgroup.org/hdf-java-html/hdfview/
◦ Panoply (3.1.4)
http://www.giss.nasa.gov/tools/panoply/download gen.html
• ライブラリ
以下のライブラリは特に断らない限り、ソースを展開したディレクトリにおいて
下記のコマンドを実行してインストール出来る。
./configure --prefix=${HOME}/local
make
make install
• GEOS (3.3.5)
http://download.osgeo.org/geos/geos-3.3.5.tar.bz2
• Szip (2.1)
http://www.hdfgroup.org/ftp/lib-external/szip/2.1/src/szip-2.1.tar.gz
• HDF4 (4.2.8)
http://www.hdfgroup.org/ftp/HDF/HDF Current/src/hdf-4.2.8.tar.gz
HDF4 のコンパイルには Szip ライブラリが必要。
./configure --with-szlib=${HOME}/local --prefix=${HOME}/local
◦ HDF5 (1.8.9)
http://www.hdfgroup.org/ftp/HDF5/current/src/hdf5-1.8.9.tar.gz
HDF5 のコンパイルには Szip ライブラリが必要。
./configure --with-szlib=${HOME}/local --enable-fortran --enable-cxx \
--prefix=${HOME}/local
◦ NetCDF (4.2.1.1)
http://www.unidata.ucar.edu/downloads/netcdf/ftp/netcdf-4.2.1.1.tar.gz
NetCDF のコンパイルには HDF4、HDF5 ライブラリが必要。
CPPFLAGS=-I${HOME}/local/include LDFLAGS=-L${HOME}/local/lib \
./configure --enable-hdf4 --enable-netcdf-4 --prefix=${HOME}/local
• Python (2.7.3)
http://www.python.org/ftp/python/2.7.3/Python-2.7.3.tgz
Python は下記の手順でインストール出来る。
export PYTHONHOME=""
./configure --prefix=${HOME}/local
make
export PYTHONHOME=${HOME}/local
make install
付 録B
62
本演習に必要なソフトウェア
• Python モジュール
以下の Python モジュールは特に断らない限り、ソースを展開したディレクトリにおいて
下記のコマンドを実行してインストール出来る。
•
•
•
•
•
◦
◦
◦
•
•
export PYTHONHOME=${HOME}/local
${HOME}/local/bin/python setup.py build
${HOME}/local/bin/python setup.py install --prefix=${HOME}/local
Numpy (1.6.2)
http://sourceforge.net/projects/numpy/files/latest/download?source=files
Matplotlib (1.1.1)
http://sourceforge.net/projects/matplotlib/files/matplotlib/matplotlib-1.1.1/
matplotlib-1.1.1.tar.gz/download
Natgrid (0.2.1)
http://sourceforge.net/projects/matplotlib/files/matplotlib-toolkits/natgrid-0.2/
natgrid-0.2.1.tar.gz/download
Basemap (1.0.5)
http://sourceforge.net/projects/matplotlib/files/matplotlib-toolkits/basemap-1.0.5/
basemap-1.0.5.tar.gz/download
Basemap のコンパイルには GEOS ライブラリが必要。build の前に以下の環境変数を定義。
export GEOS_DIR=${HOME}/local
PyHDF (0.8.3)
http://sourceforge.net/projects/pysclint/files/latest/download?source=files
PyHDF のコンパイルには HDF4 ライブラリが必要。build の前に以下の環境変数を定義。
export INCLUDE_DIRS=${HOME}/local/include
export LIBRARY_DIRS=${HOME}/local/lib
H5py (2.0.1)
http://code.google.com/p/h5py/downloads/detail?name=h5py-2.0.1.tar.gz&can=2&q=
H5py のインストールには HDF5 ライブラリが必要。build の前に以下の環境変数を定義。
export HDF5_DIR=${HOME}/local
NetCDF4-Python (1.0)
http://netcdf4-python.googlecode.com/files/netCDF4-1.0fix1.tar.gz
NetCDF4-Python のインストールには Szip、HDF5、NetCDF ライブラリが必要。
build の前に以下の環境変数を定義。
export HDF5_DIR=${HOME}/local
export NETCDF4_DIR=${HOME}/local
export SZIP_DIR=${HOME}/local
export USE_NCCONFIG=1
Scipy (0.11.0rc2)
http://sourceforge.net/projects/scipy/files/latest/download?source=files
Scipy のインストールには Numpy が必要。
IPython (0.13)
http://archive.ipython.org/release/0.13/ipython-0.13.tar.gz
その他
ツール: ROOT
ライブラリ: Tcl/Tk、Readline、cURL、SQLite、BLAS、LAPACK、ATLAS
Python モジュール: Cython、Numexpr、ScientificPython、Imaging、ConfigObj、PyMinuit2