[Exp2023]Fortran を用いた NetCDF データの解析
はじめに
ここまでに, ruby や GPhys を使って NetCDF 形式のデータの簡単な解析について説明しました. ruby や GPhys に習熟すれば, 既に説明した以上に複雑なデータ解析・処理も行うことができます. しかしながら, ruby や GPhys での複雑な処理は若干敷居が高いかもしれません. それに比べると, 惑星学実験実習の基礎 IIで Fortran を既に学んでいますので, データの入出力さえできれば Fortran でのデータ解析・処理は比較的容易です. したがって, ここでは Fortran を用いた NetCDF データの入出力について説明します. (Fortran での演算の方法については, 惑星学実験実習の基礎 II の内容を復習すると良いでしょう.)
Fortran プログラムでの NetCDF データの入出力の概要
NetCDF データはメタデータを含む独自の形式ですから, Fortran の標準的な命令である read/write では入出力することができません. したがって, NetCDF データを Fortran プログラムで入出力するためのプログラムが必要です.
NetCDF の開発元である Unidata は, Fortran のプログラムから NetCDF ファイルを入出力するためのライブラリを提供しています. (「ライブラリ」は一連のプログラムを集めたものです. 例えば, Fortran であれば read/write に相当する, NetCDF ファイルを入出力するためのサブルーチンや関数をまとめたプログラム群のことです. なお, Unidata は C, C++, Java のためのライブラリも提供しています.) その Unidata が提供する NetCDF の Fortran ライブラリを用いることで, 自分で読み込みサブルーチン・関数を新たに作ることなく Fortran プログラムで NetCDF データの入出力が可能です.
しかし, Unidata が提供するライブラリは考えられるすべての入出力ができるものの, 「よくある使い方」だけでよいユーザにとっては使いにくいかもしれません. そこでここでは, ユーザの利便性を考慮して地球対電脳倶楽部の有志が開発した NetCDF データ入出力のための Fortran ライブラリ (gtool5) を用いた入出力方法について説明します. gtool5 は, 内部で Unidata 提供のライブラリを使用しています.
gtool5 のインストール
gtool5 は下のコマンドでインストールすることができます.
$ sudo apt install gtool5
Fortran プログラムでの NetCDF ファイルの入力
ここでは, 既にダウンロードした NetCDF ファイル air.2019.nc を用いて Fortran での入力を説明します.
下のプログラムを実行してみましょう.
01 ! Sample program for gtool_history/gtool5 02 ! 03 program nc_i 04 05 use gtool_history ! モジュール指定 06 07 implicit none 08 09 character(128) :: InFileName ! 入力ファイル名 10 11 integer, parameter :: NLon = 144 ! 経度格子点数 12 integer, parameter :: NLat = 73 ! 緯度格子点数 13 integer, parameter :: NLevel = 17 ! 高度(圧力)格子点数 14 integer, parameter :: NTime = 365 ! 時間格子点数 15 16 real(4) :: x_Lon (NLon ) ! 経度配列 17 real(4) :: y_Lat (NLat ) ! 緯度配列 18 real(4) :: z_Level (NLevel) ! 高度(圧力)配列 19 real(4) :: t_Time (NTime ) ! 時間配列 20 real(4) :: xyz_Temp (NLon, NLat, NLevel) ! 温度配列 21 22 integer :: it 23 integer :: i 24 integer :: j 25 integer :: k 26 27 character(64) :: range 28 29 30 InFileName = "air.2019.nc" ! 入力ファイル名の設定 31 32 33 call HistoryGet(InFileName, 'lon' , x_Lon ) ! [入力] 経度を読む 34 call HistoryGet(InFileName, 'lat' , y_Lat ) ! [入力] 緯度を読む 35 call HistoryGet(InFileName, 'level', z_Level) ! [入力] 圧力を読む 36 call HistoryGet(InFileName, 'time' , t_Time ) ! [入力] 時間を読む 37 38 39 do it = 1, NTime 40 write( 6, * ) it, t_Time(it) 41 42 write( range, '(a,i5.5)' ) 'time=^', it ! 出力用の時刻情報文字列の 作成 43 write( 6, * ) range 44 45 ! [入力] it 番目のデータを読む 46 call HistoryGet(InFileName, 'air', xyz_Temp, range = range) 47 48 ! i=NLon/2, j=NLat/2, k=1 における温度の出力 49 i = NLon/2 50 j = NLat/2 51 k = 1 52 write( 6, * ) t_Time(it), x_Lon(i), y_Lat(j), z_Level(k), xyz_Temp(i,j,k) 53 54 end do 55 56 stop 57 end program nc_i
プログラムを実行するには,
$ gt5frt -o nc_i nc_i.f90
でコンパイルし,
$ ./nc_i
で実行します. 実行すると下のように表示されるでしょう.
$ ./nc_i
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@lon
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@lat
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@level
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@time
1 1919712.00
time=^00001
*** MESSAGE [HistoryGetReal3] *** Input air.2019.nc@air,time=1919712.
1919712.00 177.500000 2.50000000 1000.00000 300.649963
2 1919736.00
time=^00002
*** MESSAGE [HistoryGetReal3] *** Input air.2019.nc@air,time=1919736.
1919736.00 177.500000 2.50000000 1000.00000 301.649994
3 1919760.00
time=^00003
*** MESSAGE [HistoryGetReal3] *** Input air.2019.nc@air,time=1919760.
1919760.00 177.500000 2.50000000 1000.00000 301.625000
(中略)
365 1928448.00
time=^00365
*** MESSAGE [HistoryGetReal3] *** Input air.2019.nc@air,time=1928448.
1928448.00 177.500000 2.50000000 1000.00000 301.024963
[入力] プログラムの説明
上のプログラムは, air.2019.nc を読み込み, i=NLon/2, j=NLat/2, k=1 の場所の温度を表示します.
[入力] プログラムの内容の大まかな説明
ここではプログラムの内容を大雑把に説明します.
プログラムでは下のような処理をしています. (行頭の数字は行番号です.)
- 05-27: 使用するモジュールや変数を準備しています.
- 30: 入力ファイル名を指定しています.
- 33-36: 入力ファイルから経度軸, 緯度軸, 圧力(高度)軸, 時間軸のデータを入力しています.
- 39-54:
入力ファイルから時間軸方向に一つのデータずつを読み, 温度を画面に表示しています.
- 42: 何番目の時刻のデータを入出力するかを指定するための「ラベル」文字列を用意します.
- 46: 入力ファイルから時間一つ分の温度データ (3 次元) を読んでいます.
- 52: 読んだ温度のある場所の値を画面に出力しています.
例えばこのプログラムを基にして, 39-54 行目を変更することで自分がやりたいデータ処理ができるでしょう.
[入力] プログラムの詳細な説明
ここではプログラムの内容をより詳しく説明します.
プログラムでは下のような処理をしています. (行頭の数字は行番号です.)
- 05: 使用するモジュールの指定. ここでは gtool5 のモジュールを使うことを宣言しています.
- 09-27:
変数を宣言しています.
- 20: 読むデータを格納する配列を宣言しています. このとき, 次元の順番に気を付けましょう. 今回読むデータの次元は 'lon', 'lat', 'level', 'time' ですが, netcdf ファイルの内容を ncdump で確認すると, 次元は 'time', 'level', 'lat', 'lon' の順番に書かれるでしょう. しかし, Fortran の配列での順番は逆順になりますので注意しましょう. (なお, C 言語では Fortran とは逆 (ncdump の表示の順) となります.)
- 30: 入力ファイル名, 出力ファイル名を指定しています.
- 33-36: 入力ファイルから経度軸, 緯度軸, 圧力(高度)軸, 時間軸のデータを入力しています.
- 39-54:
入力ファイルから時間軸方向に一つのデータずつを読み, 温度を画面に表示しています.
- 入力 NetCDF ファイルに保存されているすべての時間のデータを一度に読み込もうとするとメモリサイズに収まらずにエラーになることが少なくありません. 時間 1 つのデータごとに読み込み, 演算して, 出力することを繰り返すことで, メモリ不足を回避することができます.
- 各行は下のような処理を行っています.
- 42:
何番目の時刻のデータを入出力するかを指定するための「ラベル」文字列を用意します.
write は端末に数値や文字を出力するときに使うことが多いでしょう. 例えば, 端末に出力する際には,
write( 6, * ) ...
のように使いました.
- "6" は出力先を表し, 端末 (画面) に対応するのでした.
- "*" は「適当に」出力することの指定でした. (詳細は惑星学実験実習の基礎 IIの資料を参照してください.)
この出力先の指定には他のものを指定することができ, 上のプログラムでは文字列 "range" を指定することで, 文字列変数に出力, つまり, 文字列変数に書き込み = 代入しています.
ここで指定しているラベルは下のような形式です.
time=^N
ここで, N は配列の要素の番号 (インデクス) です. "^" は要素番号を表します. "^" を付けない場合には, 配列に収められた値 (ここでは時間の値) を表します.
- 46: 入力ファイルから時間一つ分の温度データ (3 次元) を読んでいます.
- 52: 読んだ温度のある場所の値を画面に出力しています.
- 42:
何番目の時刻のデータを入出力するかを指定するための「ラベル」文字列を用意します.
例えばこのプログラムを基にして, 39-54 行目を変更することで自分がやりたいデータ処理ができるでしょう. もし複数の NetCDF ファイルから複数の物理量を読んでデータ処理をする必要があれば, 上で使われている HistoryGet を参考にして読み込みを追加すると良いでしょう.
Fortran プログラムでの NetCDF ファイルの入出力
ここでは, 上のプログラムに出力を追加して, Fortran での入出力を説明します.
下のプログラムを実行してみましょう.
01 ! Sample program for gtool_history/gtool5
02 !
03 program nc_io
04
05 use gtool_history ! モジュール指定
06
07 implicit none
08
09 character(128) :: InFileName ! 入力ファイル名
10 character(128) :: OutFileName ! 出力ファイル名
11
12 integer, parameter :: NLon = 144 ! 経度格子点数
13 integer, parameter :: NLat = 73 ! 緯度格子点数
14 integer, parameter :: NLevel = 17 ! 高度(圧力)格子点数
15 integer, parameter :: NTime = 365 ! 時間格子点数
16
17 real(4) :: x_Lon (NLon ) ! 経度配列
18 real(4) :: y_Lat (NLat ) ! 緯度配列
19 real(4) :: z_Level (NLevel) ! 高度(圧力)配列
20 real(4) :: t_Time (NTime ) ! 時間配列
21 real(4) :: xyz_Temp (NLon, NLat, NLevel) ! 温度配列
22 real(4) :: xyz_PotTemp(NLon, NLat, NLevel) ! 温位配列
23
24 real(4) :: Kappa = 0.29 ! R/Cp
25 real(4) :: P0 = 1000.0 ! 基準圧力
26
27 integer :: it
28 integer :: i
29 integer :: j
30 integer :: k
31
32 character(64) :: range
33
34
35 InFileName = "air.2019.nc" ! 入力ファイル名の設定
36 OutFileName = "output.nc" ! 出力ファイル名の設定
37
38
39 call HistoryGet(InFileName, 'lon' , x_Lon ) ! [入力] 経度を読む
40 call HistoryGet(InFileName, 'lat' , y_Lat ) ! [入力] 緯度を読む
41 call HistoryGet(InFileName, 'level', z_Level) ! [入力] 圧力を読む
42 call HistoryGet(InFileName, 'time' , t_Time ) ! [入力] 時間を読む
43
44
45 ! [出力] 出力ファイル作成 (開く)
46 call HistoryCreate( &
47 & file=OutFileName, title='atmospheric temperature', &
48 & source='NCEP/NCAR Reanalysis data', &
49 & institution='Department of Planetology, Kobe University', &
50 & dims=(/'lon ','lat ','level','time '/), &
51 & dimsizes=(/NLon,NLat,NLevel,NTime/), &
52 & longnames=(/'longitude','latitude ','level ', 'time '/), &
53 & units=(/'degrees_east ', &
54 & 'degrees_north ', &
55 & 'hPa ', &
56 & 'hours since 1800-01-01 00:00:0.0'/) &
57 & )
58
59 call HistoryPut('lon' ,x_Lon ) ! [出力] 次元変数出力, 経度
60 call HistoryPut('lat' ,y_Lat ) ! [出力] 次元変数出力, 緯度
61 call HistoryPut('level',z_Level) ! [出力] 次元変数出力, 高度(圧力)
62 call HistoryPut('time' ,t_Time ) ! [出力] 次元変数出力, 時間
63
64
65 ! [出力] 変数定義, 温位
66 call HistoryAddVariable( &
67 & varname='pottemp', dims=(/'lon ','lat ','level','time '/), &
68 & longname='potential temperature', units='K', xtype='real' &
69 & )
70
71 do it = 1, NTime
72 write( 6, * ) it, t_Time(it)
73
74 write( range, '(a,i5.5)' ) 'time=^', it ! 出力用の時刻情報文字列の作 成
75 write( 6, * ) range
76
77 ! [入力] it 番目のデータを読む
78 call HistoryGet(InFileName, 'air', xyz_Temp, range = range)
79
80 ! 温位の計算
81 do k = 1, NLevel
82 do j = 1, NLat
83 do i = 1, NLon
84 xyz_PotTemp(i,j,k) = xyz_Temp(i,j,k) * ( P0 / z_Level(k) )**Kappa
85 end do
86 end do
87 end do
88
89 ! [出力] it 番目のデータを出力, 温位
90 call HistoryPut('pottemp', xyz_PotTemp, range=range)
91 end do
92
93 call HistoryClose ! [出力] 出力ファイルを閉じる
94
95 stop
96 end program nc_io
プログラムを実行するには,
$ gt5frt -o nc_io nc_io.f90
でコンパイルし,
$ ./nc_io
で実行します. 実行すると下のように表示されるでしょう.
$ ./nc_io
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@lon
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@lat
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@level
*** MESSAGE [HistoryGetReal1] *** Input air.2019.nc@time
*** MESSAGE [GDNcFileOpen] *** "output.nc" is overwritten.
*** MESSAGE [HistoryCreate1] *** "output.nc" is created (origin=0. [])
1 1919712.00
time=^00001
*** MESSAGE [HistoryGetReal3] *** Input air.2019.nc@air,time=1919712.
2 1919736.00
time=^00002
(中略)
364 1928424.00
time=^00364
*** MESSAGE [HistoryGetReal3] *** Input air.2019.nc@air,time=1928424.
*** MESSAGE [HistoryClose] *** "output.nc" is closed
[入出力] プログラムの説明
上のプログラムは, air.2019.nc を読み込み, その中に保存されている温度を使って温位を計算して NetCDF ファイル (ここでは output.nc) に保存します.
- 温位についてはここでは詳しく説明しませんが, エントロピーに相当する物理量です. 大気科学では「大気塊を断熱的に基準圧力 (地表面) まで移動させたときの温度」という意味で説明されることが多いでしょう.
[入出力] プログラムの内容の大まかな説明
ここではプログラムの内容を大雑把に説明します.
プログラムでは下のような処理をしています. (行頭の数字は行番号です.)
- 05-32: 使用するモジュールや変数を準備しています.
- 35-36: 入力ファイル名, 出力ファイル名を指定しています.
- 39-42: 入力ファイルから経度軸, 緯度軸, 圧力(高度)軸, 時間軸のデータを入力しています.
- 45-69:
出力ファイルを準備しています.
- 45-57: 出力 NetCDF ファイルを作成しています.
- 59-62: 出力ファイルに経度軸, 緯度軸, 圧力(高度)軸, 時間軸のデータを出力しています.
- 65-69: 出力する温位データを登録し, メタデータを設定しています.
- 71-91:
入力ファイルから時間軸方向に一つのデータずつを読み, 温位を計算して出力しています.
- 74: 何番目の時刻のデータを入出力するかを指定するための「ラベル」文字列を用意します.
- 78: 入力ファイルから時間一つ分の温度データ (3 次元) を読んでいます.
- 81-87: 入力ファイルから読んだ温度を使って温位を計算しています.
- 90: 計算した温位を NetCDF ファイルに出力しています.
- 93: 出力 NetCDF ファイルを閉じています.
例えばこのプログラムを基にして, 71-91 行目を変更することで自分がやりたいデータ処理ができるでしょう.
[入出力] プログラムの詳細な説明
ここではプログラムの内容をより詳しく説明します.
プログラムでは下のような処理をしています. (行頭の数字は行番号です.)
- 05: 使用するモジュールの指定. ここでは gtool5 のモジュールを使うことを宣言しています.
- 09-32:
変数を宣言しています.
- 21-22: 読むデータを格納する配列を宣言しています. このとき, 次元の順番に気を付けましょう. 今回読むデータの次元は 'lon', 'lat', 'level', 'time' ですが, netcdf ファイルの内容を ncdump で確認すると, 次元は 'time', 'level', 'lat', 'lon' の順番に書かれるでしょう. しかし, Fortran の配列での順番は逆順になりますので注意しましょう. (なお, C 言語では Fortran とは逆 (ncdump の表示の順) となります.)
- 35-36: 入力ファイル名, 出力ファイル名を指定しています.
- 39-42: 入力ファイルから経度軸, 緯度軸, 圧力(高度)軸, 時間軸のデータを入力しています.
- 45-57:
出力 NetCDF ファイルを準備しています.
- NetCDF ファイルに出力する軸情報などのメタデータは下のオプションで指定します.
- file : 出力 NetCDF ファイル名
- title : データのタイトル
- source : データの出どころ
- institution : データ作成機関 / データ作成者の所属機関
- dims : NetCDF ファイルが含む軸/次元
- (/'lon ','lat ','level','time '/) は配列になっていて, このファイルが 4 つの軸, lon, lat, level, time を含むことを表しています.
- なお, 'lon ' のように "n" の後に空白を二つ入力しているのは, すべての軸の名前の長さを揃えるためです. 今回の場合は level が 5 文字の為, 他の軸の名前にも空白を入れて 5 文字に揃えています.
もし (/'lon','lat','level','time'/) のように文字列を揃えないと, コンパイル時に下のようなエラーが出ます.
Error: Different CHARACTER lengths (4/5) in array constructor at (1)
- dimsizes : dims で与えた軸/次元のそれぞれの長さ
- これも (/NLon,NLat,NLevel,NTime/) は配列になっています.
- longnames : 軸/次元の名前のより詳しい名前
- これも (/'longitude','latitude ','level ', 'time '/) は配列になっており, すべての名前の文字数を揃えるために空白を使っています.
- units : 軸/次元の単位
- これも配列になっており, すべての名前の文字数を揃えるために空白を使っています. 今回はそれぞれが長いため, "&" を使って改行しています.
- NetCDF ファイルに出力する軸情報などのメタデータは下のオプションで指定します.
- 59-62: 出力ファイルに経度軸, 緯度軸, 圧力(高度)軸, 時間軸のデータを出力しています.
- 65-69: 出力する温位データを登録し, メタデータを設定しています.
- 71-91:
入力ファイルから時間軸方向に一つのデータずつを読み, 温位を計算して出力しています.
- 入力 NetCDF ファイルに保存されているすべての時間のデータを一度に読み込もうとするとメモリサイズに収まらずにエラーになることが少なくありません. 時間 1 つのデータごとに読み込み, 演算して, 出力することを繰り返すことで, メモリ不足を回避することができます.
- 各行は下のような処理を行っています.
- 74:
何番目の時刻のデータを入出力するかを指定するための「ラベル」文字列を用意します.
write は端末に数値や文字を出力するときに使うことが多いでしょう. 例えば, 端末に出力する際には,
write( 6, * ) ...
のように使いました.
- "6" は出力先を表し, 端末 (画面) に対応するのでした.
- "*" は「適当に」出力することの指定でした. (詳細は惑星学実験実習の基礎 IIの資料を参照してください.)
この出力先の指定には他のものを指定することができ, 上のプログラムでは文字列 "range" を指定することで, 文字列変数に出力, つまり, 文字列変数に書き込み = 代入しています.
ここで指定しているラベルは下のような形式です.
time=^N
ここで, N は配列の要素の番号 (インデクス) です. "^" は要素番号を表します. "^" を付けない場合には, 配列に収められた値 (ここでは時間の値) を表します.
- 78: 入力ファイルから時間一つ分の温度データ (3 次元) を読んでいます.
- 81-87: 入力ファイルから読んだ温度を使って温位を計算しています.
- 90: 計算した温位を NetCDF ファイルに出力しています.
- 74:
何番目の時刻のデータを入出力するかを指定するための「ラベル」文字列を用意します.
- 93: 出力 NetCDF ファイルを閉じています.
例えばこのプログラムを基にして, 71-91 行目を変更することで自分がやりたいデータ処理ができるでしょう. もし複数の NetCDF ファイルから複数の物理量を読んでデータ処理をする必要があれば, 上で使われている HistoryGet を参考にして読み込みを追加すると良いでしょう.
gt5frt についての補足
上の説明ではプログラムをコンパイルするために gt5frt を使いました.
これは Fortran コンパイラ gfortran に, NetCDF ライブラリや gtool5 ライブラリを用いるためのオプションを付けた「別名」コマンドです. 具体的には下のようなオプションを指定しています.
- -I/usr/include
- -I/usr/lib/x86_64-linux-gnu/gtool5/include
- -L/usr/lib/x86_64-linux-gnu/gtool5/lib
- -L/usr/lib/x86_64-linux-gnu
- -lgtool5
- -lnetcdf
- -lnetcdff
- -fPIC
これらのオプションは, -I, -L, -l -f にオプションの「値」を組み合わせた形になっています. 興味があれば -I, -L, -l, -f の意味を調べてみると良いでしょう. (man コマンドを使うのも良いでしょうし, 「gfortran」「オプション」などをキーワードにして検索してみるのも良いでしょう.)
より詳しい情報
gtool についてのより詳しい情報は下の場所で見つけることができるでしょう.
Keyword(s):
References:[[Exp2023]スケジュール表・各回資料]