はじめに
Pythonで気象データを扱う際にnetCDFファイルを読み込むことがあります。以前はnetCDF4というライブラリを使って読み込むことが多かったのですが、現在ではxarrayが使われることも多いようです。
netCDF4ライブラリ（以下netCDF4と言う場合はファイルフォーマットではなくライブラリを指す）とxarrayとで文法は似ている点が多いのですが、一部取り扱いに違う点があります。そこで気づいた範囲について差異をまとめました。
実行環境については、Macでmatplotlib+cartopyの環境構築の記事を参照下さい。
netCDF4とxarrayの比較
前準備
netCDFファイルを手元に用意します。例えば、Z__C_RJTD_20210306180000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.binというGRIB形式のMSMのGPVファイルがあるとして、これをwgrib2を使って適当にnetCDFファイルに変換します。ここでは500hPa面のジオポテンシャル高度だけ切り出していますが、他の変数も一緒に変換しても構いません。
$ GRIB=Z__C_RJTD_20210306180000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin
$ NC=Z__C_RJTD_20210306180000_MSM_GPV_Rjp_L-pall_FH00-15.nc
$ wgrib2 $GRIB -match ":HGT:" -match ":500 mb:" -netcdf $NC
1.55:0:d=2021030618:HGT:500 mb:anl:
1.147:0:d=2021030618:HGT:500 mb:3 hour fcst:
1.239:0:d=2021030618:HGT:500 mb:6 hour fcst:
1.331:0:d=2021030618:HGT:500 mb:9 hour fcst:
1.423:0:d=2021030618:HGT:500 mb:12 hour fcst:
1.515:0:d=2021030618:HGT:500 mb:15 hour fcst:
$ ncdump -h $NC
# 以下、netCDFのヘッダ情報が表示される(略)
ファイルオープン
datasetをopenする部分が若干異なりますが、やることは同じです。以下必要に応じてコメントで変数の型(type)を書いておきます。
import netCDF4
import xarray as xr
ncfile = "Z__C_RJTD_20210306180000_MSM_GPV_Rjp_L-pall_FH00-15.nc"
ds_nc = netCDF4.Dataset(ncfile) # <class 'netCDF4._netCDF4.Dataset'>
ds_xr = xr.open_dataset(ncfile) # <class 'xarray.core.dataset.Dataset'>
変数の属性
変数にはds.variables[変数名]の形でアクセスできることは共通ですが、属性(attribute)のアクセス方法が異なります。
hgt_nc = ds_nc.variables["HGT_500mb"] # <class 'netCDF4._netCDF4.Variable'>
hgt_xr = ds_xr.variables["HGT_500mb"] # <class 'xarray.core.variable.Variable'>
# 属性へのアクセスは 変数.属性名
print(hgt_nc.units) # 'm'
print(hgt_nc._FillValue) # 9.999e+20
# 属性へのアクセスは 変数.attrs[属性名]
print(hgt_xr.attrs["units"]) # m
print(hgt_xr.attrs["_FillValue"]) # KeyErrorになる
netCDF4ライブラリの場合は元のnetCDFに書かれている属性をすべてそのまま読み出す事ができますが、xarrayの場合は値の配列(.values)が情報を持っていて不要な属性(e.g. _FillValue)については.attrsでアクセスできないようです。
値の配列
numpyのarrayとしてデータを取り出したい場合があります。その方法はnetCDF4とxarrayで異なります。
# 値をnumpyのarray(MaskedArray)で取得するには、添字を指定
values_nc = hgt_nc[:,:,:] # <class 'numpy.ma.core.MaskedArray'>
# 多次元配列の場合、コロンを1個にしても↑と同じ
values_nc = hgt_nc[:] # <class 'numpy.ma.core.MaskedArray'>
# .valuesで値の配列をnumpyのarrayで取得できる
valuse_hgt_xr.values # <class 'numpy.ndarray'>
# 添字を指定してスライスできるが、xarray.Variableのまま
hgt_xr[0,:,:] # <xarray.Variable (latitude: 253, longitude: 241)>
このようにxarrayを使う場合、xarray.Variablesのままでスライスすることができます。余談ですが、maplotlibのでcoutour図を描画する場合、2次元配列であれば、わざわざnumpy.ndarrayに変換せずともxarray.Variablesのままcontourに渡すことができます。
import matplotlib.pyplot as plt
plt.contour(hgt_xr[0,:,:]) # xarray.Variablesの2次元配列なのでOK
plt.contour(hgt_xr.values[0,:,:]) # numpy.ndarrayの2次元配列なのでOK
欠損値の扱い
numpyの配列として取り出す場合、netCDF4ではMaskedArray、xarrayではndarrayと型が異なります。それでは欠損値があるデータを読み込んだ場合には中身はどうなるでしょうか？ 試しに陸では値を持たない日本沿岸海況監視予測システムGPVの海洋最上層の水温(≃海面水温)のデータを読み込んでみます。
ncfile = "Z__C_RJTD_20210109000000_OCN_GPV_Rjp_Gll2km_Lz1-1000_Psbs_FD01.nc"
ds_nc = netCDF4.Dataset(ncfile)
ds_xr = xr.open_dataset(ncfile)
sst_nc = ds_nc.variables["WTMPSS_1mbelowsealevel"]
sst_xr = ds_xr.variables["WTMPSS_1mbelowsealevel"]
print(type(sst_nc[:])) # <class 'numpy.ma.core.MaskedArray'>
print(sst_nc[:])
# [[[297.1139221191406 297.1451721191406 297.2076721191406 ...
# 300.9264221191406 300.9264221191406 300.9576721191406]
# [297.1139221191406 297.1451721191406 297.2076721191406 ...
# 300.9264221191406 300.9264221191406 300.9264221191406]
# [297.1451721191406 297.1764221191406 297.2076721191406 ...
# 300.9264221191406 300.9264221191406 300.9264221191406]
# ...
# [-- -- -- ... 275.3951721191406 275.3951721191406 275.4264221191406]
# [-- -- -- ... 275.3639221191406 275.3639221191406 275.3951721191406]
# [-- -- -- ... 275.3014221191406 275.3326721191406 275.3326721191406]]]
print(type(sst_xr.values)) # <class 'numpy.ndarray'>
print(sst_xr.values)
# [[[297.11392 297.14517 297.20767 ... 300.92642 300.92642 300.95767]
# [297.11392 297.14517 297.20767 ... 300.92642 300.92642 300.92642]
# [297.14517 297.17642 297.20767 ... 300.92642 300.92642 300.92642]
# ...
# [ nan nan nan ... 275.39517 275.39517 275.42642]
# [ nan nan nan ... 275.36392 275.36392 275.39517]
# [ nan nan nan ... 275.30142 275.33267 275.33267]]]
netCDF4の場合は欠損値の表現にMaskedArrayのマスクの機能が使われているのに対し、xarrayの場合は欠損値の部分にnp.nanが埋められたndarrayが返ってきます。
ちなみにMaskedArrayとnp.nanが入ったndarrayは以下の方法で相互変換可能です。
ndarray_with_nan = masked_array.filled(fill_value=np.nan)
masked_array = np.ma.masked_array(ndarray_with_nan, np.isnan(ndarray_with_nan))
時刻の取得
netCDF4を使う場合、時刻はnetCDFで格納された値そのままになっています。
time_nc = ds_nc.variables["time"]
print(time_nc.units)
# seconds since 1970-01-01 00:00:00.0 0:00
print(time_nc[:])
# [1.6150536e+09 1.6150644e+09 1.6150752e+09 1.6150860e+09 1.6150968e+09
# 1.6151076e+09]
したがって、datetime.datatimeのリストに変換したい場合、少々アドホックですが例えば次のようにします。
t0 = datetime.strptime(time_nc.units[14:33], "%Y-%m-%d %H:%M:%S")
time_list = [t0 + timedelta(seconds=seconds) for seconds in time_nc[:]]
print(time_list)
# [datetime.datetime(2021, 3, 6, 18, 0), datetime.datetime(2021, 3, 6, 21, 0), datetime.datetime(2021, 3, 7, 0, 0), datetime.datetime(2021, 3, 7, 3, 0), datetime.datetime(2021, 3, 7, 6, 0), datetime.datetime(2021, 3, 7, 9, 0)]
一方、xarrayの場合、numpy.datetime64の配列として読み込んでくれます。unitsは不要なためか、.attrs["units"]でアクセスすることはできません。
time_xr = ds_xr.variables["time"]
for time in time_xr.values:
print(time, type(time)) # 2021-03-06T18:00:00.000000000 <class 'numpy.datetime64'>
time_xr.attrs["units"] # KeyError
numpy.datetime64をdatetime.datetimeに変換する方法はいろいろあるようです（参考：「numpy.datetime64」からどの様にして年月日を得るか）。最も簡単な方法はpandasを使う方法です。
import pandas as pd
time_list = [pd.to_datetime(time) for time in time_xr.values]
print(time_list)
# [Timestamp('2021-03-06 18:00:00'), Timestamp('2021-03-06 21:00:00'), Timestamp('2021-03-07 00:00:00'), Timestamp('2021-03-07 03:00:00'), Timestamp('2021-03-07 06:00:00'), Timestamp('2021-03-07 09:00:00')]
この方法だと厳密にはpandasのTimestamp型で返ってくるのですが、timedeltaの加算やstrftimeでの文字列への変換はdatetime.datetimeと同じように行えます。どうしてもdatetime.datetimeに直したい場合は、.to_pydatetime()を使います。
time_list = [pd.to_datetime(time).to_pydatetime() for time in time_xr.values]
あるいは、このためだけにPandasをimportしたくない場合は次のようにします。
time_list = [datetime.utcfromtimestamp(time.astype(datetime)*1e-9) for time in time_xr.values]
おわりに
netCDF4ライブラリはnetCDFの中身を直接読み込みます。一方、xarrayでは時刻や欠損値などをよしなに扱ってくれる代わりに、元の属性のうち不要なものは読み出せません(e.g. units, _FillValue)。この点は設計思想の違いを感じました。
このような細かい挙動の違いはあるにせよ、netCDF4をこれまで使っていた人がxarrayへ移行するのはそれほど難しくないように思います。