Python-NetCDF解析 (データ解析編)

NetCDFの読み込み編

  importするモジュール

import netCDF4

と先頭に記載する

  コマンドの引数からデータを読むために

import sys

を先頭に記載。

if len(sys.argv) < 5:                                                          
   print ("ARGV error: infile var t k")                                       
   exit(1)    

などとして、引数の数が5個以下なら、入力ファイル名や変数名、時間断面、層番号を入力させる。

infile = sys.argv[1]
var    = sys.argv[2]                                                           
tt     = int(sys.argv[3])                                                      
kk     = int(sys.argv[4])                                                      
print ("infile=", infile)                                                      
print ("var=", var)                                                            
print ("time=",tt)                                                             
print ("layer=",kk)    

として、引数をコマンドラインから取れる。そして、printで確認

  読み込み用としてopen

nc = netCDF4.Dataset(infile,'r')

  netcdfの最上位のクラスの中身を確認する

import netCDF4
import inspect

とした上で、

nc = netCDF4.Dataset(infile,"r")
for x in inspect.getmembers(nc):
  print (x)

とすると、
　nc.(個々に関数や変数を設定でき、いろいろな情報が得られる)
の全種類を閲覧することができる。これでnetcdfのデータの読み方が分かるはずです。

以下、主要な種類を上げておきます。

種類	意味
dimensions	次元に関する情報
variable	変数に関する情報
close	関数であり、netcdfを閉じる

試しに、

nc = netCDF4.Dataset(infile,"r")
print (nc.dimensions) 
print (nc.variable)
print (nc.ndims)

とすることで、情報を得られることを確かめてください。

  有る変数の情報を取得する

上記のクラスの中身の見方が分かると、マニュアルが無くても、pythonで
netcdfのデータの読み方がほぼ分かることになる。
netcdfのクラスは階層上になっており、nc.variable['SU']のように変数SUの情報だけを取得できるようになる。
例えば、変数SUが何次元で、どの次元を持つかを調べたいときには、

print (nc.variable['SU'].ndim)
print (nc.variable['SU'].dimensions)

とすれば、いい。

nc.variable[変数名].知りたい情報

のようになっている。知りたい情報をどうやって指定したらいいか分からない場合は、

a = nc.variable['SU']
for x in inspect.getmembers(a):
  print (x)

とやれば、知りたい情報の一覧が表示されるので、その情報を.の後に付ければいい。

  様々なサンプル集

ある変数Xが存在するかの確認

a = nc.variable
if "X" in a:
  print ("X is included")
else :
  print ("X is not included")

有る変数Xの次元数と次元の種類を知りたい

ndim = nc.variable['X'].ndim
dims = nc.variable["X"].dimensions
print ("ndims=",ndim)
print ("dims=",dims)

有る次元変数Yが存在するかの確認

a = nc.dimension
if "Y" in a

LAT-LON座標系だと確定している場合の格子数を取得したい

lon = nc.dimension['LON']
XNUM = len(lon)
lat = nc.dimension['LAT']
YNUM = len(lat)

変数Uが4次元（TIME,LEV,LAT,LON）で格納されており、ある時刻tの3次元分布を取得したい

U3D = nc.variable['U'][t][:][:][:]

変数Uが4次元（TIME,LEV,LAT,LON）で格納されており、ある時刻tと層番号kの2次元分布を取得したい

U2D = nc.variable['U'][t][k][:][:]

2バイトで圧縮されているデータを4バイトfloatにして取得したい

TMP = nc.variable['U'][t][k][:][:]
add_offset = nc.variable['U'].add_offset
scale_factor = nc.variable['U'].scale_factor
U2D = TMP * scale_factor + add_offset

変数Uが2バイト圧縮されているかどうかを調べたい

dtype = nc.variable['U'].dtype
if dtype == 'int16':
   print ("short compress")
elif dtype == 'float32':
   print ("float: no compress") 

そもそも、pythonでnetCDFを読むと自動的にfloatになっている！

どうやらユーザーはadd_offsetやscale_factor,missing_valueなどは読まなくていいようです。
自動的にshortから正しい実数に変換され、missing_valueの場所には、その配列の要素は、未定義となるようです。

 サンプルプログラム

main.pyとread_NC.pyの2つのファイルに分けて、netCDFを読み込むプログラムをすべて
read_NC.pyにまとめた。二つのファイルは、同じディレクトリにあることを想定しています。

main.py(324)
read_NC.py(309)

ソースコード

main.py

import sys
from read_NC import read_NC
def main():
   if len(sys.argv) < 3:
       print ("ARGV error: infile var [t] [k]")
       exit(1)
   infile     = sys.argv[1]
   varname    = sys.argv[2]
   if len(sys.argv) > 4:
       t = sys.argv[3]
   else :
       t = 0
   if len(sys.argv) > 5:
       k = sys.argv[4]
   else :
       k = 0
   DATA = read_NC(infile,varname,t,k)
   print (DATA)
####################################
if __name__ == "__main__":
   main()

read_NC.py

import netCDF4
def check_NC_type(infile,varname):
   nc = netCDF4.Dataset(infile,'r')   
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)      
   dim = nc.variables[varname].dimensions
   nc.close()
   dim_check = [0,0,0,0]
   if "TIME" in dim:
       dim_check[0] = 1
   else:
       dim_check[0] = 0
   if "LEV" in dim or "ALT" in dim:
       dim_check[1] = 1
   else:
       dim_check[1] = 0
   if "LAT" in dim:
       dim_check[2] = 1
   else:
       dim_check[2] = 0
   if "LON" in dim:
       dim_check[3] = 1
   else:
       dim_check[3] = 0
   return dim_check
def read_NC_T2D(infile,varname,t):
   nc = netCDF4.Dataset(infile,'r')    
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   DATA2D = nc.variables[varname][t][:][:]
   nc.close()
   return DATA2D
def read_NC_T3D(infile,varname,t,k):
   nc = netCDF4.Dataset(infile,'r')   
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   DATA2D = nc.variables[varname][t][k][:][:]
   nc.close()
   return DATA2D
def read_NC_2D(infile,varname):
   nc = netCDF4.Dataset(infile,'r')  
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   DATA2D = nc.variables[varname][:][:]
   nc.close()
   return DATA2D
def read_NC_1D(infile,varname):
   nc = netCDF4.Dataset(infile,'r')  
   if not varname in nc.variables:
       print ("not found ", varname)
       exit(1)
   DATA1D = nc.variables[varname][:]
   nc.close()
   return DATA1D
def read_NC(infile, varname,t=0,k=0):
   type = check_NC_type(infile,varname)
   ndim = sum(type)
   if type   == [1,0,1,1] :
       print ("2D evolution")
       DATA = read_NC_T2D(infile,varname,t)
   elif type == [1,1,1,1] :
       print ("3D evolution")
       DATA = read_NC_T3D(infile,varname,t,k)
   elif type == [0,1,1,1] :
       print ("3D static")
       DATA = read_NC_T2D(infile,varname,k)
   elif type == [0,0,1,1] :
       print ("2D static")
       DATA = read_NC_2D(infile,varname)
   elif ndim == 1:
       print ("1dimensional")
       DATA = read_NC_1D(infile,varname)
   return DATA

解説

main.py は引数として、入力ファイル、変数、（時間指定）（層番号指定）
として、最低でも2つの引数を、必要であれば、時間と高さの指定を行い、4つまでの引数を受け付けます。
main.py は、read_NC.pyを呼び出して、DATAという2次元配列を取得することを目的とします。
この2次元配列は、X-Y平面を想定しています。抜き出したい変数がどんな次元でも抜き出せるように
read_NC.pyは工夫しています。ただし、次元変数としては、TIME,ALTもしくはLEV, LAT, LONという
次元変数名であることを想定しています。この次元変数とマッチするかどうかを確かめて、read_NC.pyの
中で、呼び出す関数を選択しています。read_NC.py の引数は、入力ファイルと変数名を
必須として、時間と層番号については任意選択できるようにしています。引数で
t=0,k=0と設定しているのは、引数設定が無い場合には、デフォルトで0が割当てられる
デフォルト引数というものを設定しています。

main.pyの一番下にある文は、端末から入力された場合には＿＿main＿＿という名前空間が
割り当てられるので、その場合、mainを実施するという意味である。
各サブルーチンもこうした＿＿main＿＿という名前空間を持たせた場合の実験ができる。