memoring

GrADS形式で書かれたFortranで使うようなリトルエンディアンのヘッダなし4バイト浮動小数点バイナリ。

• 読み込み

f = open("ファイル名", "rb")
data = numpy.fromfile(f, dtype="<f4").reshape((tn,zn,yn,xn))

• 書き込み

data = np.array(data, dtype=np.float32)
data.tofile("ファイル名")

Fortranとは配列の順序が逆になることに注意。

ctlファイルがあるならpygradsも使える。

from grads.ganum import GaNum
ga = GaNum(Bin='grads -b')                             #X-windowはいならいのでバッチモードで
ga.open('z.1979.ctl')                                           #ga> open z.1979.ctl
z = ga.exp('hgt')                                                    #hgtという変数をnumpy ndarrayで読み込む

ga('set lev 1000 100')                                        #gradsコマンドが使える。

npy形式のファイルをoutputできるlibnpyというライブラリがあるらしい。こちらから。

program test_fnpy

    use fnpy
    implicit none

    integer, parameter :: DP = selected_real_kind(15)
    complex, parameter :: I = cmplx(0.0,1.0)

    real     :: a(4)   = [ 1, 2, 3, 4 ]
    real     :: b(3,2) = reshape([ 1, 2, 3, 4, 5, 6 ], [3, 2])
    real(DP) :: c(2,3) = reshape([1, 2, 3, 4, 5, 6], [2, 3])
    integer  :: n, d(2,5) = reshape([(n, n=1, 10)], [2,5])
    complex  :: e(3,1) = reshape([1+3*I, 2+2*I, 3+I], [3,1])
    complex(DP) :: f(2) = [ 1 + 2*I, 2 + I ]

    print *, "Creating files fa.npy, fb.npy, fc.npy, fd.npy"
    call save_single ("fa.npy", shape(a), a)
    call save_single ("fb.npy", shape(b), b)
    call save_double ("fc.npy", shape(c), c)
    call save_integer("fd.npy", shape(d), d)
    call save_complex_single("fe.npy", shape(e), e)
    call save_complex_double("ff.npy", shape(f), f)

end program test_fnpy

NumPyのみで使うならこの方が単純で速い。

#書き込み
numpy.save("test.npy", data)
#読み込み
data = numpy.load("test.npy")

savezを使うと複数の配列をひとつのファイルに保存できる。

#書き込み
lonarray = np.arange(0, 360, 2.5)
latarray = np.arange(-90, 90.5, 2.5)
np.savez('hogehoge.npz', lon=lonarray, lat=latarray)
#読み込み
data = np.load('hogehoge.npz')
lonarray = data['lon']
latarray = data['lat']

配列が辞書的に名前で保存される。

• 読み込み

data = np.loadtxt("asci.csv")

genfromtxtはloadtxtよりオプションが多く欠損も扱える。

data = np.genfromtxt("ascii.csv")                #読み込み data[行,列]となる。
data = np.genfromtxt("ascii.csv", dtype='|S5')   #文字列データの場合
data = np.genfromtxt("ascii.csv", skip_header=1) #ヘッダ1行を飛ばす

• 書き込み

np.savetxt("ascii.txt", data)
np.savetxt("ascii.txt", data, fmt="%12.6G")      #フォーマットの指定

SciPy.io.netcdfやPyNio、netcdf4-pythonなどがある。PyNioはgribやhdfも読める。netcdf4-pythonは時間を扱うクラスなどユーティリティが充実している。基本的な使い方はほぼ同じなのでnetcdf4-pythonを例にメモする。

• netcdf4-python

本家サイト
2011年Unidata Workshopでのスライド

どのライブラリもNetCDFのファイルオブジェクトを作成しその属性に次元や変数、値を格納する。ncdumpを使えばnetcdfの中身を見ることができる。

$ ncdump olr.day.mean.nc
netcdf olr.day.mean {
dimensions:
        lon = 144 ;
        lat = 73 ;
        time = UNLIMITED ; // (13791 currently)
        nmiss = 7 ;
variables:
        float lon(lon) ;
                lon:units = "degrees_east" ;
                lon:long_name = "Longitude" ;
                lon:actual_range = 0.f, 360.f ;
                lon:standard_name = "longitude" ;
                lon:axis = "X" ;
...
// global attributes:
                :title = "Daily Mean Interpolated OLR" ;
...
data:

 lon = 0, 2.5, 5, 7.5, 10, 12.5, 15, 17.5, 20, 22.5, 25, 27.5, 30, 32.5, 35, 
...

pythonの場合はdimensions、global atributesなどがファイルオブジェクトの属性となり、その属性に辞書型でデータが格納される。

netcdfクラスのファイルオブジェクトを用いて操作する。

from netCDF4
nc = netCDF4.Dataset("filename.nc", "r")
nc.close()

rだと読み込み専用、wだと書き込み専用。

• 読み込み

上記の例だとdimensions属性に入っている。

>>> print nc.dimensions
OrderedDict([(u'lon', <netCDF4.Dimension object at 0x227fb00>), (u'lat', <netCDF4.Dimension object at 0x227fb90>), (u'time', <netCDF4.Dimension object at 0x227fbd8>), (u'nmiss', <netCDF4.Dimension object at 0x227fc20>)])

例えば緯度lonの次元数は

>>> len(f.dimensions['lon'])
144

で取得できる。python4-netcdfライブラリではnetCDF4.Dimensionクラスのオブジェクトが格納されている。

• 書き込み

次元名と大きさで作成。

lon_dim = nc.createDimension('lon', 144)

• 書き込み

上記の例ではvariables属性に格納される。

>>> print nc.variables
OrderedDict([(u'lon', <netCDF4.Variable object at 0x2242ad0>), (u'lat', <netCDF4.Variable object at 0x2242a50>), (u'time', <netCDF4.Variable object at 0x2242bd0>), (u'olr', <netCDF4.Variable object at 0x2242c50>), (u'info', <netCDF4.Variable object at 0x2242cd0>)])

名前とnetCDF4.Variableオブジェクトが関連付けられている。numpy.ndarrayのように扱える。

>>> olr = nc.variables['olr']
>>>  olr.dimensions, olr.dtype, olr.shape, olr.ndim  #次元名、変数型、要素数、次元数
((u'time', u'lat', u'lon'), dtype('int16'), (13791, 73, 144), 4)
>>> olr.units   #ncdumpで名前=値で入っていた情報には属性名としてアクセス
u'W/m^2'
>>> print olr[::2. : , :]   #numpy.ndarrayのスライスと同じようにデータは読み取れる
[[[205.449996948 205.449996948 205.449996948 ..., 205.449996948
   205.449996948 205.449996948]
  [204.25 204.299987793 204.359985352 ..., 204.0 204.149993896 204.25]
  [204.5 204.699996948 204.899993896 ..., 203.949996948 204.149993896
   204.359985352]
  ..., 

netcdf4-pythonの場合はscale_factorやadd_offsetも加味した値がデータに入るみたい。scipyの場合はファイルの値がそのまま入るのでolr.scale_factorなどの属性名で取得してscaleやoffsetを調整しなければならない。

netcdf4-pythonにはdatetimeオブジェクトとnetcdfでの時間を変換するメソッドが用意されている。データとしては時間は実数型として入っていて扱いにくいのでdatetime型オブジェクトに変換すると便利。

num2date(times, units, calender='standard')	数値をdatetimeに変換する。unitsは単位と開始時間、calenderには'gregorian'や'365_day'など指定できる
date2num(dates, units, calender='standard')	datetimeを数値に変換する。最小単位は秒
date2index(dates, time)	netcdfのtime変数においてdatetimeに相当するインデックスを返す

• 例

>>> units = nc.variables['time'].units  #規約に従ったファイルなら時間の単位情報は文字列として刻み幅、開始日付が入る
>>> print units
hours since 1-1-1 00:00:0.0
>>>  time = nc.variables['time'][:]                            #実数で入っている時間データ
>>> dtime = netCDF4.num2date(time, units=units)   #数字を渡してもよいし、numpyのunicersal functionのように配列を渡してもよい
>>> dtime
array([1974-06-01 00:00:00, 1974-06-02 00:00:00, 1974-06-03 00:00:00, ...,
       2012-03-01 00:00:00, 2012-03-02 00:00:00, 2012-03-03 00:00:00], dtype=object)
>>> time = nc.variables['time']
>>> nt = netCDF.date2index(datetime(1979,1,1), time)
>>> olr = print nc.variables['olr'][nt, :, :]          #1979.1.1のolrが得られる

データが複数のファイルに分割されている場合、命名が規則的ならばまとめて読み込める。ただし個々に扱うより遅くなるらしい。

>>> nc = netCDF4.MFDataset('mftest*.nc)
>>> print f.variables['x'][:]
[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]

PyNioで読める。

gradsで読める形式ならpygradsを通して読むことができる。gradsでやればいいんだけどctlファイルをかくだけでgrads的に読めるのはとても直感的でわかりやすい。ctlファイルさえあればバイナリでもgribでもnetcdfでも読めるし、特定の範囲だけ読み出したい場合なども使いやすいし早い。

from grads.ganum import GaNum

ga = GaNum(Bin='grads -b')  #gradsを起動する。描画せずにデータのハブとしてのみ使うのでバッチモードにする。
ga('open test.ctl')                   #引数にコマンドをとってgradsを操作する。
ga('set lat 30 90')
ga('set t 1 31')

data = ga.exp('hgtprs')         #gradsで指定した範囲のデータがnumpy配列でexportされる。

気象関係の多くのフォーマットに対応している。scipy.io.netcdfとよく似ており、netcdf、grib、hdfも扱えるのでこちらの方がいいだろう。NCAR command languageがインストールされている必要がある。

import numpy as np
import Nio

f = Nio.open_file("test.grib", "r")   #ファイル形式は拡張子で判断する。拡張子がない場合はformat="grib"などのオプションをつける。

f.close()

以下のファイルオブジェクトの属性に辞書の要領でアクセスする。

dimensions	次元名と次元を辞書で格納
variables	変数名とNioVariableオブジェクトを辞書で格納
attributes	global attributeがattribute名をキーにして辞書で格納

• 例

>>> print f.dimensions
{'lat': 73, 'lon': 144, 'time': 768}
#lat次元の要素数
>>> print f.dimensions['lat']
73
#global attributeのConventionsの中身
>>> print f.attributes['Convntions']

NioVariableオブジェクトに変数情報が格納されている。Nio.NioVariableは以下の属性を持つ。

rank	次元
shape	次元の要素数がタプルで格納
dimensions	次元名がタプルで格納
attributes	変数のattributeをキーにして辞書で格納

変数の値をarrayで取得するには以下のようにする。

data = f.variables['slp'][:]
#または
data = f.variavles['slp'].get_value()

求めたい場所、時間が何番目に当たるかを知っているなら（GrADS等で調べる）スライスで求めれば良い。要素のインデックスを知らなくても範囲を直接指定してスライスを得ることができる。

import Nio
f = Nio.open_file("olr.day.mean.nc", "r")
olr = f.variables['olr']
#先頭の時間要素、緯度0度、経度0-20E
input = olr['time|i0 lat|0 lon|0:20']

次元名|範囲をスペースで区切った文字列で指定する。次元名を省略する場合はデータの並びの次元となる。以下のように様々な指定ができる。

• 先頭の時間、緯度15N-15S、経度0-20E

input = olr['time|i0 lat|15:-15 lon|0:20']

範囲はデータに入っている順に指定する。このデータの場合緯度次元が90Nから南に変数が入っているので-15:15とするとエラーになる。

• 先頭の時間、緯度15N-15S、経度0、20E

input = olr['time|i0 lat|15:-15 lon|0,20']

離散的に指定する場合はカンマで区切る。この場合データの順に指定する必要はない。指定した順で格納されてデータが返ってくる。

• 先頭の時間、緯度15N-15S、経度0-20Eまで5度きざみ

input = olr['time|i0 lat|15:-15 lon|0:20:5']

• 先頭の時間、緯度15N-15S、経度0-20Eまで3点きざみ

input = olr['time|i0 lat|15:-15 lon|0:20:i3']

iを数字の前につけると要素数（格子）単位を表す。

• 先頭の時間、緯度15N-15S、経度0-20Eまで1.25度ごとに内挿

input = olr['time|i0 lat|15:-15 lon|0:20:1.25i']

数字の後ろにiをつけると内挿する。

• 先頭の時間、緯度15S-15N、経度0-20E

input = olr['time|i0 lat|15:-15:i-1 lon|0:20']

負の値は逆順を表す（スライスと同じ）。iをつけてi-1としていることに留意。つけないと1度単位を表してしまう。