Last update: 01 May 2018
追記メモ(2018年5月1日)
2018年5月現在、Windows環境でのSEIBの実行には、フリーのgfortranコンパイラであるMinGW-w64をお使い頂くのがお勧めです。以下のサイトから入手できます。コンパイラコマンドは「g95 」ではなくて「gfortran」です。MinGW-w64は、インストーラーがgfortranコマンドのあるディレクトリをPATHに自動追加してくれませんので、スタートメニュのリンクからコマンドプロンプトを起動する必要があります。必要に応じて、PATHに追加してしまってもよいかもしれません。
https://sourceforge.net/projects/mingw-w64/
なお、TDM-GCCのgfortranは、実行時にopen文・write文でSegmentation faultが出てしまうため（これは、SEIBのバグによるものでは無く、このコンパイラのバグに起因しています）、ご使用になれません。
(1) G95のインストール
サイトhttp://www.g95.org/のトップページ→Downloads→Binaries→Information for Windows users→項目"g95-MinGW.exe"内のhttp://ftp.g95.org/g95-MinGW.exeをクリック、インストーラーをダウンロードし、実行。インストール先のフォルダの問いについては、適当に邪魔にならないところを指定（"c:\g95"とか）。「Install MinGW utilities and libs?」に対してはYes。基本的に全てYesとOkでインストールを進める。
G95のインストールについての更に詳しい解説は、次のページをご覧下さい。http://d.hatena.ne.jp/arakik10/20090213/p1
何らかの理由でファイルがサイトから削除されているようです（2018年5月時点に確認）。以下のリンク集に、コピーを保存している第三者のサイトが紹介されていますが、ウィルス・スパイウェアの類で無いかどうか十二分にご確認の上で、利用されることをお勧めします。
http://www.filewatcher.com/m/g95-MinGW.exe.5509899-0.html
(3) Fortran実行環境の構築
適当な場所に、fortranという名前のフォルダを作成。このフォルダ内に、コマンドプロンプトのショートカットを作成。このショートカットを右クリックし、プロパティを選択、「ショートカット」タブ内の「作業フォルダー」を空欄にする（これで、ショートカットの置かれているフォルダが、作業フォルダになる）。
次に、ここでFortranコードのコンパイルと実行ができる事を確認する。テキストファイルを新規作成し、以下のコードをコピペし、ファイル名をtest.f90と変更し、これを上記のフォルダに入れる。
Program sample_code write(*,*) 'hello world!' stop End Program
先に作成したショートカットをダブルクリックしてコマンドプロンプトを開き、ここで「g95 test.f90」と入力すると、このフォルダ中に実行ファイルa.exeができる。「a.exe」と入力して、正しく実行されることを確認する。なお、コマンドプロンプトにおける基本コマンドは以下の通り。
| exit | コマンドプロンプトから抜ける |
| dir | コマンド内のファイルとフォルダの一覧を表示(Unixのlsに相当) |
| cd | カレントディレクトリを移る |
| del | ファイルの消去 |
(5) SEIB-DGVMの実行
SEIB-DGVMの最新版Distribution Packageをダウンロードし、解凍する。解凍されたフォルダ内のフォルダ"code"を中身ごと、上のフォルダ"Fortran"に移動させる。コマンドプロンプトにおいて「cd code」に続いて「g95 start_point.f90」と入力すると、コンパイルが始まる。この環境においては、コンパイルにかなり時間がかかるので、本格的な開発にはLinux上でIntel Fortran Compilerを使うのがお勧め。コンパイルが終了したら、「a.exe」で実行される。実行終了後、フォルダ"code"内にシミュレーションの結果ファイル群ができていることを確認する。
(6) 結果ファイルの可視化
このシミュレーションで、どのような仮想林分が発達したのか、三次元で可視化するためにはPOV-Rayを使います。また、物質やエネルギー循環などの詳細を可視化するためにはSEIB-Viewerを使います。具体的な手順については、SEIB-DGVMのDistribution Package内のreadme_j.txt内に説明があります。本格的な解析・開発のためには、コードoutput.f90を手直しして、オリジナルな結果ファイルを出力させと良いでしょう。
Last update: 14 July 2008
気象データなど、座標軸上に配置されたデータセットは、しばしばnetCDF(network Common Data Form)というフォーマットで配布されています。netCDFには、データ本体に加えて、座標系や変数名などの情報が含まれていますので、アプリケーションによっては、いろいろ数字を指定しなくてもファイルを読むだけで可視化できたりと、メリットも大きいようです。しかし、私のように「ちょっとデータをFortranやExcelで処理したい」などといったライトユーザーにとっては、ヘッダの構造など学ばなければならない事が多くて、なかなか悩ましい存在なのです。というわけで、統計解析や作図などに広く使われるフリーソフトのRを用いて、netCDFに含まれるデータを読み出したり可視化したりする方法を同じ研究室の伊勢さんより教わりましたので、備忘録も兼ねて技術メモを書いてみた次第です。
サンプルデータとしてCPC Merged Analysis of Precipitation (CMAP)の世界降水量データを使います。このページ[外部リンク]の「Download/PlotData」欄の一番最初に出てくるprecip.mon.mean.ncというやつをダウンロードしてきて下さい。これは、複数の人工衛星データを組み合わせて、1979年からの月平均降水量を全球で推定したものです。現在でもデータセットの更新が続けられておりますので、最新のデータが必要になった場合には、再度ダウンロードしてみると良いでしょう。
もちろんRもインストールしなければなりません。最新版Rの入手法やインストール法については、こちら[外部リンク]を参考にして下さい。また、竹中明夫さんのサイト[外部リンク]にRのチュートリアルがあります。こちらも大変参考になります。
続いて、RでNetCDF形式を扱うための拡張ライブラリをインストールします。WindowsでRを立ち上げて、メニューの「パッケージ」--＞「インストール」とクリックし、ダウンロードするサーバを選んで、「ncdf」というライブラリを選択します。
以下に紹介するコードを一気に読み込んで使うには、例えば"netcdf_reader.txt"というテキストファイルにコードを保存し、これをワーキングディレクトリに置いて
source("netcdf_reader.txt")
と入力します。コードを部分的に使いたいときは、該当部分をコピー＆ペーストします。最初は、一行ごとにコピー＆ペーストして、Rの動作を確認なさると良いでしょう。
### パート１: netCDFの読み込み、変数名の確認、変数の取り込み
library(ncdf) #NetCDFライブラリの呼び出し
setwd('c:/R/') #作業ディレクトリの指定(環境に合わせて適切に設定して下さい)
dat <- open.ncdf("precip.mon.mean.nc") #NetCDFファイルのオープン
print(dat) #データの大きさを確認
dat$var #座標値や変数名を確認
x <- get.var.ncdf(dat,"lon") #軸lonの座標値をxに読み出す
y <- get.var.ncdf(dat,"lat") #軸latの座標値をyに読み出す
z <- get.var.ncdf(dat,"precip") #変数precipをzに読み出す
tmp1 <- z[73:144,,] #ここからの4行で、経度座標を変換
tmp2 <- z[ 1: 72,,] #
for(i in 1:72) z[i,,] <- tmp1[i,,] #
for(i in 73:144) z[i,,] <- tmp2[i-72,,] #
x = x - 180 #上に合わせて経度の表記を変更
まずはパート１のprint(dat)までを実行すると、データの次元や大きさが確認できます。すなわちデータがlon、lat、timeという３つの軸上に配置されている事、そして各軸に144、72、352個の点が打たれている事が分かります。省略されていますが、それぞれ経度、緯度、時間の意味でしょう。続いて、dat$varを実行すると、これら軸上の点に与えられた座標値、そして変数名($name直下に記述あり)がprecipであることが分かります。次の３行では、経度の座標値・緯度の座標値・データ本体を、それぞれ配列変数x、y、zに読み込みます。
すでにRにデータが取り込まれていますので、あとは自在に加工できます。例えばdim(z)とコマンドを打つと、zが144×72×352の配列サイズを有する三次元データであることが分かります。このデータは、西経180度を基点に東方向を正に座標値が与えられておりますので、パート１の残りの部分では、これを通常の経度の表現、すなわち東経0を起点にするように変更を行います。また経度の表記を、東経は正の値、西経は負の値と改めます。
パート２では、このデータファイルの最初の12ヶ月分の数字をテキスト形式で書き出します。できたファイルをテキストエディタで開くと、144個の数字が含まれる行が864行(72×12)できている事が分かります。
### パート２: データをテキストで書き出す
out <- file("precip.txt", "w") #ファイルを書き込みモードで開く
for (mon in 1:12) { #最初の12ヶ月
for (lat in 1:72 ) { #緯度
for (lon in 1:144) { #経度
if (is.na(z[lon,lat,mon])) {z[lon,lat,mon]=-9999} #欠測値を-9999で置換
if (lon<144) {writeLines(paste(z[lon,lat,mon]),out,sep=",")}
else {writeLines(paste(z[lon,lat,mon]),out,sep="\n")}
}}}
close(out) #ファイルを閉じる
Rのimageコマンドを使えば、このデータをプロットすることも出来ます。なお、以下のコードでは、プロットに世界地図を重ね書きするため「maps」という拡張ライブラリを使っていますので、最初に先ほどの要領でインストールしておいて下さい。
### パート３: データを地図に表示する library(maps) #mapライブラリの呼び出し frame() #画面のクリア br <- c(seq(from=0,to=5.5,by=0.5),100) #値の分割点を設定 col <- rainbow(14)[3:14] #各値の色設定(色の数=値の分割点+1) image(x,y=y[60:5],z[,60:5,1],breaks=br,col=col) #マッピング実行 map(add=T) #ここに地図を重ね書き
参考資料
http://www.esa-soilmoisture-cci.org/node/136#data-format
「How can I open the ECV SM data files in R?」
# パッケージのインストール
# これで上手くいかなければ、「パッケージ」→「CRANミラーサイトの設定」
install.packages("ncdf4")
# ライブラリの呼び出し
library(ncdf4)
# データファイルのオープン
dat <- nc_open("sample.nc")
# 開いたデータファイルを閉じる
nc_close(dat)
# 変数の読み込み
dat <- nc_open(file_name)
z <- ncvar_get(dat,"sm")
文責：佐藤永