データの読み書き

データ読み書きの基本はNetCDFを参照。

データの読み書き

NetCDF†

NCLのファイル読み書きをNetCDFを基本に説明する。NCLはNetCDFの変数を読み込んだ時に自動的に格子情報やattributionを読み込んでくれる。NCLの利点のひとつである。

↑

読む†

まずは開きたいファイルをNCLに伝える作業を行う（Fortranでいうところのopen）。

addfile

f = addfile(filename, X)

〔入力変数〕

filename

読み込むファイル名（操作したいファイル名）ファイル名には拡張子を必ず付けること（実際のファイルに付いていない場合でも付ける）。"ファイル名" もしくはファイル名が入ったcharacter型変数で指定する。絶対パスでも相対パスでも構わない。

X

X = "r" とすると読み込み "w" とすると書き込み "c"とすると新しいファイル形成

〔出力変数〕

f

ファイルを示す変数(file型)。

対応するファイル形式
ファイル形式	拡張子	X
NetCDF	".nc", ".cdf", ".netcdf"	"r", "w", "c"
GRIB1,2	".gr", ".gr1", ".grb", ".gr2" など	"r" のみ
HDF4,5	".hdf", ".hd", ".h5"	"r", "w", "c"

複数のファイルを一度に開きたいときはaddfilesを用いる。

addfiles

f = addfiles(filenames, X)

〔入力変数〕

filenames

読み込むファイル名が入った1次元配列。自分で作ってもよいし，systemfuncなどを使って得てもよい。

X

X = "r" とすると読み込み "w" とすると書き込み "c"とすると新しいファイル形成

〔出力変数〕

f

ファイルを示す変数(list)

無事にファイルを開けたら変数を取り出す。ファイル内の変数を指定するさいには->という記号を使う。

var = f->X

var = f[:]->X

var = f->$X$

など。逆に変数varの値をXとしてファイルfに書き込みたい時には，

f->X = var

のようにする。

NetCDFでの出力のやり方は4byteバイナリ：例2も参照。

f

ファイルを示す変数。前述のaddfile等で開いたもの。fが複数のファイルからなる場合，すべてから変数を取り出したければf[:]とする。fのうち一つのファイルだけから取り出したければf[i]のようにできる。

X

取り出す変数名。ファイル内の変数名を直接書いても良いし，変数名が文字列として入ったNCL上の変数でも良い。後者については，例えばファイル内の変数名がNCL上の変数Xに入っている場合，f->$X$のように$で囲んで記述する。Xは配列でもよく，f->$X(i)$のように記述できる。

var

NCLスクリプト上での変数

addfilesで開いたファイルから変数を読む場合，すべてのファイルで変数名が共通で，最も左の次元以外のサイズが同じでなくてはならない。
変数は，defaultでは最も左の次元が結合された配列として読み込まれる。つまり，1ファイルに[12]×[180]×[360]のデータが入っており，これを30ファイルから同時に読み込むと，[12*30]×[180]×[360]となる。タイムステップごとに別のファイルにlat×lonとして入っているのをtime×lat×lonとして読みたいときや，アンサンブルメンバー別のファイルを読み込んで左端にアンサンブル次元を増やしたいときなど，複数ファイルを別次元として扱いたい場合，読み込む前にListSetTypeで設定変更をすればよい。

ListSetType

ListSetType(f, X)

f

file属性の1次元配列。addfilesの出力変数

X

X = "join" とするとファイル数の長さを持つ次元が左端に加わる。 X = "cat" とすると最も左の次元を結合して読み込む。

また，ファイル毎に異なるoffsetやscaleの値をもつファイルを一度に読み込む場合には，これらの値を取っておかないと，正しい値を復元できないので注意。
なお，取り出す変数の変数名はncdumpやncl_filedumpなどのコマンドで確認してもよいし，getfilevarnamesで変数名を読み込んでもよい。以下に，ファイル内の変数に関する情報を取得する関数をいくつか紹介する。

getfilevarnames
変数名を得る

getfilevartypes
変数の型を得る

getfilevardimnames
変数の次元名を得る

getfilevardimsizes
変数の次元の大きさを得る

例1：例として，precip.mon.mean.ncというNetCDFファイルを読むとしよう。ファイル内の変数等の情報は，例えばコマンドライン上で，

$ ncdump -h precip.mon.mean.nc

とすると

netcdf precip.mon.mean {
dimensions:
        lon = 144 ;
        lat = 72 ;
        time = UNLIMITED ; // (449 currently)
variables:
      ... (途中省略) ...
        float precip(time, lat, lon) ;
                precip:long_name = "Average Monthly Rate of Precipitation" ;
                precip:valid_range = 0.f, 70.f ;
                precip:units = "mm/day" ;
                precip:add_offset = 0.f ;
                precip:scale_factor = 1.f ;
                precip:missing_value = -9.96921e+36f ;
      ... (以下省略)

のように出る。このとき，NCLスクリプト

in = addfile("precip.mon.mean.nc","r")  ; ファイルを開く
x  = in->precip                         ; ファイル内の変数precipをNCL変数xに入れる

printVarSummary(x)                      ; 変数xについての情報を表示

を実行すると，xにprecipのデータが入って，

Variable: precip
Type: float
Total Size: 18620928 bytes
            4655232 values
Number of Dimensions: 3
Dimensions and sizes:   [time | 449] x [lat | 72] x [lon | 144]
Coordinates:
            time: [1569072..1896312]
            lat: [88.75..-88.75]
            lon: [1.25..358.75]
Number Of Attributes: 15
  long_name :   Average Monthly Rate of Precipitation
  valid_range : (  0, 70 )
  units :       mm/day
  add_offset :   0
  scale_factor :         1
  missing_value :       -9.96921e+36
 ... (途中省略) ...
  _FillValue :  -9.96921e+36

のように出る。

↑

書く†

NetCDFでの書き出しも，読み込みと同様にすることでできる。つまり，まずは書き出したいファイルをaddfileを用いて指定する。
例えば，ファイル名をoutfileとすると，新しいファイルの作成にはオプション"c"を指定することに注意して，

out = addfile(outfile,"c")

とする。スクリプト内の変数uをuという変数名で，変数vをvという変数名で書き出したいときには，

out->u = u
out->v = v

とすれば書き込みできるし，スクリプト内の変数名とファイルでの変数名を変えることも，もちろん可能である。

out->ucomp = u
out->vcomp = v

以上が最も簡易的なやり方である。しかし，ファイルや変数，座標に対して細かい指定をすることが必要な場合があり，その場合には以下の例2のような方法を使う。また，こちらのやり方のほうが速いらしい。

例2：時間time[長さNT]，緯度lat[長さNY]，経度lon[長さNX]の座標をもつ3次元配列u，vがあるとする。これをwind.ncというファイルに書き込む。まずは，ファイルを開く。

system("rm -f wind.nc")           ;; すでにそのファイルがあれば消す
out = addfile("wind.nc","c")

ファイルの各変数が持っている座標変数の名前やサイズの配列を作っておく。

dimNames = (/"time","lat","lon"/)
dimSizes = (/NT,NY,NX/)
dimUnlim = (/True,False,False/)   ;; timeには配列に無限長(unlimited)を持たせている

setfileoptionでNetCDFの書き込みオプション"DefineMode"をTrueに指定する。

setfileoption(out,"DefineMode",True)

ファイルそのものにつく属性等を指定する。

fatts = True
fatts@creation_date = systemfunc("date")  ;; dateコマンドの戻り値
fileattdef(out,fatts)

ファイルに含まれる座標変数を指定する。

filedimdef(out,dimNames,dimSizes,dimUnlim)

ファイルに含まれる全ての変数を指定する。

filevardef(out,"time",typeof(time),getvardims(time))
filevardef(out,"lat",typeof(lat),getvardims(lat))
filevardef(out,"lon",typeof(lon),getvardims(lon))
filevardef(out,"ucomp",typeof(u),getvardims(u))   ;; 変数uをucompという名前で書くことにする
filevardef(out,"vcomp",typeof(v),getvardims(v))

変数のattributeを指定する。

filevarattdef(out,"time",time)
filevarattdef(out,"lat",lat)
filevarattdef(out,"lon",lon)
filevarattdef(out,"ucomp",u)
filevarattdef(out,"vcomp",v)

値を書き込む(値のみの代入を行う)。

out->time  = (/time/)
out->lat   = (/lat/)
out->lon   = (/lon/)
out->ucomp = (/u/)
out->vcomp = (/v/)

なお，デフォルトでは2GB以上のファイルを書き出すことができないと思われる。大きなNetCDFファイルを作成したい場合には，あらかじめ，setfileoptionでNetCDFの"Format"オプションを"LargeFile"に指定する。

setfileoption("nc","Format","LargeFile")

↑

GRIB†

基本的にNetCDFと同じ。ただし読み込みしかできない。変数名はコマンドラインからncl_filedumpで確認してもよいし，getfilevarnamesで知ることもできる。

↑

4byteバイナリ(GrADSの形式)†

setfileoption
開くファイルの種類、データの形式、エンディアンを指定する。 GrADSのデフォルトは、データ形式はダイレクトアクセスで、エンディアンはPCに依存となっている。（nclもデフォルトのエンディアンはPC依存）

fbindirread
ダイレクトアクセス(access='direct')の場合はこの関数で読む

fbinrecread
シーケンシャル(access='sequential')の場合はこちらで読む GrADSの場合levごとに書き出されていると思うので注意。

例1： CTLファイルが

DSET ^hoge.grd
OPTIONS BIG_ENDIAN
XDEF 360 LINEAR 0 1.
YDEF 181 LINEAR -90 1.
ZDEF 1   LEVELS 0
TDEF 12 LINEAR 01JAN2000 1mn 
VARS 1
hoge 0 99 
ENDVARS

のとき1ステップ目(01JAN2000)を持って来たい時は

setfileoption("bin", "ReadByteOrder", "BigEndian")
hoge= fbindirread("hoge.grd",0,(/1,181,360/),"float")

のようにする。

Attributesで未定義値(undef)を設定するとよい

例2：例1と同じファイルからすべてのデータを読み，NetCDFファイルで書き出す。

xdef   = 360      ; 経度方向のサイズ
ydef   = 181      ; 緯度方向のサイズ
tdef   = 12       ; 時間方向のサイズ
msgval = -1.e21   ; 未定義値が-1.e21だとする

;;; varの座標情報となる緯度，経度，時間の配列をつくる 
lat        = fspan(-90,90,ydef)
lat@units  = "degree_N"
lon        = fspan(0,359,xdef)
lon@units  = "degree_E"
time       = ispan(0,tdef-1,1)
time@units = "months since 2000-01-01 00:00:00"
;;; 座標変数自体にも座標をつけたほうがお行儀が良い（別に必要ないことも多い）
lat!0      = "lat"
lat&lat    = lat
lon!0      = "lon"
lon&lon    = lon
time!0     = "time"
time&time  = time

;;; 読む
setfileoption("bin", "ReadByteOrder", "BigEndian")
var  = fbindirread("./hoge.grd",0,-1,"float")
var := reshape(var,(/tdef,ydef,xdef/))

;;; varに座標情報などを付加する
var!0 = "time"
var!1 = "lat"
var!2 = "lon"
var&time  = time
var&lat   = lat
var&lon   = lon
var@_FillValue = msgval 
var@units = "degC"
var@longname = "temperature"

;;; NetCDFで出力
system("rm -f ./hoge.nc")
out = addfile("./hoge.nc","c")
att = True
att@title = "temperature data hogehoge"
att@source_file   = "hoge.grd"
att@creation_date = systemfunc("date")
fileattdef(out,att)             ; ファイル全体についての情報を付加
filedimdef(out,"time",-1,True)  ; 時間次元をunlimitedに
out->hoge = var

↑

ASCIIデータ†

ASCIIデータを読むにはasciireadを使う。

asciiread

x = asciiread(filename, dim, type)

〔入力変数〕

filename

読み込むファイル名

dim

読み込むデータの次元を示す1次元配列またはスカラー。例えば，ファイル内のデータを5×3の配列として読みたければ (/5,3/) とする。

type

変数の型を示す文字列

〔出力変数〕

x

読み込んだ値。typeで指定した型のdimで指定した大きさをもつ配列。

読みたいASCIIデータのファイルにヘッダーやフッターの行が含まれていて，それを除外して値を読みたい場合，readAsciiTableが使えるかもしれない。

ASCIIデータを書くにはasciiwriteを使う。

asciiwrite

asciiwrite(filename, x)

〔入力変数〕

filename

書き込むファイル名

x

出力したい値

この関数では，xがいかなる形の配列であっても，値は各行に一つずつ書き込まれる。表のようなフォーマットで出力したいときにはwrite_tableが使えるかもしれない。

↑

CSV†

単に数値のみから構成されるCSVを読む方法は基本的にはASCIIと同じで，asciireadなどを使えばよい。
文字を含むCSVを読む場合には，まずasciireadなどを用いて1行ずつ文字列として読み込んだあとで，文字列を変換する関数を用いる。
空白が含まれていないCSVファイルの場合は，1行ごとに読んだ文字列をstr_splitでデリミタを","として分割すればよい(例1)。
しかし，一般的にCSVファイルは空白を含むことも多い。このような場合は，str_split_csvを用いるのが便利である(例2)。

例1： CSVファイル(test1.csv)

2015,1,10.0,A
2015,2,10.5,B
2015,3,11.0,C
2015,4,11.5,D
2015,5,12.0,E

をstr_splitを用いて読むスクリプト例

in    = asciiread("test1.csv",-1,"string")
ROW   = dimsizes(in)
dum   = str_split(in(0),",")
COL   = dimsizes(dum)
var   = new((/ROW,COL/),"string")
var(0,:) = dum
do i = 1, ROW-1
  var(i,:) = str_split(in(i),",")
end do
delete([/in,dum/])

X = tofloat(var(:,2))
Y = var(:,3) 

print(X)
print(Y)

この結果として，

Variable: X
Type: float
Total Size: 20 bytes
             5 values
Number of Dimensions: 1
Dimensions and sizes:   [5]
Coordinates:
Number Of Attributes: 1
  _FillValue :  9.96921e+36
(0)     10
(1)     10.5
(2)     11
(3)     11.5
(4)     12

Variable: Y
Type: string
Total Size: 40 bytes
             5 values
Number of Dimensions: 1
Dimensions and sizes:   [5]
Coordinates:
Number Of Attributes: 1
  _FillValue :  missing
(0)     A
(1)     B
(2)     C
(3)     D
(4)     E

を得る。

例2： CSVファイル(test2.csv)

2015,1,10.0,A
2015,2,,B
2015,3,11.0,C
2015,4,11.5,
2015,5,12.0,E

をstr_split_csvを用いて読むスクリプト例

var  = asciiread("test2.csv",-1,"string")
var := str_split_csv(var,",",0)

X = tofloat(var(:,2))
Y = var(:,3) 

print(X)
print(Y)