Fortran で NetCDF データを読み込む

エンジニアであれば、Fortran で NetCDF データを読み込むなんてこともよくありますよね。だから今日はそのやり方を簡単に説明します！

前提

Fortran っていうのは、数値計算とかに使う計算が速い言語のことです。NetCDF は大気海洋系の人たちがよく使う、多次元配列データを複数まとめてメタデータといっしょに記録できるバイナリ形式の共通フォーマットです。

OS X ユーザーだったら、とりあえず Homebrew を使って

brew install gfortran
brew install netcdf --enable-fotran

とかしておきましょう。ポイントとしては、NetCDF をインストールしてから gfortran のバージョンが変わるとインクルードに失敗することがあるので、そういうときは面倒だから NetCDF 入れ直すといいです。

例

今回は例として NOAA の “NCEP/NCAR Reanalysis 1: Surface Flux” から、Daily の “Downward solar radiation flux” を使いましょう。適当に2006年のデータをダウンロードしておきました。

ncdump

NetCDF をインストールすると ncdump コマンドが使えるようになります。ダウンロードした拡張子が .nc のファイルのあるところにワーキングディレクトリを移動して、試しに

ncdump -c dswrf.sfc.gauss.2006.nc

としてみましょう。ごちゃっといろいろ表示されると思います。ncdump は NetCDF ファイルをダンプするコマンドで、-c オプションを使うことでヘッダ情報と次元情報を見ることができます。オプションなしだとデータも出ます。-h だとヘッダだけ見られます。

netcdf dswrf.sfc.gauss.2006 {
dimensions:
    lon = 192 ;
    lat = 94 ;
    time = UNLIMITED ; // (365 currently)
variables:
    float lat(lat) ;
        lat:units = "degrees_north" ;
        lat:actual_range = 88.542f, -88.542f ;
        lat:long_name = "Latitude" ;
        lat:standard_name = "latitude" ;
        lat:axis = "Y" ;
    float lon(lon) ;
        lon:units = "degrees_east" ;
        lon:long_name = "Longitude" ;
        lon:actual_range = 0.f, 358.125f ;
        lon:standard_name = "longitude" ;
        lon:axis = "X" ;
    double time(time) ;
        time:units = "hours since 1-1-1 00:00:0.0" ;
        time:long_name = "Time" ;
        time:actual_range = 17575512., 17584248. ;
        time:delta_t = "0000-00-01 00:00:00" ;
        time:avg_period = "0000-00-01 00:00:00" ;
        time:standard_name = "time" ;
        time:axis = "T" ;
    short dswrf(time, lat, lon) ;
        dswrf:long_name = "mean Daily Downward Solar Radiation Flux at surface" ;
        dswrf:unpacked_valid_range = -120.f, 1300.f ;
        dswrf:actual_range = -3.5f, 568.f ;
        dswrf:units = "W/m^2" ;
        dswrf:add_offset = 3156.5f ;
        dswrf:scale_factor = 0.1f ;
        dswrf:missing_value = 32766s ;
        dswrf:precision = 1s ;
        dswrf:least_significant_digit = 0s ;
        dswrf:GRIB_id = 204s ;
        dswrf:GRIB_name = "DSWRF" ;
        dswrf:var_desc = "Downward Solar Radiation Flux" ;
        dswrf:dataset = "NCEP Reanalysis Daily Averages" ;
        dswrf:level_desc = "Surface" ;
        dswrf:statistic = "Mean" ;
        dswrf:parent_stat = "Individual Obs" ;
        dswrf:valid_range = -32765s, -18565s ;

// global attributes:
        :Conventions = "COARDS" ;
        :title = "mean daily NMC reanalysis (2006)" ;
        :history = "created 2006/01/03 by Hoop (netCDF2.3)" ;
        :description = "Data is from NMC initialized reanalysis\n",
            "(4x/day).  It consists of T62 variables interpolated to\n",
            "pressure surfaces from model (sigma) surfaces." ;
        :platform = "Model" ;
        :references = "http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html" ;
}

-h したのをそのまま貼り付けました。

netcdf ファイル名 {
dimensions:
    次元名 = 次元の長さ ;
    ...
variables:
    型 変数名(次元名) ;
        変数名:属性名 = 属性値;
        ...
    ...
// global attributes:
    :グローバル属性名 = グローバル属性値 ;
    ...
}

こういう感じで表されています。わかりやすいですね。次元の中に長さが UNLIMITED となっているものがありますが、これは無制限に続けることができるということです。次元と変数で同じ名前のものがありますが、これはある次元で何番目だったら実際にはこの値というのが分かるようになっています。変数は次元を持ち、その個数に応じてN次元の配列になります。変数毎の属性はけっこう大事で、単位 units はもちろん、いろいろ重要な情報が記録されているので、必ず一通り確認しましょう。

Fortran

実際に Fortran で読み込んでみます。これには、NetCDF の Fortran 90 インターフェースを使います。基本的にはドキュメントにまとまっているので、必要に応じて確認しましょう。

先に実際のプログラムを載せておきます。

program netcdf_read
    use netcdf
    implicit none

    ! NetCDF file
    character (len=*), parameter :: FILE_NAME = "dswrf.sfc.gauss.2006.nc"
    integer :: ncid

    ! Dimensions
    integer, parameter :: NDIMS = 3
    integer :: lon_dim, lat_dim, time_dim
    real, allocatable :: lats(:), lons(:), times(:)

    ! Variable
    character(len=*), parameter :: VAR_NAME="dswrf"
    integer :: var_varid
    real, allocatable :: vars(:, :)
    real :: add_offset, scale_factor

    integer :: start(NDIMS), count(NDIMS)

    integer :: lat, lon, time
    integer :: err

    ! Open NetCDF file
    call check( nf90_open(FILE_NAME, nf90_nowrite, ncid) )

    ! Get dimension data
    call get_dimension(ncid, 'lon', lon_dim, lons)
    call get_dimension(ncid, 'lat', lat_dim, lats)
    call get_dimension(ncid, 'time', time_dim, times)

    ! Allocate variable
    allocate(vars(lat_dim, lon_dim), stat=err)
    if (err /= 0) print *, "Allocation request denied"
    ! Get variable id
    call check( nf90_inq_varid(ncid, VAR_NAME, var_varid) )
    ! Get add_offset and scale_factor
    call check( nf90_get_att(ncid, var_varid, 'add_offset', add_offset) )
    call check( nf90_get_att(ncid, var_varid, 'scale_factor', scale_factor) )

    ! Initialize count and start for reading data
    count = (/ lon_dim, lat_dim, 1 /)
    start = (/ 1, 1, 1 /)

    do time = 1, time_dim
        start(3) = time
        ! Read data
        call check( nf90_get_var(ncid, var_varid, vars, start = start, &
            count = count) )
        ! Restore values to original ones
        vars = vars * scale_factor + add_offset
        
        do lat = 1, lat_dim
            do lon = 1, lon_dim
                ! You can do anything here
                ! write(*, *) times(time), lats(lat), lons(lon), vars(lon, lat)
            end do
        end do
    end do

    ! Close NetCDF file
    call check( nf90_close(ncid) )

    ! If we got this far, everything worked as expected. Yipee! 
    print *, "*** SUCCESS reading example file ", FILE_NAME, "!"

contains
    subroutine get_dimension(ncid, name, dim, dims)
        integer,           intent( in) :: ncid
        character(len=*),  intent( in) :: name
        integer,           intent(out) :: dim
        real, allocatable, intent(out) :: dims(:)

        integer :: varid, dimid

        call check( nf90_inq_dimid(ncid, name, dimid) )
        call check( nf90_inquire_dimension(ncid, dimid, len=dim) )
        allocate(dims(dim), stat=err)
        if (err /= 0) print *, name, ": Allocation request denied"
        call check( nf90_inq_varid(ncid, name, varid) )
        call check( nf90_get_var(ncid, varid, dims) )
    end subroutine get_dimension

    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 netcdf_read

これを分割して説明を加えていきます。

program netcdf_read
    use netcdf
    implicit none

use netcdf することで NetCDF ライブラリの関数を使うと宣言しています。これはコンパイル時に必要です。また必ず implicit none してください。これがない Fortran は認めません。

ここからは変数宣言です。

    ! NetCDF file
    character (len=*), parameter :: FILE_NAME = "dswrf.sfc.gauss.2006.nc"
    integer :: ncid

ファイル名を定数にしています。定数を大文字にするのはよい習慣だと思います。ncid というのは、NetCDF ファイルを開いたあとに使うファイルの識別子です。

    ! Dimensions
    integer, parameter :: NDIMS = 3
    integer :: lon_dim, lat_dim, time_dim
    real, allocatable :: lats(:), lons(:), times(:)

次元に関する変数宣言をしています。後で使います。配列を allocatable にしているのは、実際の長さをデータに応じて変えるためです。( times は NetCDF で double になっているので real(8) にすべきですが、特に差がないので `real(4) にしています。 )

    ! Variable
    character(len=*), parameter :: VAR_NAME="dswrf"
    integer :: var_varid
    real, allocatable :: vars(:, :)
    real :: add_offset, scale_factor

NetCDF 内の変数についての変数宣言です。ここでも識別子となる var_varid があります。vars が2次元の配列なのは、実際には3次元だけど一度に読み込むのは2次元分にしたいからです。add_offset や scale_factor は後で説明します。

    integer :: start(NDIMS), count(NDIMS)

    integer :: lat, lon, time
    integer :: err

ファイルを読み込むループに使う変数とかです。

ここから実際の処理に移ります。

    ! Open NetCDF file
    call check( nf90_open(FILE_NAME, nf90_nowrite, ncid) )

nf90_open 関数で NetCDF ファイルを開きます。check サブルーチンは下の方で定義されていて後で説明しますが、エラーハンドリングのために用意しています。nf90_ というのが NetCDF の関数で、返り値として成否を返すので、これを check サブルーチンで確かめています。nf90_nowrite としているのは予期せず書き込んでしまうことを防ぐため、読み込み専用にしたいからです。これでファイルが開かれ、ncid にファイルと対応した識別子が代入されます。

    ! Get dimension data
    call get_dimension(ncid, 'lon', lon_dim, lons)
    call get_dimension(ncid, 'lat', lat_dim, lats)
    call get_dimension(ncid, 'time', time_dim, times)

ここでは次元の情報を後で定義するサブルーチン get_dimension で得ています。同じ処理が続くのでサブルーチンにしてあります。後で説明しますが、これで各次元の名前から長さと実際の値を取り出しています。

    ! Allocate variable
    allocate(vars(lat_dim, lon_dim), stat=err)
    if (err /= 0) print *, "Allocation request denied"

ここで allocatable にしてあった多次元配列 vars を実際に初期化します。次元の長さが分からないと配列の大きさも決まらないので、ここで初期化するのがよいでしょう。

    ! Get variable id
    call check( nf90_inq_varid(ncid, VAR_NAME, var_varid) )

ファイルと同じように、変数にも識別子があります。それをここで nf90_inq_varid を用いて取得します。

    ! Get add_offset and scale_factor
    call check( nf90_get_att(ncid, var_varid, 'add_offset', add_offset) )
    call check( nf90_get_att(ncid, var_varid, 'scale_factor', scale_factor) )

nf90_get_att は属性値を得るための関数です。変数の識別子と属性名を使って値を取り出します。

add_offset と scale_factor は、NetCDF データの容量削減と精度向上のために、予め変数の値を一定のルールで操作されているときに記録されている属性です。var = (var - add_offset) * 1/scale_vactor されているので、逆に var = var * scale_factor + add_offset してやる必要があります。そのためにここで値を得ておきます。

    ! Initialize count and start for reading data
    count = (/ lon_dim, lat_dim, 1 /)
    start = (/ 1, 1, 1 /)

ファイルを読み込むとき、一気に読み込むとメモリを多量に使ってしまうことになりがちです。今回は real としているので、仮にすべて読み込んだとしても25 MB程度に収まるはずですが、これより大きなデータを扱うことも多いでしょうから、今回も time 次元について分割して読み込んでいきます。

次に出てくる変数の値を読み込むところで、start と count を使います。変数の次元毎に読み込みの開始と量を設定できます。

    do time = 1, time_dim
        start(3) = time
        ! Read data
        call check( nf90_get_var(ncid, var_varid, vars, start = start, &
            count = count) )
        ! Restore values to original ones
        vars = vars * scale_factor + add_offset

時間次元のループの中で nf90_get_var 関数を使って、ファイル識別子と変数識別子から、vars 変数に値を読み込みます。start と count が使われています。その前に start(3) = time とすることで、読み込み場所をカウンタ変数 time によって設定できています。さらにここで vars = vars * scale_factor + add_offset としていて、配列すべてについて add_offset と scale_factor を考慮しています。

        do lat = 1, lat_dim
            do lon = 1, lon_dim
                ! You can do anything here
                ! write(*, *) times(time), lats(lat), lons(lon), vars(lon, lat)

いよいよここで実際の値を処理することができます。

            end do
        end do
    end do

    ! Close NetCDF file
    call check( nf90_close(ncid) )

    ! If we got this far, everything worked as expected. Yipee! 
    print *, "*** SUCCESS reading example file ", FILE_NAME, "!"

最後に nf90_close でファイルを閉じて、処理が終了したことを出力します。

この下は、このプログラム中で使うサブルーチンです。

contains
    subroutine get_dimension(ncid, name, dim, dims)
        integer,           intent( in) :: ncid
        character(len=*),  intent( in) :: name
        integer,           intent(out) :: dim
        real, allocatable, intent(out) :: dims(:)

        integer :: varid, dimid

get_dimension サブルーチンでは、次元名を受け取って次元の長さと各値を返します。nf90_inq_dimid で次元名から次元の識別子を取り出し、nf90_inquire_dimension で次元の長さを得ます。あとはそれを使って allocate して、変数を取り出すだけです。

        call check( nf90_inq_dimid(ncid, name, dimid) )
        call check( nf90_inquire_dimension(ncid, dimid, len=dim) )
        allocate(dims(dim), stat=err)
        if (err /= 0) print *, name, ": Allocation request denied"
        call check( nf90_inq_varid(ncid, name, varid) )
        call check( nf90_get_var(ncid, varid, dims) )
    end subroutine get_dimension

    subroutine check(status)
        integer, intent(in) :: status

check サブルーチンは、nf90_ で始まる関数が返すステータスが nf90_noerr かどうかを確認し、エラーであったときには nf90_strerror 関数でエラーに対応する文字列を得て出力し、プログラムを終了します。

        if (status /= nf90_noerr) then 
            print *, trim(nf90_strerror(status))
            stop "Stopped"
        end if
    end subroutine check  
end program netcdf_read

ここまでがプログラムの説明です。わかりやすいですね。

コンパイル

Homebrew でデフォルトの場所にインストールされた NetCDF を使うには

gfortran ファイル名.f90 -I/usr/local/include -L/usr/local/lib -lnetcdff -o ファイル名

みたいな感じにします。インクルードする netcdf.mod がある場所を -I オプションで、リンクする libnetcdff.a がある場所を -L オプションで指定し、-lnetcdff でリンクします。それぞれの場所は環境に合わせて変えてください。

まとめ

これで Fortran を使って NetCDF ファイルを読み込むことができましたね。今回のケースは変数などの名前が分かっているという前提ですが、そうでないときにもちゃんと読み込むことができます。あと書き込むことももちろんできますが、そういうのは今回は取り扱っていません。ドキュメント読めば何とかなります。

add_offset とか scale_factor とか、気をつけないと全然違う値になります。ちょっと面倒ですね。そういうのも勝手に何とかなっていたらいいですよね。というかそもそもいろいろ面倒ですね。

そこでオススメするのが Python の netcdf4 モジュールです。だいたいうまくやってくれるし、データが numpy で取り扱えるので、そのまま matplotlib でプロットもできます。記述量も減って簡単になるし、良いことしかないですね。ということで、Fortran で読み込むのやめて Python で読み込むといいですよ。