| はじめに：netcdfの基本的な情報の取得 |
まず、netcdfファイルの中身がどのようなものかを調べるところから始めます。
ncdump -h xxx.nc
として、中に書き込まれている変数と軸情報を取得します。たとえば気象庁MSMをもとに作られたnetcdfファイルだと、
ncdump -h MSM2013102421P.nc の結果 netcdf MSM2013102421P { dimensions: lon = 241 ; lat = 253 ; p = 16 ; time = 12 ; variables: とでてきます。これは、lon軸に241個、lat軸に253個、p軸に16個、time軸に12個の長さを持っていることを示しています。さらに下の方を読んでいくと、 variables: (中略) float z(time, p, lat, lon) ; z:long_name = "geopotential height" ; z:units = "m" ; z:standard_name = "geopotential_height" ; と書いてあります。これは、変数[z](ジオポテンシャル高度)はtime,p,lat,lonの4つの軸で定義されていることを示しています。 ncdump -c xxx.nc とすると軸情報が詳しく出てきます。たとえばこの場合、lon, lat, timeは等間隔ですが、pは不当感覚であることがわかります。 |
| ncviewを使ったクイックルック | クイックルックならncviewが便利です。 ncview xxx.nc とすればGUIのインターフェースが立ち上がり、変数情報をクリックするとすぐに図を書いてくれます。見やすくするのは難しいので、クイックルック専用で使っています。。。 |
| GrADSで読む |
通常は、gradsを起動した後、.ncファイルを引数にsdfopenすればよい。 しかし、COARDS規約でない.ncファイルはsdfopenできない（エラーメッセージSDF file has no discernable X coordinate. To open this file with GrADS, use a descriptor file with an XDEF entry.が出る）ので、ctlファイルを作成します。上記のMSMのnetcdfファイルの場合、この場合つくるべきctlファイルは、 DSET ^MSM2013102421P.nc DTYPE netcdf TITLE MSM pressure level UNDEF 0.0 XDEF 241 linear 120.0 0.125 YDEF 253 linear 22.4 0.1 ZDEF 16 levels 1000 975 950 925 900 850 800 700 600 500 400 300 250 200 150 100 TDEF 12 21Z24OCT2013 180mn VARS 2 z=>z 16 t,z,y,x geopotential height w=>w 16 t,z,y,x wind vertical velocity ENDVARS のようなものになります。 DSETには、読むべきnetcdfのファイルを指定し、^をつけると相対パスで指定できます。 DTYPEには、netcdfです。 UNDEFには、欠損値の数値を入れます。ncdump -cしたときにUNDEFが書かれていればその数値を指定します。 XDEFにはX方向の軸を指定します。ncdump -cで出てくるx軸の情報に従って書き込みます。ZDEFのように値を241個並べてもよいのですが、格子間隔が線型である場合にはこのように簡略化して書けます。最初の241は配列長、linearは1格子の長さが等間隔であることを指定しています。120.0はx=1の位置、0.125は1格子の長さです。 YDEFには、XDEFと同様に、Y方向の軸を指定します。 このファイルの場合Z軸が等間隔でないので、Z軸を直接書き込みます。最初の16は配列長、levelsは1格子の長さが等間隔でないことを示します。 TDEFにはT方向の軸を指定します。21Z24OCT2013はt=1の時刻、180mnは1あたりの時間（分）です。mn（分）以外にhr（時）,yr（年）などが指定できます。 VARS 2はこのctlファイルで扱う変数の数を2個に指定する、という意味です。 z=>z がnetcdfに書かれている変数情報をGrADSに渡すための情報です。矢印の出所がGrADSで扱う変数名、向いた先がnetcdfにかかれているデータの変数です。したがってgh=>zとすると、GrADSではghという変数で字おポテンシャル高度を扱うことになります。次の16はz方向の配列長、t,z,y,xは配列の次元です。この順番で並べる必要があります。書きたい変数1個につき、1行書き入れます。 最後にENDVARSをいれます。あとはこのctlファイルをopenします。 |
| GMTで読む | シンプルにやるには、grd2xyzを使います。varを変数名, timeをnetcdfの時刻番号（最初の時刻が0であることに注意。）として grd2xyz xxx.nc?var[time] をすると(x座標, y座標, 値)のテーブルが出力されます。 |
| fortranで読む |
最小限のプログラムは、 program netcdf_example implicit none include '/usr/include/netcdf.inc' integer, parameter :: nx=241, ny=253, nz=16 , nt=12 ! 次元長 integer :: ncid, status, rhid real(4) :: hz(nx,ny,nz,nt) ! データ character(30) :: filename, var integer :: i,j,k,l var='z' ! 変数設定名。''にはnetcdfに入っている変数名を書く。 filename='MSM2013102421P.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, hz)! 変数読み込み。3個目の引数には自分で変数名を設定 status = nf_close(ncid) ! ファイルのclose ! do l=1,nt ! do k=1,nz ! do j=1,ny ! do i=1,nx ! write(*,*) i,j,k,l,hz(i,j,k,l) ! end do ! end do ! end do ! end do stop end program です。ポイントは、3行目のinclude文、途中のnf_*のライブラリの設定です。 コンパイルする時には gfortran readnetcdf.f90 -L/usr/local/netcdf-4.3.0/lib -lnetcdf -lnetcdff としてnetcdfのライブラリの場所を-Lで、ライブラリ名を-lで指定します。やりやすいミスとして、fortranで定義する配列hz(xsz,ysz,zsz,tsz)の配列x,y,z,tの順序は、ncdumpしたときにでてくるfloat h(time, p, lat, lon)の並び順t,z,y,xの逆順でないといけません。（GrADSと逆。） |