15
16

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

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

Last updated at Posted at 2021-03-07

はじめに

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へ移行するのはそれほど難しくないように思います。

15
16
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
15
16

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?