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")
文責：佐藤永