FortranでNetCDFを読み書きする

FortranでNetCDFを読み書きする。

基本

netCDFとは

自己記述的なデータフォーマットの一種。気象・気候分野の配布データとしてよく使われている。
Fortran, C, C++, pythonなどのAPIが用意されている
ここではFortran90での読み書きを対象とする。
netcdf-fortranの本家サイト

インストール

ソースからコンパイルしても良いし、apt-getやanacondaでインストールしてもよいと思う。

コンパイル

インストールするとnc-configコマンドも同時にインストールされる。
パスの指定等はこのコマンド用いると便利。

　 $ `nc-config --fc` ncread.f90 `nc-config --fflags` `nc-config --flibs`

netCDFモジュールを読む。

use netcdf

netCDFの関数の引数を受け取るsubroutineを用意しておく。

test.f90

program main 
use netcdf
implicit none
!! main code
contains
subroutine check(status)
 integer, intent (in) :: status
 if(status /= nf90_noerr) then 
   print *, trim(nf90_strerror(status))
   stop "Stopped"
 end if
end subroutine check
end program main

読む

基本的にncid：netcdfのファイル番号。
varid:netcdf内の変数の番号
を指定して読めば良い。

ncidを得る 。

character(len=1024) :: filename
integer :: ncid
call check( nf90_open(filename, nf90_nowrite, ncid) )

varid を得る（要ncid)

integer :: ncid, varid
call check ( nf90_inq_varid( ncid, "varname", varid) )

変数を読む。

integer :: start_nc(4), count_nc(4)
integer :: nx,ny,nz
start_nc=[1,1,1,1]
count_nc=[nx,ny,nz,1]
! data に読み取ったデータを入れる
! dataの型に注意
call check ( nf90_get_var(ncid, varid, data, start=start_nc, count=count_nc) )

ここでは読み取るデータの型に注意する。
ncdump -h などで事前に必要なデータの型を確認すること

ncdump での変数	Fortranの変数
double	real(8)
float	real(4)
int	integer(4)
short	integer(2)

となる。

また、多次元配列の場合ncdump等のコマンドで表示されるものと、逆に読まないといけない。これはFortranの配列アクセスが列優先で、c言語などと異なるからである。
例えばncdump で表示されるのが、

float air(time, lat, lon) ;

である場合には

integer,parameter :: ntime=12, nlat=73, nlon=144
real(4) :: air(nlon,nlat,ntime)

のようにする。ntime,nlat,nlonはそれぞれtime,lat,lonの要素数

さらに、次元の情報を得るやり方。

  ! 次元の大きさを得る
  call check( nf90_inq_dimid(ncid, 'lon', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nx) )

サンプルコード

対象は経度、緯度、高度（気圧）、時間の4次元あるデータを読み、各気圧層の最大最小を出力するプログラム。
サンプルデータは以下から得た。
[https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html]

サンプルデータの中身を確認する

$ ncdumo -h air.mon.nc

サンプルコード本体

read.f90

program read
  use netcdf
  implicit none
  ! dimension
  integer :: nx, ny, nz
  character( len=1024 ) :: f_in
  real(4), allocatable :: data(:,:,:)
  ! for netcdf
  integer :: ncid, varid, dimid
  integer :: start_nc(4), count_nc(4)
  ! doloop
  integer :: k, n
  !
  !
  f_in = "air.mon.mean.nc"
  call check( nf90_open( trim( f_in ), nf90_nowrite, ncid) )
  ! 次元の大きさを得る
  call check( nf90_inq_dimid(ncid, 'lon', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nx) )
  call check( nf90_inq_dimid(ncid, 'lat', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=ny) )
  call check( nf90_inq_dimid(ncid, 'level', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nz) )
  allocate( data(nx,ny,nz) )
  ! varid を得る
  call check( nf90_inq_varid(ncid, "air", varid))
  n = 1
  ! n-step目の変数を読む
  start_nc = [1,1,1,n] 
  count_nc = [nx,ny,nz,1]
  call check( nf90_get_var(ncid, varid, data, start=start_nc, count=count_nc) )
  do k = 1, nz
     write(*,*) k, maxval(data(:,:,k)), minval(data(:,:,k))
  end do
  ! file を閉じる
  call check( nf90_close(ncid) )
contains
  subroutine check( status )
    integer, intent (in) :: status
    if(status /= nf90_noerr) then 
       print *, trim(nf90_strerror(status))
       stop "Stopped"
    end if
  end subroutine check
end program read

書き込む

書き込むファイルを開く(ncidを得る)。

特にこだわりがなければnetcdf4形式で良いと思う。

  call check( nf90_create( trim(f_out), NF90_HDF5, ncido) ) 

変数を定義する(define mode)

次元の定義

次元を定義しその次元のdimidを得る。次の変数の定義で用いる。

call check( nf90_def_dim(ncido, 'lat', ny, lat_dim_id) )
! 1次元のみNF90_UNLIMITEDとして次元長を指定しなくても良い。RECORD次元という。
call check( nf90_def_dim(ncido, 'time', NF90_UNLIMITED, time_dim_id) )

変数の定義

 ! 1次元量の場合
 call check( nf90_def_var(ncido, 'lat', NF90_REAL, lat_dim_id, lat_id) )
 !　多次元量は配列でdimidを指定する。
! deflate_levelを指定すると圧縮がかかる。netcdf4のみ
 call check( nf90_def_var(ncido, "air", NF90_REAL, (/ lat_dim_id, level_dim_id, time_dim_id/), &
& varido, deflate_level = 1 ) )

変数にattributionを付ける

attribution(属性)を付ける。変数の説明(単位、名前 etc.）を付加できる。
ncid、 varidを指定し書き込む。下の例は"units" として 'degrees_north'を書き込む例。

call check( nf90_put_att(ncido, lat_id, 'units', 'degrees_north') )

変数定義モードを終える。

  call check( nf90_enddef(ncido) )

変数を書き込む。

ncidとvaridを指定して、変数を入れれば良い。
時間ステップ毎に入れる際などはstartとcountを指定する。

 call check( nf90_put_var(ncido, lat_id, lat ) )
 start_nco = [1,1,n]
 count_nco = [ny,nz,1]
 call check( nf90_put_var(ncido, varido, out, start=start_nco, count=count_nco) )

書き込むファイルを閉じる。

 call check( nf90_close(ncido) )

サンプルコード

対象は経度、緯度、高度（気圧）、時間の4次元あるデータを読み、経度方向に平均をとって、緯度、高度（気圧）、時間の3次元のデータに変換する。

write.f90

program read
  use netcdf
  implicit none
  ! dimension
  integer :: nx, ny, nz, nt
  ! file
  character( len=1024 ) :: f_in, f_out
  real(4), allocatable :: data(:,:,:)
  real(4), allocatable :: out(:,:)
  real(4), allocatable :: lat(:), level(:)
  real(8) :: time(1)
  ! for netcdf input
  integer :: ncid, varid, dimid
  integer :: start_nc(4), count_nc(4)
  ! for netcdf output
  integer :: ncido, varido
  integer :: lat_dim_id, level_dim_id, time_dim_id
  integer :: lat_id, level_id, time_id
  integer :: start_nco(3), count_nco(3)
  ! doloop
  integer :: k, n
  !
  !
  f_in = "air.mon.mean.nc"
  f_out = "air.mon.zonmean.nc"
  !
  call check( nf90_open( trim( f_in ), nf90_nowrite, ncid) )
  ! 次元の大きさを得る
  call check( nf90_inq_dimid(ncid, 'lon', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nx) )
  call check( nf90_inq_dimid(ncid, 'lat', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=ny) )
  allocate(lat(ny))
  call check( nf90_get_var(ncid, dimid, lat) )
  call check( nf90_inq_dimid(ncid, 'level', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nz) )
  allocate(level(nz))
  call check( nf90_get_var(ncid, dimid, level) )
  call check( nf90_inq_dimid(ncid, 'time', dimid) )
  call check( nf90_inquire_dimension(ncid, dimid, len=nt) )
  ! 時間次元だけ読み込む
  !
  allocate( data(nx,ny,nz) )
  ! varid を得る
  call check( nf90_inq_varid(ncid, "air", varid))
  ! 読む設定はここまで。
  ! 書く設定
  ! 
  call check( nf90_create( trim(f_out), NF90_HDF5, ncido) )
  ! 次元を定義する
  call check( nf90_def_dim(ncido, 'lat', ny, lat_dim_id) )
  Call check( nf90_def_dim(ncido, 'level', nz, level_dim_id) )
  call check( nf90_def_dim(ncido, 'time', NF90_UNLIMITED, time_dim_id) )
  ! 変数を定義する
  call check( nf90_def_var(ncido, 'lat', NF90_REAL, lat_dim_id, lat_id) )
  call check( nf90_def_var(ncido, 'level', NF90_REAL, level_dim_id, level_id) )
  call check( nf90_def_var(ncido, 'time', NF90_DOUBLE, time_dim_id, time_id) )
  call check( nf90_def_var(ncido, "air", NF90_REAL, (/ lat_dim_id, level_dim_id, time_dim_id/),&
       & varido, deflate_level = 1 ) )
  ! 変数に attributionをつける　
  call check( nf90_put_att(ncido, lat_id, 'units','degrees_north') )
  call check( nf90_put_att(ncido, level_id, 'units','millibar') )
  call check( nf90_put_att(ncido, time_id, 'units', 'hours since 1800-01-01 00:00:0.0' ) )
  ! define モードを終了する。
  call check( nf90_enddef(ncido) )
  call check( nf90_put_var(ncido, lat_id, lat ) )
  call check( nf90_put_var(ncido, level_id, level ) )
  allocate( out(ny,nz) )
  !
  do n = 1,nt 
  ! n-step目の変数を読む
     start_nc = [1,1,1,n]
     count_nc = [nx,ny,nz,1]
     call check( nf90_get_var(ncid, varid, data, start=start_nc, count=count_nc) )
     call check( nf90_get_var(ncid, dimid, time, start=[n], count=[1]) )
     out = sum(data,dim=1)/nx
     start_nco = [1,1,n]
     count_nco = [ny,nz,1]
     call check( nf90_put_var(ncido, varido, out, start=start_nco, count=count_nco) )
     call check( nf90_put_var(ncido, time_id, time, start=[n], count=[1]) )
  end do
  ! file を閉じる
  call check( nf90_close(ncid) )
  call check( nf90_close(ncido) )
  ! deallocate
  deallocate( data, out, lat, level )

contains

  subroutine check( status )
    integer, intent (in) :: status
    if(status /= nf90_noerr) then 
       print *, trim(nf90_strerror(status))
       stop "Stopped"
    end if
  end subroutine check

end program read

#おまけ

上記はcdo(Climate Data Operators)を使うと一瞬である。
cdoのページ[https://code.mpimet.mpg.de/projects/cdo/]
気候の計算で使うような処理の多くはこのコマンドで扱える。

$ cdo zonmean air.mon.mean.nc air.mon.zonmean.nc