pythonを使ったNetCDFの書き込みについて
サンプルを示しながら、解説する
サンプルプログラムのダウンロードはこちら.
sample.py(216)
このプログラムは、とあるプログラムから書き込みの部分だけを切り取ったものです。
そのままでは動きません。TDATAやIDATAに具体的な値を入れると動きます。
------
import netCDF4 from numpy import dtype def main():
下に記載しているのが、ハッシュへのハッシュを使い、変数名に対して、まとめてunitやlong_nameを定義しています。
# define netcdf attribute ####
varinfo = {'IWV':{'units':'mm','long_name':'Integrated Water Vapor','missing_value':-999}}
varinfo['LWP']= {'units':'mm','long_name':'Liquid Water Path','missing_value':-999}
varinfo['BRT']= {'units':'K','long_name':'Black Body Temperature','missing_value':-999}
varinfo['CMP.TPC']= {'units':'K','long_name':'Temperature Profile','missing_value':-999}
varinfo['HPC']= {'units':'g/kg','long_name':'Humidity Profile','missing_value':-999}
詳しい理由はわからないが、missing_valueとしては整数しか与えられないようです。-999.9とすると
cannot be safely cast というエラーが出ます。私は、-999としました。
変数のアトリビュートを追加しております。
# for surface observation ## #
varinfo['SP']= {'units':'hPa','long_name':'surface pressure','missing_value':-999}
varinfo['ST']= {'units':'K','long_name':'surface temperature','missing_value':-999}
varinfo['SRH']= {'units':r'%','long_name':'surface relative humidity','missing_value':-999}
varinfo['SU']= {'units':'m/s','long_name':'x-comp of surface wind','missing_value':-999}
varinfo['SV']= {'units':'m/s','long_name':'y-comp of surface wind','missing_value':-999}
varinfo['RR']= {'units':'mm/h','long_name':'rainrate','missing_value':-999}
NetCDFを書き込みとしてopenします。
# open NetCDF as output ########### nc = netCDF4.Dataset(outfile,"w",format="NETCDF4")
UNLIMITED の場合、Noneとします。そうじゃない場合、具体的な時間数を入れます
nc.createDimension("TIME",None)
ここでは、次元変数のTIMEに対して、変数を定義しています。CやFortranでは次元を定義するときに、
次元の型（doubleなのか、floatなのか）を定義しましたが、
Pythonでは、上のように、次元の定義は、次元数のみだけ。変数として定義するときに、型を定義します。
time = nc.createVariable('TIME', dtype('float64').char, ('TIME'))
time.long_name = "TIME"
time.units = "seconds since 1970-01-01 00:00:00"
nc.createVariableでは、第1引数に変数名、第2引数に型、第3引数にその変数を構成する次元 を記載します。
上の例では、次元変数のTIMEを変数として定義しているので、次元は１つでTIMEとなります。
なお、dtypeについては、以下のページを参照してください。
https://note.nkmk.me/python-numpy-dtype-astype/
timeという変数に対して、.を付けて、アトリビュートを足していきます。ここでは、long_nameとunitsを足しています。
次元変数なので、missing_valueは無いはずなので、足していません。
以下は、一次元配列の場合の書き込みと2次元配列の書き込みを説明します。
他のサイトでは、1次元配列の例はあるのですが、2次元配列の場合の例が全く無く、困りました。
if datatype == 1: ## 1 dimensional
var = nc.createVariable(vname, dtype('float32').char, ('TIME'))
var.long_name= varinfo[vname]['long_name']
var.units = varinfo[vname]['units']
var.missing_value = varinfo[vname]['missing_value']
上記に記載しているのは、変数名がvnameに入っていて、vnameが例えば、IWVの場合には、
プログラムの最初に書いた、long_nameやunitsやmissing_valueが入るようにしました。
vnameだけを変更するだけで、同じプログラムを使いまわせるからです。この変数はdoubleではなく、
floatの精度で十分だったため、float32を選びました。IWV(可降水量)は時間の1次元変数なので、
TIMEが次元として張っています。
次に、TDATAもデータ数だけUNIXTIMEを入れておき、その時刻に対応するデータをIDATAとして保存しておきます
（ほかのプログラムで。ここではTDATA,IDATAはもう収録されていると仮定しています）
time[:] = TDATA[:]
var[:] = IDATA[:]
これを、先ほどnc.createVariableで生成したオブジェクトに1次元配列として入れるだけ。最後にnc.closeすると書き込み完了。
次に2次元配列の場合の例を示します。他のサイト（日本語では他になかった。）では紹介していない事項です。
elif datatype == 2: ## 2 dimensional
まずは、2次元なので、既に定義した時間ともう一つの次元を定義します。ここでは、高さ(LEV)を定義します。
nc.createDimension("LEV",NX)
LEV = nc.createVariable('LEV', dtype('float64').char, ('LEV'))
LEV.long_name = "altitude"
LEV.units = "m"
次元変数にもlong_nameやunitsを記載します。今度の次元は、LEVとなります。LEVについては、NX個の層があるとしています。
時間が12個で高さが50個の次元があるとすると、収録しようとするデータは2次元配列で、12×50の配列になっている必要があります。
次に収録しようとする変数を定義します。違うのは、最後のTIME,LEVだけ。2次元で、それぞれの次元がTIME,LEVであるということ。
左側にある次元が上位の次元。右側が早く回る下位の次元。高さ方向にデータが50個入ると、次の時間の50個が入っているという順番。
var = nc.createVariable(vname, dtype('float32').char, ('TIME','LEV'))
var.long_name= varinfo[vname]['long_name']
var.units = varinfo[vname]['units']
var.missing_value = varinfo[vname]['missing_value']
createVariable以外は同じ。
最後に書き込み。これは、IDATAは2次元なのに対して、createVariableのオブジェクトのvarは1次元配列だということに注意！！
これが本家のマニュアルにも書いてなくて、試行錯誤でようやく成功しました。3次元配列でも、1次元としてNetCDFのオブジェクトには移し替えることが必要。
time[:] = TDATA[:]
LEV[:] = X_ARRAY[:]
var[:] = IDATA[:][:]
書き込みは、最後にcloseするだけ。
nc.close()
C言語やFortranとの違い
- 編集モード、書き込みモード というものがない！！
- したがって、enddefもしくはredefが無く、シンプル。
- ncidに対応して、上記ではnetcdf4.datasetから出てくるオブジェクト(nc)が対応している。
- varidに対応して、nc.createVariableから出てくるオブジェクト(var)が対応している。
- オブジェクト指向のプログラムで、ncやvarの下に.を付けていろんなことができるようになっている。
- 処理スピードだけが問題ではあるが、CやFortranに比べて、圧倒的にプログラミングが楽。これはいい。