Diff of Matlab/Basic2


#br
*matlab講座その2:ROMS outputをmatlabで図化する [#d92af33d]
#br

* ROMS outputファイルを読み込む [#nd8175b4]

例えば,読み込みたいファイルを「roms_his.0124.nc」,1ファイルに12レコード格納されている場合を考える.まず,ncdumpでファイルにどんな変数が格納されているのか見てみる.

 $ ncdump roms_his.0124.nc | less

これでnetcdfファイルの中身を確認する.止める時は「:q」.
このファイルをmatlabで読み込むには,以下のようにする.

 nc=netcdf('roms_his.0124.nc','r');

変数「nc」に読み込むべきファイル(roms_his.0124.nc)が受け渡される.「'r'」はread onlyという意味.ファイルを開いたら,所定の作業が終了した時点できちんと閉じるように.

 close(nc);

この2つのコマンドの間にデータを読み込む.例えば水温tempを読みたい場合は,

 nc=netcdf('roms_his.0124.nc','r');
 dat=nc{'temp'}(:);
 close(nc);

tempは4次元の配列(temp(1:tmax, 1:kmax, 1:jmax, 1:imax))なので,上記のように一度に全てを読み込むと,読み込む時間もかかるし,コンピュータのメモリーを大量に消費してしまうので得策ではない.そういうときはsqueezeコマンドを併用すると,必要な要素だけをピックアップして読み込むことができる.

 tind=3;
 dat=squeeze(nc{'temp'}(tind, :, :, :));

つまり,4次元配列tempのうち,レコード番号3番のみを読むことになる.配列datの大きさは以下のように確認する.

 size(dat)

仮に,imax=320, jmax=256, kmax=32だとすると,sizeコマンドの出力は,

 >> size(dat)
 >> ans = 32, 256, 320

のようになる.応用例として,レコード番号10,表層(k=kmax=32)の水温(つまりSST)を読みたいなら,

 tind=3; kind=32;
 dat=squeeze(nc{'temp'}(tind, kind, :, :););
 dat=squeeze(nc{'temp'}(tind, kind, :, :));
 size(dat)
 ans = 256, 320

となる.

** ROMS gridファイルも読み込んでおく. [#p71d3420]

図化するためには空間座標データも必要になる.手っ取り早いのは,計算に用いたグリッドファイルから座標を得ることである.グリッドファイルを「roms_grd.nc」とすると,

 nc=netcdf('roms_grd.nc','r');
 lon=nc{'lon_rho'}(:);
 lat=nc{'lat_rho'}(:);
 mask=nc{'mask_rho'}(:);
 close(nc);

というような感じで,水平座標lon, latと,マスクを読み込む.図化に際しては,mask=0 (陸グリッド)の部分を表示したくないこともある.そういう場合は,予めmask=0の値を「NaN」に置き換えておくと,陸グリッドを白抜きで図化することができる.

 mask(mask==0)=NaN;

上のコマンドはmatlabっぽい処理の仕方で,配列maskの全要素に対して,maskの値がゼロのものにNaNを代入するというものである.

** 2次元pcolorプロットを作成する. [#d9bd7502]

ncviewで見るようなカラーマップを作成する.上記の手順に従って,座標の値,マスク,SSTが読み込めているとする.

 figure;
 pcolor(lon, lat, dat.*mask);
 hold on;
 shading flat;
 caxis([15 25]);
 colormap(jet);
 colorbar;

これできれいなSSTのマップが緯度経度座標軸を持った2次元プロットとして表示される.
まず,1行目「figure」でグラフを描くべきfigureウィンドウを作成する.以降のプロットはこのfigureウィンドウ内に描画されることになる.2行目はx座標をlon,y座標をlat,データとしてdat.*maskとしたマップを作成するコマンド.

3行目は上記の例では不要だが,1つのfigureウィンドウ内で複数のプロットを作成する際に,今まで描画したプロットを消さないために付けるコマンドである.4行目はpcolorのshadingを制御する関数で、指定なしの場合はグリッド線が描画されてしまう。5行目はデータのレンジ(水温15度〜25度)までをスケーリングして色を付けなさいという意味で,省略するとデータの最小値〜最大値の間でスケーリングが行われる.1つ飛ばして,7行目はカラーバーを付けるコマンドで,6行目はカラーマップのデザインを指定する部分であり,hotとかcoolとか色々なカラースキームがmatlab標準で用意されている.ncviewのカラースキームも使うことができて,その場合は以下のようにする.

 colormap(colormaps_rainbow)

 colormap(colormaps_bright)

また,2行目のpcolor関数の引数において,

 dat*mask

 dat.*mask

の違いに気をつけること.前者はスカラー量同士,あるいはスカラーと行列の掛け算であれば計算できるが,行列同士の乗算には使えない.その場合は後者を使うこと.ポイントは「*」の前にピリオド「.」をいれるということだけである.除算の場合も同様で,行列同士の除算では

 dat./mask
 dat./20

のようにしないとエラーになる.「dat.*mask」の意味するところは,海グリッドではmask=1なので,dat=datのまま,しかし陸グリッドではdat*NaN=NaNに置き換わるということである.

** 相対渦度の2次元pcolorプロットを作成する. [#q33b19a4]

海洋の乱流構造を理解する上で,相対渦度(curl u)を見る必要が出てくる(curlはベクトルuの回転rotationを取る演算子).海洋において渦度というと,通常はポテンシャル渦度を差し,1次オーダーでは相対渦度と惑星渦度の和を水深で割った量,(curl u + f) / h になる(ここでhは水深).このうち,curl uが相対渦度で,惑星渦度 f はコリオリ周波数である.ただ,相対渦度はロスビー数の低い流れ(地球自転効果が卓越する流れ)では,渦度は非常に小さな値を取るので,これをコリオリ周波数 f で無次元化した無次元相対渦度(vortical Rossby numberとも言われる)で表示することが多い.ここではこの無次元相対渦度の可視化を行う.

まずはグリッドファイルから必要な変数を読み込む.先程のSSTより多くの変数が必要になることに注意.

 nc=netcdf('roms_grd.nc','r');
 lon=nc{'lon_rho'}(:);
 lat=nc{'lat_rho'}(:);
 mask=nc{'mask_rho'}(:);
 f = nc{'f'}(:);
 pm = nc{'pm'}(:);
 pn = nc{'pn'}(:);
 close(nc);

pmとpnはグリッド間隔dx,dyの逆数であり,渦度の計算に必要で,fはコリオリ周波数で無次元化に用いる.
続いて,表層における流速の水平2成分をROMS outputファイルから読み込む.

 tind=12; kind=32;
 nc=netcdf('roms_his.0124.nc','r');
 u = squeeze( nc{'u'}(tind,kind, :, :));
 v = squeeze( nc{'v'}(tind,kind, :, :));
 close(nc);

渦度の計算には,関数「vorticity」を使う.

 chi = vorticity(u,v,pm,pn);

これで(有次元:1/s)の相対渦度の鉛直成分(curl u)が求まる.chiはROMSのpsiポイント(i方向に-1/2グリッド,j方向に-1/2グリッドだけrhoポイント(グリッド中央)から離れた点)で定義されているので,rhoポイントに戻す.

 chi = psi2rho(chi);

さらに,fで無次元化するので,

 chi = chi ./ f;

とする.ROMSがメソスケール渦を解像する程度の格子間隔(5km以上程度)であれば,この無次元化されたchiは -1 < chi < 1程度のレンジに入ることが多く,解像度を上げると,chiの変化も大きくなる.
プロットはSSTと同様に以下のようにすれば良い.

 figure;
 pcolor(lon, lat, chi.*mask);
 hold on;
 shading flat;
 caxis([-1 1]);
 colormap(colormaps_rainbow);
 xlabel('lon');
 ylabel('lat');
 title('normalized relative vorticity');
 colorbar;

SSTと違ってchiは正負の値を取るので,カラーマップはjetやbrightよりもrainbowの方が好ましい.

** プロットを画像ファイルとして保存する. [#y8df42cd]

このようにして作成した図面は,簡易的に表示するだけなら「jpg」や「png」などのラスター形式(表示や印刷の汎用性は高いが,拡大するとシャギーが目立つ)で,論文や報告書,印刷物に載せることまで考えるなら「eps」や「pdf」などのベクトル形式(汎用性はやや低いが,拡大しても品質が落ちない)で保存することになる.手順は,

 figure;
 plot(x,y,dat);
 print('-djpeg90', 'output.jpg');

のように,直近で作成したfigureウィンドウに表示されているものをファイルに落とす(dampする)ことができる.上の例では,jpegフォーマット,圧縮率90%で,出力ファイル名はoutput.jpgを指定している.

 print('-dpng', '-r200', 'output.png');  % PNGフォーマット,解像度200(大きいほど高解像度)
 print('-depsc2', 'output.eps');        % EPSフォーマット,カラー,レベル2
 print('-dpdf', 'output.pdf');                % PDFフォーマット

など.印刷品質を考えるとEPS出力が好ましく,印刷物上でカラーをより綺麗に見せるには

 print('-depsc2', '-cmyk', 'output.eps'); 

とすると良い.

PDF出力は汎用性が高くてファイルサイズもEPSよりは小さくなるので使いたいところだが,A4などの用紙サイズにあわせてしまうのか,figureの周りに余白が入ってしまう.今のところ,この余白を簡単に除去する方法が分からないので,実に使いづらい.誰か解決方法を知っていたら教えて下さい.