MatlabでGrADSデータを読み書きする方法

GrADSデータの構造
GrADSデータの読み込み
GrADSデータの書き込み

GrADSデータの構造

GrADSのデータファイル(拡張子はgrdが一般的)は，4バイト実数(単精度)のバイナリ―データの羅列であり，必ずctlファイルとセットである．

grdファイルにデータは， x (経度) → y (緯度) → z (高度) → v (変数)→ t (時間) の順に入っている．

以下，説明の簡略化のため，次のように定義する．

xsz: 経度方向のグリッド数
ysz: 緯度方向のグリッド数
zsz: 高さ方向のグリッド数
vsz: 変数の数
tsz: 時間ステップ数

このときgrdファイルのサイズは， 4 byte×xsz×ysz×zsz×vsz×tszとなる．

3次元データの例

xsz=4, ysz=3, tsz=2の3次元データの場合，以下の順番でgrdファイルにデータが入っている．

3次元データの順番

4次元データの例

xsz=4, ysz=3, vsz=2, tsz=2の4次元データの場合は，以下の順番でgrdファイルにデータが入っている．

4次元データの順番

GrADSデータの読み込み

GrADSのgrdファイルは，4バイト実数(float32)で読みこめばよい．以下は，xytの3次元データtest.grdを読み込むプログラムです．

fid=fopen('test.grd');
data=fread(fid,'float32');
fclose(fid);
data=reshape(data,xsz,ysz,tsz);

これを関数(rdgrd.m)にすると以下になります．

function data=rdgrd(filename)
  fid=fopen(filename);
  data=fread(fid,'float32');
  fclose(fid);
return

GrADSデータの書き込み

GrADSのgrdファイルは，4バイト実数で出力すれば良い． (ctlファイルは別途用意する必要がある．) 以下は，xytの3次元データtest.grdを書き出すプログラムである．

data=reshape(data,xsz,ysz,tsz);
fid=fopen('test.grd','wb'); % b denotes binary
fwrite(fid,data,'float32');
fclose(fid);

これを関数(wrtgrd.m)にすると以下になります．

function wrtgrd(data,filename)
  fid=fopen(filename,'wb');
  fwrite(fid,data,'float32');
  fclose(fid);
return

さらにこの関数にエラー処理等をつけて， wrtgrd.mは完成です．

function [count]=wrtgrd(data,outname,flag)
% function [count]=wrtgrd(data,outname,flag)
% flag==1 write (default)
% flag~=1 append

if(nargin<3 || flag==1)
  fid=fopen(outname,'wb');
else
  fid=fopen(outname,'ab');
end
if(fid<=0); error('disable to open'); end

szdm=size(data);
count=fwrite(fid,data,'float32');
fclose(fid);
if count~=prod(szdm); error('disable to write'); end

Homeに戻る