Backup of Matlab/Basic2(No. 3)

 

matlab講座その2：ROMS outputをmatlabで図化する †

 
↑

ROMS outputファイルを読み込む †

例えば，読み込みたいファイルを「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, :, :););
size(dat)
ans = 256, 320

となる．

↑

ROMS gridファイルも読み込んでおく． †

図化するためには空間座標データも必要になる．手っ取り早いのは，計算に用いたグリッドファイルから座標を得ることである．グリッドファイルを「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プロットを作成する． †

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プロットを作成する． †

海洋の乱流構造を理解する上で，相対渦度（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の方が好ましい．

↑

プロットを画像ファイルとして保存する． †

このようにして作成した図面は，簡易的に表示するだけなら「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フォーマット，カラー，レベル２
print('-dpdf', 'output.pdf');                % PDFフォーマット

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

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

とすると良い．

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