PythonでnetCDFを読み込む際のnetCDF4とxarrayの違い

# はじめに
Pythonで気象データを扱う際にnetCDFファイルを読み込むことがあります。以前は[netCDF4](https://unidata.github.io/netcdf4-python/)というライブラリを使って読み込むことが多かったのですが、現在では[xarray](http://xarray.pydata.org/en/stable/)が使われることも多いようです。

netCDF4ライブラリ（以下netCDF4と言う場合はファイルフォーマットではなくライブラリを指す）とxarrayとで文法は似ている点が多いのですが、一部取り扱いに違う点があります。そこで気づいた範囲について差異をまとめました。

実行環境については、[Macでmatplotlib+cartopyの環境構築](https://qiita.com/wm-ytakano/items/8445db3456e4c35dbb54)の記事を参照下さい。

# netCDF4とxarrayの比較
## 前準備
netCDFファイルを手元に用意します。例えば、`Z__C_RJTD_20210306180000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin`というGRIB形式のMSMのGPVファイルがあるとして、これをwgrib2を使って適当にnetCDFファイルに変換します。ここでは500hPa面のジオポテンシャル高度だけ切り出していますが、他の変数も一緒に変換しても構いません。

```shell
$ 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)を書いておきます。

```python
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)のアクセス方法が異なります。

```python
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で異なります。

```python
# 値を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に渡すことができます。

```Python
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の海洋最上層の水温(≃海面水温)のデータを読み込んでみます。

```Python
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は以下の方法で相互変換可能です。

```python
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で格納された値そのままになっています。

```Python
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のリストに変換したい場合、少々アドホックですが例えば次のようにします。

```Python
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"]`でアクセスすることはできません。

```Python
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」からどの様にして年月日を得るか](https://funmatu.wordpress.com/2018/03/10/%E3%80%8Cnumpy-datetime64%E3%80%8D%E3%81%8B%E3%82%89%E3%81%A9%E3%81%AE%E6%A7%98%E3%81%AB%E3%81%97%E3%81%A6%E5%B9%B4%E6%9C%88%E6%97%A5%E3%82%92%E5%BE%97%E3%82%8B%E3%81%8B/)）。最も簡単な方法はpandasを使う方法です。

```Python
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()`を使います。

```Python
time_list = [pd.to_datetime(time).to_pydatetime() for time in time_xr.values]
```

あるいは、このためだけにPandasをimportしたくない場合は次のようにします。

```Python
time_list = [datetime.utcfromtimestamp(time.astype(datetime)*1e-9) for time in time_xr.values]
```

# おわりに
netCDF4ライブラリはnetCDFの中身を直接読み込みます。一方、xarrayでは時刻や欠損値などをよしなに扱ってくれる代わりに、元の属性のうち不要なものは読み出せません(e.g. units, _FillValue)。この点は設計思想の違いを感じました。

このような細かい挙動の違いはあるにせよ、netCDF4をこれまで使っていた人がxarrayへ移行するのはそれほど難しくないように思います。