[1.0] はじめに
前回，みなさんには地球惑星科学分野における「知の集積」の役割を担う数値モデルの一つのDCPAM を動かしてもらいました．
ここではその計算データを Dennou Ruby Project で開発されているツール群を用いて解析し，描画することで(目で見て)分かりやすい形にしてみましょう．ちなみに，人の目で見て分かりやすい形にすることを「可視化」と言います．
[1.0.1] 事前準備
今回も実技時間内での作業は 1 つのアカウントで行い，適宜入力する人を交代して進めてください．他のアカウントでも同じ作業を行ないたい場合には授業時間外にやってください.
まず，インストールされているパッケージに更新が無いか確認を行ってください．
(この作業の目的は，第 8 回 「0. 今後, 情報実習を受講するにあたってのお願い 」を参照)
$ sudo apt-get update |
$ sudo apt-get upgrade |
[1.1] Dennou Ruby Project のツールとは
Dennou Ruby Project のツールはオブジェクト指向スクリプト言語の Ruby を使って作られています.
Ruby は日本人のまつもとひろゆき氏によって開発されており,日本語の書籍も多く出版されています. 詳しくはRuby のホームページなどを参照してください. Ruby はオブジェクト指向に基づいて設計された言語です. 非常に粗っぽい説明をすれば, プログラムの世界におけるオブジェクトとは,「データ」と手続きを一体にしたようなものです. それぞれのデータのタイプに応じて手続きを定めるというのを基本にしてプログラミングを行ないます. なお, ここまで「手続き」と呼んでいたものは Ruby の世界ではメソッドと呼ばれます.
ruby を使ったスクリプトの例としてシェルスクリプトの課題 (INEX 第 3 回, 04/28) の解答となる ruby スクリプトを挙げます. こちらを見て各自で作成したシェルスクリプトと比較してみてください. だいぶ様子が違うのではないでしょうか.
今回の実習で使う Dennou Ruby Project のツールのうち, 中心となるのはGPhys というものです. GPhys については GPhys のホームページを参照してください. ここにある GPhys のチュートリアルを一通りやると GPhys の基礎を学ぶことができます.
ここでは, GPhys のチュートリアルの中の以下の部分を実際にやって, GPhys と Ruby に慣れることにしましょう.
- 「1. はじめる前に 」の中の「irb について」, 「サンプルデータのダウンロード」, 「netCDF データの内容確認」
- 「2. 可視化の第一歩」の中の「最速描画」, 「いろんな種類の図を描いてみる」, 「描画オプション (1)」
- 「3. 可視化の第二歩」の中の「いろんな断面」
[1.2] 計算結果の netCDF ファイルの中身を確認する
解析・描画の作業を開始する前に, 使用する netCDF ファイルの中身を見ておくことにします. まず実験を行なったディレクトリに行ってください. そして, 以下のようにしてファイルの中身を見てください.
$ cd ~/dcpam/dcpam5-exp/Earth-exp/ $ ncdump Temp.nc | less |
まずデータの軸変数を見ておきましょう. Temp.nc の variables: の項目を見てください. 軸に関する変数の中には以下のものがあります.
- lon : 経度
- lat : 緯度 DCPAM はスペクトル法で計算を行なっていますが, データファイルには格子点値に変換されたものが出力されます. 格子点の間隔はおよそ何km になっているかわかりますか?
- sig : σ=p/ps (ps は地表面気圧). DCPAM の鉛直座標は気圧に基づいた ものになっています. 対流圏界面は σ=0.15 付近にあります. 対流圏 (高さおよそ 10km) の 中に鉛直格子点はいくつあるでしょうか?
- time : 時刻
次に, time の値を確認してください.
time = 0, ......, 1095 |
[1.3] 作図スクリプトを作成して大気の状態を眺めてみよう
前回, irb を使って図を 1 枚かきました. 今回はスクリプトを使って解析描画を行なってみましょう. zonalmean_3dim_snapshot.rb は前回 irb を使って行なったものと同様の内容をスクリプトに保存したものです. まず, このスクリプトをダウンロードして実行できるか確認しましょう.- まず zonalmean_3dim_snapshot.rb をダウンロードする.
$ wget http://www.ep.sci.hokudai.ac.jp/~inex/y2017/0721/practical/zonalmean_3dim_snapshot.rb
- zonalmean_3dim_snapshot.rb を実行できるようにパーミッションの設定を行ないましょう. パーミッションについて忘れた方は( 最低限 Unix / Linux [I] 4. パーミッション) を参照してください.
$ chmod 755 zonalmean_3dim_snapshot.rb
- 描画を行ってみる.
$ ./zonalmean_3dim_snapshot.rb
図1: 2 日目の気温の経度平均平均分布 - 終了する場合は，"q" を入力
$ less zonalmean_3dim_snapshot.rb |
いまは 2 日目の瞬間値の図 (スナップショットと言います) を書いてみましたが,これからは zonalmean_3dim_snapshot.rb を編集して, 時間を変更した図も作ってみましょう.
以下の三つについて図を作ってください．
- [図2] データに格納されている最後の時刻 (1095 日目) のスナップショット(時間平均していないもの)
- [図3] 最後の 1 年 (731 日目 ~ 1095 日目) の時間平均
- [図4] 最後の年の北半球が夏の間 (881 日目 ~ 971 日目) だけの時間平均
gphys = GPhys::IO.open('Temp.nc', 'Temp').cut('time'=>0..1095)
:
GGraph.tone( gphys.mean('lon').cut('time'=>731..1095).mean('time') )
|
図2: 1095 日目の温度の経度平均分布
図1 の 2 日目の図とは異なり温度分布が形成されていることがわるでしょう. 地表付近では, 低緯度で高温・高緯度で高温となっていることがわかると思います. 赤道と極の温度差はおよそ 70K になっていますね. また, 熱帯域では地表から上空に向かって高度を高くしていくとσ =0.1 付近 (気圧にしておよそ 100hPa, 高度にしておよそ 15km) で 温度が極小になっています. これが対流圏界面と呼ばれる領域です.
図3: 3 年目の 1 年平均をとった温度の経度平均分布
1 年平均をしたために, 図2 に比べて南北対称なパターンになっていますね.
図4: 3 年目の夏季平均をとった温度の経度平均分布
図4 では, 図3 と比べて南極の上空における低温領域が目立ちます. この場合で南極上空で温度が下がる理由はわかりますか?
今度は, 違う変数の図を書いてみましょう. 東西風のデータファイル U.nc (東西風の変数名: U )を使って次の図を描いてみてください.
- [図5] データに格納されている最後の時刻のスナップショット (時間平均していないもの)
- [図6] 3 年目 (731 日目から 1095 日目まで) の 1 年平均
- [図7] 最後の年の夏の間だけの時間平均
図5: 1095 日目の東西風の経度平均分布
図5 では, 緯度 40 度, σ =0.2 を中心とした領域で東西風が正 (西風. 西から東に向かって吹く風) で強くなっています. これがジェット気流と呼ばれるものです.
図6: 3 年目の1 年平均をとった東西風の経度平均分布
図7: 3 年目の夏季平均をとった東西風の経度平均分布
なお, 作成した図をファイルに保存するには, スクリプト中において
DCL.gropn(2) |
$ convert -rotate 90 dcl.pdf Temp_zonalmean.png |
また, dcl.pdf は上書きされてしまうので画像を保存しておきたい場合には, 毎回 convert してファイル名を変えるようにしてください.
[1.4] 水に関する量を作図してみよう
降水量と蒸発量のデータを使って簡単なデータ処理をやってみましょう.
まず, horizontalmap_2dim_snapshot.rb をダウンロードして実行してみてください. このスクリプトは 3 年目の最後の瞬間の降水量の経度緯度分布の図(図8)を描くものです.
$ wget http://www.ep.sci.hokudai.ac.jp/~inex/y2017/0721/practical/horizontalmap_2dim_snapshot.rb |
図8: 1095 日目の降水量の経度緯度分布
図8 は緯度経度分布を示す図になっているので, 海岸線もプロットされています.おおまかには, 降水量は赤道で大きくなっています. これは赤道上で上昇流が生じ雲が形成されるためです. また, 日本付近でも降水量が大きくなっています. これは低気圧に伴う雨です. 図8 はとある瞬間の図に過ぎませんが, 以下で作成する時間平均図でも同様の傾向が現れているのかよく見てみてください.
また降水量に関する別の図も作ってみましょう. 以下の平面分布の図を作ってください.
- [図9] 3 年目の 1 年平均
- [図10] 最後の年の夏の間だけの時間平均
図9: 3 年目の1 年平均をとった降水量の経度緯度分布
図10: 3 年目の夏季平均をとった降水量の経度緯度分布
次に，降水量の東西平均 (経度方向の平均) の図も作ってみましょう. そのためのスクリプト zonalmean_2dim_snapshot.rb をダウンロードして実行してください.
図11: 1095 日目の降水量の経度平均分布
- [図12] 3 年目の 1 年平均
- [図13] 最後の年の夏の間だけの時間平均
図12: 3 年目の1 年平均をとった降水量の経度平均分布
上で描いた経度緯度分布 (図9) では海陸分布の影響が強く出ていますが, 図12 ではその影響を目立たなくなります. その分, 降水量の緯度変化のおおまかな様子が よくわかるようになります. 赤道付近と中緯度域 (緯度にして 40 度から 50 度)の領域で降水量が大きくなるということがわかると思います.
図13: 3 年目の夏季平均をとった降水量の経度平均分布
更に，緯度方向にも平均をして降水量の水平平均 (経度方向にも緯度方向にも平均) の図を作ってみましょう. そのためのスクリプト globalmean_2dim_timeseries.rb をダウンロードして実行してください.
図14: 3 年目の水平平均をとった降水量の時間変化
では, ここまでに用いたスクリプトを改変して, 表面からの蒸発量の図を描いてみましょう. 蒸発量のデータファイルは SurfH2OVapFlux.nc です. このファイルの中には SurfH2OVapFlux という変数名で蒸発量が格納されています. 蒸発量に関して, 以下の図を何枚か作ってみてください.
- [図15] 3 年目の 1 年平均の経度緯度分布
- [図16] 最後の年の夏の間だけの時間平均の経度緯度分布
- [図17] 3 年目の最後の瞬間の経度緯度分布
- [図18] 3 年目の 1 年平均の経度平均分布
- [図19] 最後の年の夏の間だけの時間平均の経度平均分布
- [図20] 3 年目の 1 年間の水平平均の時間変化
図15: 3 年目の 1 年平均をとった蒸発量の経度緯度分布
降水の経度緯度分布(図9) とは異なり熱帯の海洋上で蒸発量が大きくなっていることがわかります. おおざっぱには, 蒸発量は表面温度が高いほど大きくなるためです.
図16: 3 年目の夏季平均をとった蒸発量の経度緯度分布
図17: 1095 日目の蒸発量の経度緯度分布
図18: 3 年目の1 年平均をとった蒸発量の経度平均分布
図19: 3 年目の夏季平均をとった蒸発量の経度平均分布
図20: 3 年目の水平平均をとった蒸発量の時間変化
最後に, 降水と蒸発の差 (ここでは P-E と表します) を計算してそれを図にしてみましょう. ここでは, 水平分布の図を作ってみます. horizontalmap_2dim_snapshot.rb を更に改変しましょう. 降水量の読みこみ先を gphys_prcp に変更した上で, 蒸発量をの読みこみ先を gphys_h2ovapflux とする設定を追加しましょう. 更に, スクリプト内で gphys_prcp - gphys_h2ovapflux を計算して, 3 年目の最後の瞬間の P-E の平面分布の図を描けるようにしてください.
図21: 1095 日目の降水と蒸発の差の経度緯度分布
[1.5] 課題に向けての準備
今日はレポート課題が出題されます. 内容は P-E に関する以下の図を作成し, それらの図に関する記述を行なうというものです.
- 3 年目の 1 年平均の経度平均分布
- 3 年目の 1 年間の水平平均の時間変化