Tips: NetCDF

NetCDFライブラリのインストール
- Fortran
- Matlab
ファイル読み込みの初歩
GrADSで読めるNetCDFファイルの新規作成
- Matlab
- Matlab (mexncの場合)

NetCDFライブラリのインストール

NetCDF形式のデータを扱うためには，NetCDFのインストールが必要です．ここではFortranとMatlabの例を紹介します．

Fortran用のNetCDFのインストール

Fortran用のNetCDFを，環境：Vine 4.1 & Intel fortran 8.1にインストールした例を紹介します．

NetCDFのwebサイト http://www.unidata.ucar.edu/software/netcdf/で， netcdf.tar.Z (ver. 3.6.0)をダウンロードする．
展開する．(以下の作業はrootで行う)
```
# tar -xzvf netcdf.tar.Z
```
できたディレクトリに移動する．
```
# cd netcdf-3.6.0-p1/src/
```
環境変数の設定をする．ここでは http://www.unidata.ucar.edu/software/netcdf/docs/other-builds.htmlの「Linux x86 with Intel ifort 8.1 fortran compiler」を参考にした．
```
# export FC=ifort
# export F90=ifort
# export FFLAGS="-O -mp"
# export CPPFLAGS="-DNDEBUG -DpgiFortran"
```
configureを実行する．/usr/local/netcdf-3.6.0はインストールする場所である．
```
# ./configure --prefix=/usr/local/netcdf-3.6.0
```

makeする．

# make
# mkdir /usr/local/netcdf-3.6.0
# make install

再び，環境変数の設定をする．/etc/profileに以下を書き込む．

PATH="$PATH:/usr/local/netcdf-3.6.0/bin/"
LD_LIBRARY_PATH="$LD_LIBRARY_PATH:/usr/local/netcdf-3.6.0/lib"

Matlab用のNetCDFのインストール

Matlab用のNetCDFであるMEXNCをインストールした例を紹介します．

SOURCEFORGE.NETのサイト http://mexcdf.sourceforge.net/index.htmlで， Matlabのバージョンに合ったMEXNCをダウンロードする．

好みのディレクトリで展開する．

$ mkdir ~/Matlab
$ cd ~/Matlab
$ tar -xzvf mexnc-2.0.30.tar.gz

Matlabを起動して，展開したディレクトリへパスを通す．
```
>> addpath ~/Matlab/mexnc/
>> savepath
```

ファイル読み込みの初歩

NetCDFファイルをFortranとMatlabで読み込んだ例を紹介します．ここでは AvisoのFTPサイトにある，海表面高度偏差データファイル dt_ref_global_merged_msla_h_19921014_19921014_20051108.nc.gzを読み込みます．

変数情報の取得

FortranとMatlabに共通する手順として，NetCDFファイルの変数情報の取得があります．

まずファイルを解凍します．

$ gunzip dt_ref_global_merged_msla_h_19921014_19921014_20051108.nc.gz

次に，ファイルについての情報を得るために，

$ ncdump -c dt_ref_global_merged_msla_h_19921014_19921014_20051108.nc | less

とします．すると出力の前半の内容は以下の通りになっています．

netcdf dt_ref_global_merged_msla_h_19921014_19921014_20051108.nc {
dimensions:
        NbLatitudes = 915 ;
       NbLongitudes = 1080 ;
        GridDepth = 1 ;
        LatLon = 2 ;
        
variables:
        int LatLon(LatLon) ;
                LatLon:_FillValue = 2147483647 ;
                LatLon:long_name = "No sense but necessary for some automatic tools" ;
                LatLon:units = "count" ;
        double NbLatitudes(NbLatitudes) ;
                NbLatitudes:_FillValue = 1.84467440737096e+19 ;
                NbLatitudes:long_name = "Latitudes" ;
                NbLatitudes:units = "degrees_north" ;
        double NbLongitudes(NbLongitudes) ;
                NbLongitudes:_FillValue = 1.84467440737096e+19 ;
                NbLongitudes:long_name = "Longitudes" ;
                NbLongitudes:units = "degrees_east" ;
        double LatLonMin(LatLon) ;
                LatLonMin:_FillValue = 1.844674407370955e+19 ;
                LatLonMin:long_name = "Latitude/Longitude of south/ouest corner" ;
                LatLonMin:units = "degree" ;
        double LatLonStep(LatLon) ;
                LatLonStep:_FillValue = 1.844674407370955e+19 ;
                LatLonStep:long_name = "latitude/longitude steps" ;
                LatLonStep:units = "degree" ;
        float Grid_0001(NbLongitudes, NbLatitudes) ;
              Grid_0001:_FillValue = 1.844674e+19f ;
              Grid_0001:long_name = "H" ;
              Grid_0001:units = "cm" ;
              Grid_0001:data = "1992-10-14 00:00:00.000000 UTC"

海表面高度偏差を得るのに必要な情報は，色で強調した2箇所です．

まず，上の部分から，緯度方向のグリッド数(NbLatitudes)が915，経度方向のグリッド数(NbLongitudes)が1080であることが判ります．

次に下の部分から，海表面高度偏差の値はデータの型が単精度(float)で，変数名が「Grid_0001」でグリッド数が(1080×915)，欠損値が1.844674e+19f，正確な変数名は「H」，単位が「cm」であること等が判ります．

FortranでのNetCDFファイル読み込み

上述のデータをFortranで読み込みます．

最小限のNetCDFライブラリ

Fortranの読み込みプログラム例netcdf_example.f90を示します．必要最小限のコマンドのみ実行しています．重要なことはinclude文によるnetcdf.incの読み込みと， NetCDFライブラリの使い方です．

program netcdf_example
  implicit none
  include '/usr/local/netcdf-3.6.0/include/netcdf.inc'
!
  integer, parameter :: xsz=1080, ysz=915  ! 次元長
  integer :: ncid, status, rhid
  real :: sla(ysz,xsz) ! データ
  character :: filename*60, var*10

  var='Grid_0001' ! slaの変数設定名
  filename='dt_ref_global_merged_msla_h_19921014_19921014_20051108.nc'

  status = nf_open(filename, nf_nowrite, ncid) ! ファイルのopenとNetCDF ID(ncid)の取得

  status = nf_inq_varid(ncid, var, rhid) ! 変数のID(rhid)取得

  status = nf_get_var_real(ncid, rhid, sla)! 変数読み込み

  status = nf_close(ncid) ! ファイルのclose

  stop
end program

このプログラム中で使用されているNetCDFライブラリは次の4つです．

nf_open
ファイルをopenします．2つ目の引数nf_nowriteは読み込み専用を意味します．
nf_inq_varid
読み込む変数のID(rhid)を取得します． 2つ目の引数varには，先ほどのncdumpで得た変数名Grid_0001を入れます．
nf_get_var_real
実際に変数の値をすべてを読み込みます．この場合，配列slaに海表面高度偏差の値が格納されます．先ほどのncdumpで得た情報から，「Grid_0001」の変数の型は単精度(float)であったので， nf_get_var_realを使用しています．例えば倍精度(double)であった場合には， nf_get_var_doubleを使用します．
nf_close
ファイルをcloseします．

NetCDFライブラリについては netCDFユーザガイド日本語版＠地球流体電脳倶楽部が詳しいです．

このプログラムをコンパイルするには，

$ ifort netcdf_example.f90 -L/usr/local/netcdf-3.6.0/lib -lnetcdf

とします．指定するディレクトリはインストール場所により異なります．

MatlabでのNetCDFファイル読み込み

上述のデータをMatlabで読み込みます．

Matlab付属のNetCDFライブラリ

最近のMatlabではmexncを使用しなくても，付属のNetCDFライブラリでデータの読み書きができます．

上述のデータをMatlab付属のNetCDFライブラリで読み込むと以下のようになります．

filename='dt_ref_global_merged_msla_h_19921014_19921014_20051108.nc';

var='Grid_0001'; % slaの変数設定名

[sla]=ncread(filename,var); % 変数読み込み

mexncライブラリを用いた読み込み

mexncを用いたMatlabでの読み込みプログラム例netcdf_example.mを示します．必要最小限のコマンドのみ実行しています．重要なことはmexncコマンドの使い方です．

filename='dt_ref_global_merged_msla_h_19921014_19921014_20051108.nc';

[ncid,status] = mexnc('OPEN',filename,nc_nowrite_mode); % fileのopenとNetCDF ID(ncid)の取得
if(status~=0); error('file not open'); end

var='Grid_0001'; % slaの変数設定名
[varid,status] = mexnc('INQ_VARID',ncid,var); % 変数のID(varid)取得

[sla,status] = mexnc('GET_VAR_FLOAT', ncid,varid); % 変数読み込み

[status] = mexnc('CLOSE',ncid); % fileのclose

このプログラム中でmexncコマンドは4回使用されています．最初の引数の違いにより，それぞれ次の動作を行います．

mexnc('OPEN',filename,nc_nowrite_mode)
ファイルをopenします．3つ目の引数nf_nowrite_modeは読み込み専用を意味します．
mexnc('INQ_VARID',ncid,var)
読み込む変数のID(rhid)を取得します． 3つ目の引数varには，先ほどのncdumpで得た変数名Grid_0001を入れます．
mexnc('GET_VAR_FLOAT', ncid,varid)
実際に変数の値をすべてを読み込みます．この場合，配列slaに海表面高度偏差の値が格納されます．先ほどのncdumpで得た情報から，「Grid_0001」の変数の型は単精度(float)であったので， 'GET_VAR_FLOAT'を使用しています．例えば倍精度(double)であった場合には， 'GET_VAR_DOUBLE'を使用します．
mexnc('CLOSE',ncid)
ファイルをcloseします．

mexncコマンドがどのような引数をとるのかについてのより詳しい情報は，Matlab上で

>> help mexnc-doc

とすることで表示されます．

GrADSで読めるNetCDFファイルの新規作成

Matlab

ここではMatlabを用いて，GrADSのsdfopenで読めるNetCDFファイル作成の例を紹介します．

Matlab付属のNetCDFライブラリを用いた場合のプログラム例grads_nc.mを示します．

lat=[10:5:50];  ysz=length(lat); % 緯度と経度は適当です
lon=[100:5:200];  xsz=length(lon);
tsz=5; % 時間ステップは5とします
data=randn(xsz,ysz,tsz); % このランダムデータを出力します．
filename='test.nc'; % 出力するNetCDFファイル名です．

% ここからNetCDFライブラリを使用します
nw = netcdf.create(filename, '64BIT_OFFSET'); % ファイル作成

dimid1=netcdf.defDim(nw,'longitude',xsz); % 次元の作成
dimid2=netcdf.defDim(nw,'latitude',ysz);
dimid3=netcdf.defDim(nw,'time',tsz);

varid=netcdf.defVar(nw,'longitude','double',[dimid1]); % 変数とその属性の作成
netcdf.putAtt(nw,varid,'long_name','longitude');
netcdf.putAtt(nw,varid,'units','degrees_east');

varid=netcdf.defVar(nw,'latitude','double',[dimid2]);
netcdf.putAtt(nw,varid,'long_name','latitude');
netcdf.putAtt(nw,varid,'units','degrees_north');

varid=netcdf.defVar(nw,'time','double',[dimid3]);
netcdf.putAtt(nw,varid,'long_name','time');
netcdf.putAtt(nw,varid,'units','days since 0000-01-01'); % dailyデータにします

varid=netcdf.defVar(nw,'temperature','double',[dimid1 dimid2 dimid3]);
netcdf.putAtt(nw,varid,'long_name','sea surface temperature');
netcdf.putAtt(nw,varid,'units','Kelvin');

varid2=netcdf.getConstant('NC_GLOBAL'); % グローバル属性は任意です．
netcdf.putAtt(nw,varid2,'title','sea surface temperature data');
netcdf.putAtt(nw,varid2,'date',date);

netcdf.close(nw); % ファイルを閉じます

ncwrite(filename,'longitude',lon); % 実際の値の書き込み
ncwrite(filename,'latitude',lat);
ncwrite(filename,'time',[0:4]);
ncwrite(filename,'temperature',data);

これでgrads上で以下のようにすれば．ファイルが開けるはずです．

ga-> sdfopen test.nc

Matlab (mexncの場合)

mexncライブラリを用いた場合のプログラム例grads_nc.mを示します．

lat=[10:5:50];  ysz=length(lat); % 緯度と経度は適当です
lon=[100:5:200];  xsz=length(lon);
tsz=5; % 時間ステップは5とします
data=randn(xsz,ysz,tsz); % このランダムデータを出力します．
filename='test.nc'; % 出力するNetCDFファイル名です．

% ここからNetCDFライブラリを使用します
[ncid,status] = mexnc('CREATE',filename,nc_clobber_mode); % ファイル作成

[dimid,status] = mexnc('DEF_DIM',ncid,'longitude',xsz);
[dimid2,status] = mexnc('DEF_DIM',ncid,'latitude',ysz);
[dimid3,status] = mexnc('DEF_DIM',ncid,'time',tsz);

[varid,status] = mexnc('DEF_VAR',ncid,'longitude','double',1,dimid);
status = mexnc('PUT_ATT_TEXT',ncid,varid,'units','char',... 
 length('degrees_east'),'degrees_east');

[varid2,status] = mexnc('DEF_VAR',ncid,'latitude','double',1,dimid2);
status = mexnc('PUT_ATT_TEXT',ncid,varid2,'units','char',... 
 length('degrees_north'),'degrees_north');

[varid3,status] = mexnc('DEF_VAR',ncid,'temperature','double',3,[dimid dimid2 dimid3]);
status = mexnc('PUT_ATT_TEXT',ncid,varid3,'units','char',... 
 length('Kelvin'),'degree Celsius');
status = mexnc('PUT_ATT_TEXT',ncid,varid3,'long_name','char',... 
 length('temperature'),'temperature');

[varid4,status] = mexnc('DEF_VAR',ncid,'time',xtype,1,dimid3);
status = mexnc('PUT_ATT_TEXT',ncid,varid4,'units','char',... 
 length('days since 0000-01-01'),'days since 0000-01-01'); % dailyデータにします

status = mexnc('ENDDEF',ncid);

[varid,status] = mexnc('INQ_VARID',ncid,'longitude'); 
status = mexnc('PUT_VAR_DOUBLE',ncid,varid,lon);

[varid2,status] = mexnc('INQ_VARID',ncid,'latitude');
status = mexnc('PUT_VAR_DOUBLE',ncid,varid2,lat);

[varid3,status] = mexnc('INQ_VARID',ncid,'temperature');
status = mexnc('PUT_VAR_DOUBLE',ncid,varid3,data);

[varid4,status] = mexnc('INQ_VARID',ncid,'time');
status = mexnc('PUT_VAR_DOUBLE',ncid,varid4,[0:4]);

status = mexnc('CLOSE',ncid); % ファイルを閉じます

これでgrads上で以下のようにすれば．ファイルが開けるはずです．

ga-> sdfopen test.nc

Homeに戻る