xarray を用いたデータ解析

# xarray を用いたデータ解析

以前の記事

+ [多次元データ解析ライブラリ　xarray2](https://qiita.com/fujiisoup/items/ffd04f44bfeb876a6623)
+ [多次元データ解析ライブラリ　xarray](https://qiita.com/fujiisoup/items/0e9b2e74902551aa2855)

で紹介したxarrayを用いたデータ解析について、もう少し実用的な面から紹介します。


## 著者について

2017年から、pydata xarray の開発メンバー。
ここ最近はもっぱらユーザーですが、たまにバグ修正のPRを送ったりしています。

本業で行うほぼ全てのデータ解析に xarray を使っています。

## 想定している読者について

numpy についてある程度知っている（np.ndarray を使った操作・演算ができる）ことを想定しています。

# 本記事の構成

この記事では、もう一度 xarray のデータモデルを紹介したあと、基本的な操作であるインデクシングについて述べます。
このインデクシングが xarray の本質だと言っても良いかもしれません。

その後、インデクシングの逆操作である concat を紹介します。

これらを踏まえた上で、実際にデータを作ることを考えた時、どのようなデータ構成にしておけば xarray での解析が容易になるか紹介します。

データの保存・読み込みや、dask を使った out-of-memory 処理についても少し述べます。

なお、本記事はGoogle collaboratoryを使って執筆しています。
[こちら](https://drive.google.com/file/d/1747FtkVJ3KnOszugTGgnmMPIMLPSFg0N/view?usp=sharing)からコピーしてもらうことで、誰でも実行・トレースできます。

# xarray のデータモデル

## 概要

xarray を一言で述べると、 **座標軸付きの多次元配列** です。numpy の nd-array と、pandas の pd.Series を合わせたものだと考えてもよいかもしれません。

使い方に慣れてくると、データ解析の途中で座標のことを考えなくてよくなるので非常に便利です。

一方で、使い方が少し独特であり、学習コストが高いとよく言われます。
この記事は、その学習コストを少しでも低くすることを目指しています。

xarray には大きく2つのクラスがあります。1つは xarray.DataArray で、もう一つは xarray.Dataset です。
それらについて紹介します。

なお、これは [多次元データ解析ライブラリ　xarray](https://qiita.com/fujiisoup/items/ffd04f44bfeb876a6623) で述べた内容とほぼ同一です。

## xarray.DataArray

**座標軸付きの多次元配列** である xarray.DataArray を紹介します。オブジェクトの作り方については後に述べますが、まずは xarray.tutorial にあるデータを用いて説明します。


```python
import numpy as np
import xarray as xr
```


```python
data = xr.tutorial.open_dataset('air_temperature')['air']
data
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lat: 25, lon: 53)&gt;
[3869000 values with dtype=float32]
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



ここで、 `data` は3次元の xr.DataArray インスタンスです。なお、xarray のオフィシャルな略語は `xr` です。

`data` を表示すると、上記のようなサマリが表示されます。一番上の行の `(time: 2920, lat: 25, lon: 53)` は、格納されている配列の次元名およびその要素数を表しています。

このように、xarray では **各次元に名前をつけて扱います**（オフィシャルなドキュメントでは *label* と呼ばれています）。

１軸目が `time` 軸、２軸目が `lat` 軸、３軸目が `lon` 軸です。
このように名前をつけておくことで、軸の順番を気にしなくてよくなります。

さらに、それぞれの軸には座標値が付加されています。座標のリストが `Coordinates` セクションに表示されます。
このデータでは、`time`, `lat`, `lon` の全ての軸に座標がついています。

また、その他のデータも添付できます。それらは `Attributes` セクションにあります。

インデクシングの節で詳しく述べますが、座標軸がついていることで所望の範囲のデータを直感的に取得することができます。例えば 2014年8月30日のデータが欲しい場合、 


```python
data.sel(time='2014-08-30')
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 4, lat: 25, lon: 53)&gt;
array([[[270.29   , 270.29   , ..., 258.9    , 259.5    ],
        [273.     , 273.1    , ..., 262.6    , 264.19998],
        ...,
        [299.5    , 298.79   , ..., 299.     , 298.5    ],
        [299.29   , 299.29   , ..., 299.9    , 300.     ]],

       [[269.9    , 270.1    , ..., 258.6    , 259.5    ],
        [273.     , 273.29   , ..., 262.5    , 265.     ],
        ...,
        [299.19998, 298.9    , ..., 298.5    , 298.     ],
        [299.4    , 299.69998, ..., 299.19998, 299.6    ]],

       [[270.4    , 270.69998, ..., 261.29   , 261.9    ],
        [273.     , 273.6    , ..., 266.4    , 268.6    ],
        ...,
        [297.9    , 297.6    , ..., 298.29   , 298.29   ],
        [298.69998, 298.69998, ..., 299.1    , 299.4    ]],

       [[270.     , 270.4    , ..., 266.     , 266.19998],
        [273.     , 273.9    , ..., 268.1    , 269.69998],
        ...,
        [298.5    , 298.29   , ..., 298.69998, 299.     ],
        [299.1    , 299.19998, ..., 299.5    , 299.69998]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2014-08-30 ... 2014-08-30T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



のようにできます。ここで`.sel`メソッドは座標軸を用いてデータを選ぶために用います。`.sel(time='2014-08-30')` は、`time`軸を参照して`'2014-08-30'` であるデータを選んでくることに相当します。

### DataArray に関する演算

DataArray は np.ndarray のような多次元配列でもあるので、np.ndarrayに対するような演算を行うことができます。


```python
data * 2
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lat: 25, lon: 53)&gt;
array([[[482.4    , 485.     , 487.     , ..., 465.59998, 471.     ,
         477.19998],
        [487.59998, 489.     , 489.4    , ..., 465.59998, 470.59998,
         478.59998],
        [500.     , 499.59998, 497.78   , ..., 466.4    , 472.78   ,
         483.4    ],
        ...,
        [593.2    , 592.39996, 592.8    , ..., 590.8    , 590.2    ,
         589.39996],
        [591.8    , 592.39996, 593.58   , ..., 591.8    , 591.8    ,
         590.39996],
        [592.58   , 593.58   , 594.2    , ..., 593.8    , 593.58   ,
         593.2    ]],

       [[484.19998, 485.4    , 486.19998, ..., 464.     , 467.19998,
         471.59998],
        [487.19998, 488.19998, 488.4    , ..., 462.     , 465.     ,
         471.4    ],
        [506.4    , 505.78   , 504.19998, ..., 461.59998, 466.78   ,
         477.     ],
...
        [587.38   , 587.77997, 590.77997, ..., 590.18   , 589.38   ,
         588.58   ],
        [592.58   , 594.38   , 595.18   , ..., 590.58   , 590.18   ,
         588.77997],
        [595.58   , 596.77997, 596.98   , ..., 591.38   , 590.98   ,
         590.38   ]],

       [[490.18   , 488.58   , 486.58   , ..., 483.37997, 482.97998,
         483.58   ],
        [499.78   , 498.58   , 496.78   , ..., 479.18   , 480.58   ,
         483.37997],
        [525.98   , 524.38   , 522.77997, ..., 479.78   , 485.18   ,
         492.58   ],
        ...,
        [587.58   , 587.38   , 590.18   , ..., 590.58   , 590.18   ,
         589.38   ],
        [592.18   , 593.77997, 594.38   , ..., 591.38   , 591.38   ,
         590.38   ],
        [595.38   , 596.18   , 596.18   , ..., 592.98   , 592.38   ,
         591.38   ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00</pre>



直感的にそうあるべきなように、掛け算はアレイデータに対してのみ行われ、座標軸データは変更されずにそのまま保持されます。

### アレイ以外で保持できるデータ

アレイデータ以外に、Coordinate を保持できることを述べました。Coordinateデータは以下のように分類できます。

- dimension coordinate
- non-dimension coordinate
- scalar coordinate

#### Dimension coordinate

これは、前節で Coordinate として紹介したものです。dimensionの名前(前節のデータでは`time`, `lon`, `lat`)と同じ名前を持っています。

#### Non-dimension coordinate

dimension とは異なる名前の座標も保持できます。どういう場面でそれが必要・便利になるかは後ほど説明しますが、基本的には **保持されているけど演算されないデータ** だと理解すれば良いと思います。

#### Scalar coordinate

座標軸になりうるスカラーを保持できます。これもよく使うデータタイプですが、ひとまずは Non-dimension coordinate のように演算されないけど保持されているデータだと思っておけば良いと思います。

### 基本的なメソッド

np.ndarray にある基本的なメソッドを有しています。

+ `.shape`: 形状を返します
+ `.size`: 合計の大きさを返します
+ `.ndim`: 次元数を返します

そのほか、

+ `.sizes`: 次元名と大きさの連想配列（辞書型）を返します

型に関するもの

+ `.dtype`: 配列の型を返します
+ `.astype()`: 配列の型を変更します
+ `.real`: 全データの実部を返します
+ `.imag`: 全データの虚部を返します

### DataArrayからnp.ndarrayに変換する。

DataArrayの中には、np.ndarray が収納されています。
（補足ですが、実際はnp.ndarray 以外の多次元配列も格納できます。）

`.data` とすることで、中身の配列にアクセスできます。
（なお、`.values` では中身の配列がどんなオブジェクトであれ np.ndarray に変換して返してくれます。）


```python
type(data.data), data.shape
```




    (numpy.ndarray, (2920, 25, 53))



### 座標軸情報を取り出す

座標軸オブジェクトは、辞書のように`[]`内に軸名を渡すことで取り出すことができます。


```python
data['lat']
```




<pre>&lt;xarray.DataArray &#x27;lat&#x27; (lat: 25)&gt;
array([75. , 72.5, 70. , 67.5, 65. , 62.5, 60. , 57.5, 55. , 52.5, 50. , 47.5,
       45. , 42.5, 40. , 37.5, 35. , 32.5, 30. , 27.5, 25. , 22.5, 20. , 17.5,
       15. ], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
Attributes:
    standard_name:  latitude
    long_name:      Latitude
    units:          degrees_north
    axis:           Y</pre>



取り出したものも DataArray オブジェクトになります。


## xarray.Dataset

DataArrayを複数まとめたオブジェクトがDatasetです。


```python
data = xr.tutorial.open_dataset('air_temperature')
data['mean_temperature'] = data['air'].mean('time')
data
```




<pre>&lt;xarray.Dataset&gt;
Dimensions:           (lat: 25, lon: 53, time: 2920)
Coordinates:
  * lat               (lat) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0
  * lon               (lon) float32 200.0 202.5 205.0 ... 325.0 327.5 330.0
  * time              (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Data variables:
    air               (time, lat, lon) float32 241.2 242.5 ... 296.19 295.69
    mean_temperature  (lat, lon) float32 260.37564 260.1826 ... 297.30502
Attributes:
    Conventions:  COARDS
    title:        4x daily NMC reanalysis (1948)
    description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
    platform:     Model
    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...</pre>



この例では、`data` は `air` と `mean_temperature` の2つのDataArrayを保持しています。
それらは `Data variables` というセクションに書かれています。

+ `air` は (`time`, `lon`, `lat`) 次元を持つ３次元配列
+ `mean_temperature` は (`lon`, `lat`) 次元を持つ２次元配列

です。

保持するDataArrayは、座標軸を共有できます。そのため、以下のように複数のDataArrayに対して同時にインデクシングすることが可能です。


```python
data.sel(lat=70, method='nearest')
```




<pre>&lt;xarray.Dataset&gt;
Dimensions:           (lon: 53, time: 2920)
Coordinates:
    lat               float32 70.0
  * lon               (lon) float32 200.0 202.5 205.0 ... 325.0 327.5 330.0
  * time              (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Data variables:
    air               (time, lon) float32 250.0 249.79999 ... 242.59 246.29
    mean_temperature  (lon) float32 264.7681 264.3271 ... 253.58247 257.71475
Attributes:
    Conventions:  COARDS
    title:        4x daily NMC reanalysis (1948)
    description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
    platform:     Model
    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...</pre>



保持されている DataArray は、異なる次元（Dimension）を持っていてもかまいません。軸が共通かどうかは、軸の名前で判断します。そのため、共通でないけど同じ名前の座標軸をもたせるということはできません。

Datasetから DataArray を取得するには、辞書のようにDataArrayの名前をキーとして渡せばよいです。


```python
da = data['air']
da
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lat: 25, lon: 53)&gt;
array([[[241.2    , 242.5    , ..., 235.5    , 238.59999],
        [243.79999, 244.5    , ..., 235.29999, 239.29999],
        ...,
        [295.9    , 296.19998, ..., 295.9    , 295.19998],
        [296.29   , 296.79   , ..., 296.79   , 296.6    ]],

       [[242.09999, 242.7    , ..., 233.59999, 235.79999],
        [243.59999, 244.09999, ..., 232.5    , 235.7    ],
        ...,
        [296.19998, 296.69998, ..., 295.5    , 295.1    ],
        [296.29   , 297.19998, ..., 296.4    , 296.6    ]],

       ...,

       [[245.79   , 244.79   , ..., 243.98999, 244.79   ],
        [249.89   , 249.29   , ..., 242.48999, 244.29   ],
        ...,
        [296.29   , 297.19   , ..., 295.09   , 294.38998],
        [297.79   , 298.38998, ..., 295.49   , 295.19   ]],

       [[245.09   , 244.29   , ..., 241.48999, 241.79   ],
        [249.89   , 249.29   , ..., 240.29   , 241.68999],
        ...,
        [296.09   , 296.88998, ..., 295.69   , 295.19   ],
        [297.69   , 298.09   , ..., 296.19   , 295.69   ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



## DataArray のインスタンス化

ここでは、DataArray オブジェクトを作る方法について説明します。

最も簡単な方法は、`np.ndarray` とそれぞれの軸名を与える方法です。


```python
xr.DataArray(np.arange(12).reshape(3, 4), dims=['x', 'y'])
```




<pre>&lt;xarray.DataArray (x: 3, y: 4)&gt;
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
Dimensions without coordinates: x, y</pre>



上記の例では、3x4 要素を持つ DataArrayを作成しています。1つ目の次元の名前が `x` で、2つ目の次元が `y` です。

座標軸を与えるには、`coords` キーワードに対して辞書型で指定します。


```python
xr.DataArray(np.arange(12).reshape(3, 4), dims=['x', 'y'], 
             coords={'x': [0, 1, 2], 'y': [0.1, 0.2, 0.3, 0.4]})
```




<pre>&lt;xarray.DataArray (x: 3, y: 4)&gt;
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
Coordinates:
  * x        (x) int64 0 1 2
  * y        (y) float64 0.1 0.2 0.3 0.4</pre>



non-dimension coordinate を含めるためには、それが依存する軸名を指定しないといけません。そのために、辞書型の引数（コロンの右側）をタプルにします。1つ目の要素が依存する軸名で、2つ目が配列本体です。


```python
xr.DataArray(np.arange(12).reshape(3, 4), dims=['x', 'y'], 
             coords={'x': [0, 1, 2], 'y': [0.1, 0.2, 0.3, 0.4], 
                     'z': ('x', ['a', 'b', 'c'])})
```




<pre>&lt;xarray.DataArray (x: 3, y: 4)&gt;
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
Coordinates:
  * x        (x) int64 0 1 2
  * y        (y) float64 0.1 0.2 0.3 0.4
    z        (x) &lt;U1 &#x27;a&#x27; &#x27;b&#x27; &#x27;c&#x27;</pre>



ここで、non-dimension coordinate `z` は `x`軸上に定義されたものです。

scalar coordinate も `coords` に渡します。


```python
xr.DataArray(np.arange(12).reshape(3, 4), dims=['x', 'y'], 
             coords={'x': [0, 1, 2], 'y': [0.1, 0.2, 0.3, 0.4], 'scalar': 3})
```




<pre>&lt;xarray.DataArray (x: 3, y: 4)&gt;
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
Coordinates:
  * x        (x) int64 0 1 2
  * y        (y) float64 0.1 0.2 0.3 0.4
    scalar   int64 3</pre>



scalar coordinate は、どの次元にも属していないことがわかります。

軸情報の他に、attribute や DataArray の名前も渡すことができます。ただし、attribute は演算の途中で失われてしまうので、まとめてファイルに保存したいものなどを入れる程度にしたほうがよいかもしれません。

このあたりは後で説明する、保存のためのデータ構造の章で説明したいと思います。

# インデクシング `.isel`, `.sel`, `.interp`, `.reindex`

インデクシングは、最も基本的かつ本質的な操作です。それをシンプルにするために xarray の開発がスタートしたと言っても過言ではないと思います。

xarray では、位置ベースインデクシングと座標軸ベースのインデクシングが可能です。

## 位置ベースインデクシング `.isel`

これは、np.ndarray などの一般的な配列に対するインデクシングと同様のものです。
`data[i,j,k]` のようにカギカッコの中に整数を引数として与えます。
引数は、 **配列の中の要素の位置** を示すものです。


```python
da[0, :4, 3]
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (lat: 4)&gt;
array([244.     , 244.2    , 247.5    , 266.69998], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5
    lon      float32 207.5
    time     datetime64[ns] 2013-01-01
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



上記の書き方だと、各次元の位置（例えば`time`軸が一番最初の軸で、というようなもの）を覚えておく必要があります。

解析ごとにそれらを覚えておくというのは意外と大変です。以下のような書き方をすると、軸の名前だけを覚えておけばよいことになります



```python
da.isel(time=0, lat=slice(0, 4), lon=3)
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (lat: 4)&gt;
array([244.     , 244.2    , 247.5    , 266.69998], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5
    lon      float32 207.5
    time     datetime64[ns] 2013-01-01
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



このように、`.isel()`メソッドの引数にキーワード形式で軸名を入れ、それに対応するインデックスをイコールの後に与えます。

なお、範囲を指定するときは明示的にスライス（`slice(0, 4)`）を使う必要があります。

## 座標ベースインデクシング `.sel`

実際のデータ解析では、興味のある範囲は座標で指定されることが多いと思います。xarray では `.sel` メソッドを用いることで実現できます。


```python
da.sel(time='2013-01-01', lat=slice(75, 60))
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 4, lat: 7, lon: 53)&gt;
array([[[241.2    , 242.5    , ..., 235.5    , 238.59999],
        [243.79999, 244.5    , ..., 235.29999, 239.29999],
        ...,
        [272.1    , 270.9    , ..., 275.4    , 274.19998],
        [273.69998, 273.6    , ..., 274.19998, 275.1    ]],

       [[242.09999, 242.7    , ..., 233.59999, 235.79999],
        [243.59999, 244.09999, ..., 232.5    , 235.7    ],
        ...,
        [269.19998, 268.5    , ..., 275.5    , 274.69998],
        [272.1    , 272.69998, ..., 275.79   , 276.19998]],

       [[242.29999, 242.2    , ..., 236.09999, 238.7    ],
        [244.59999, 244.39   , ..., 232.     , 235.7    ],
        ...,
        [273.     , 273.5    , ..., 275.29   , 274.29   ],
        [275.5    , 275.9    , ..., 277.4    , 277.6    ]],

       [[241.89   , 241.79999, ..., 235.5    , 237.59999],
        [246.29999, 245.29999, ..., 231.5    , 234.5    ],
        ...,
        [273.29   , 272.6    , ..., 277.6    , 276.9    ],
        [274.1    , 274.     , ..., 279.1    , 279.9    ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 62.5 60.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2013-01-01T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



`.isel` の場合と同様に、キーワードに軸名、イコールの後にほしい座標の値を与えます。
`slice` オブジェクトも使うことができます。その場合は、`slice`の引数に実数を入れればよいです。

ただし、`.sel`メソッドにはいくつか注意点があります。

### おおよその座標値を用いたい時

`.sel` メソッドにはオプション `method`, `tolerance` があります。

最もよく使う `method` は `nearest` だと思います。これは、最も近い座標値を用いてインデクシングしてくれます。


```python
da.sel(lat=76, method='nearest')
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lon: 53)&gt;
array([[241.2    , 242.5    , 243.5    , ..., 232.79999, 235.5    , 238.59999],
       [242.09999, 242.7    , 243.09999, ..., 232.     , 233.59999, 235.79999],
       [242.29999, 242.2    , 242.29999, ..., 234.29999, 236.09999, 238.7    ],
       ...,
       [243.48999, 242.98999, 242.09   , ..., 244.18999, 244.48999, 244.89   ],
       [245.79   , 244.79   , 243.48999, ..., 243.29   , 243.98999, 244.79   ],
       [245.09   , 244.29   , 243.29   , ..., 241.68999, 241.48999, 241.79   ]],
      dtype=float32)
Coordinates:
    lat      float32 75.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



### 座標軸ベースのインデクシングが働かない時

座標軸は自由度があるので、うまくインデクシングできない場合があります。例えば以下のようなケースが該当します。
#### 座標軸が単調増加・単調減少でないとき  
`da.sortby('lat')` のように`sortby`メソッドでその軸方向にソートできます。
#### 座標軸に Nan が含まれている時  
`da.isel(lat=~da['lat'].isnull())` のように、`Nan` でない座標を返す`isnull`メソッドを併用することで、その座標値をスキップした配列を作れます
#### 座標軸に重複値があるとき  
`da.isel(lat=np.unique(da['lat'], return_index=True)[1])` とすると、重複座標に対応する部分をスキップした配列を作れます



## 補間 `.interp`

インデクシングと同じようなシンタックスで補間が可能です。


```python
da.interp(lat=74)
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lon: 53)&gt;
array([[242.23999321, 243.29999998, 243.9799988 , ..., 232.79998779,
        235.41999512, 238.87998962],
       [242.69999081, 243.25999451, 243.53999329, ..., 231.60000001,
        233.1599945 , 235.75999145],
       [243.21998903, 243.07599789, 242.97999266, ..., 232.69998783,
        234.45999444, 237.49999702],
       ...,
       [245.72999275, 245.38999009, 244.68999648, ..., 242.74999088,
        243.20999146, 244.00999451],
       [247.42999566, 246.58999336, 245.48999023, ..., 242.4899933 ,
        243.38999027, 244.58999329],
       [247.00999749, 246.28999329, 245.32999587, ..., 240.84999084,
        241.00999144, 241.74999085]])
Coordinates:
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
    lat      int64 74
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



## 条件によるインデクシング `xr.where` 

複雑な条件を満たすデータだけを選びたいときがあります。例えば上記の例だと、

+ `time` がある曜日の時だけ選びたい
+ `air` がある値以上のものだけを選びたい

とかです。

選択結果が一次元のとき、例えばある曜日を選択するときなどは、Bool型の一次元配列を`.isel`メソッドに渡すことで実現できます。




```python
import datetime

# 週末だけを選びます
is_weekend = da['time'].dt.dayofweek.isin([5, 6])
is_weekend
```




<pre>&lt;xarray.DataArray &#x27;dayofweek&#x27; (time: 2920)&gt;
array([False, False, False, ..., False, False, False])
Coordinates:
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00</pre>




```python
# bool型を.isel メソッドに渡します。
da.isel(time=is_weekend)
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 832, lat: 25, lon: 53)&gt;
array([[[248.59999, 246.89   , ..., 245.7    , 246.7    ],
        [256.69998, 254.59999, ..., 248.     , 250.09999],
        ...,
        [296.69998, 296.     , ..., 295.     , 294.9    ],
        [297.     , 296.9    , ..., 296.79   , 296.69998]],

       [[245.09999, 243.59999, ..., 249.89   , 251.39   ],
        [254.     , 251.79999, ..., 247.79999, 250.89   ],
        ...,
        [296.79   , 296.19998, ..., 294.6    , 294.69998],
        [297.19998, 296.69998, ..., 296.29   , 296.6    ]],

       ...,

       [[242.39   , 241.79999, ..., 247.59999, 247.29999],
        [247.18999, 246.09999, ..., 253.29999, 254.29999],
        ...,
        [294.79   , 295.29   , ..., 297.69998, 297.6    ],
        [296.38998, 297.19998, ..., 298.19998, 298.29   ]],

       [[249.18999, 248.     , ..., 244.09999, 242.5    ],
        [254.5    , 253.     , ..., 251.     , 250.59999],
        ...,
        [294.88998, 295.19998, ..., 298.     , 297.29   ],
        [296.1    , 296.88998, ..., 298.5    , 298.29   ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-05 ... 2014-12-28T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



多次元配列のデータをもとに選択を行いたいとき、選んだ後の配列の形状を定義できないので、選ばない要素をNanで置き換えるとよいでしょう。
そういうときには `where` 関数が使えます。


```python
da.where(da > 270)
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lat: 25, lon: 53)&gt;
array([[[      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        ...,
        [296.6    , 296.19998, 296.4    , ..., 295.4    , 295.1    ,
         294.69998],
        [295.9    , 296.19998, 296.79   , ..., 295.9    , 295.9    ,
         295.19998],
        [296.29   , 296.79   , 297.1    , ..., 296.9    , 296.79   ,
         296.6    ]],

       [[      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
...
        [293.69   , 293.88998, 295.38998, ..., 295.09   , 294.69   ,
         294.29   ],
        [296.29   , 297.19   , 297.59   , ..., 295.29   , 295.09   ,
         294.38998],
        [297.79   , 298.38998, 298.49   , ..., 295.69   , 295.49   ,
         295.19   ]],

       [[      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        [      nan,       nan,       nan, ...,       nan,       nan,
               nan],
        ...,
        [293.79   , 293.69   , 295.09   , ..., 295.29   , 295.09   ,
         294.69   ],
        [296.09   , 296.88998, 297.19   , ..., 295.69   , 295.69   ,
         295.19   ],
        [297.69   , 298.09   , 298.09   , ..., 296.49   , 296.19   ,
         295.69   ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



ここでは、 `da > 270` に当てはまらない要素を `np.nan` で置き換えています。 


## 発展的なインデクシング

これまで基本的なインデクシングを紹介しましたが、xarray では発展的なインデクシングも可能です。

この図のように、多次元配列の中から別の配列を切り出すことができます。
オフィシャルドキュメントでは Advanced indexingと呼んでいます。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104970/71014c6d-2701-e0a2-68a4-ef05502d997f.png)


これを実現するためには、引数として渡す配列を DataArrayとして定義します。そのとき、それらの配列の次元を新しく作成する配列の次元にしておきます。

上の図に沿うと、新しい次元は `z` なので、`z` 上のDataArrayを定義します。


```python
lon_new = xr.DataArray([60, 61, 62], dims=['z'], coords={'z': ['a', 'b', 'c']})
lat_new = xr.DataArray([16, 46, 76], dims=['z'], coords={'z': ['a', 'b', 'c']})
```

これらを `.sel`メソッドの引数として渡すと `(lon, lat)` がそれぞれ `[(60, 16), (61, 46), (62, 76)]` のもの（に最も近いもの）を選んできてくれます。


```python
da.sel(lon=lon_new, lat=lat_new, method='nearest')
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, z: 3)&gt;
array([[296.29   , 280.     , 241.2    ],
       [296.29   , 279.19998, 242.09999],
       [296.4    , 278.6    , 242.29999],
       ...,
       [298.19   , 279.99   , 243.48999],
       [297.79   , 279.69   , 245.79   ],
       [297.69   , 279.79   , 245.09   ]], dtype=float32)
Coordinates:
    lat      (z) float32 15.0 45.0 75.0
    lon      (z) float32 200.0 200.0 200.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
  * z        (z) &lt;U1 &#x27;a&#x27; &#x27;b&#x27; &#x27;c&#x27;
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



新たに得られた配列は、もはや`lon`, `lat`軸には依存していなくて、 `z` に依存していることがわかります。

また、同様の方法を補間にも用いることができます。


```python
da.interp(lon=lon_new, lat=lat_new)
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, z: 3)&gt;
array([[nan, nan, nan],
       [nan, nan, nan],
       [nan, nan, nan],
       ...,
       [nan, nan, nan],
       [nan, nan, nan],
       [nan, nan, nan]])
Coordinates:
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
    lon      (z) int64 60 61 62
    lat      (z) int64 16 46 76
  * z        (z) &lt;U1 &#x27;a&#x27; &#x27;b&#x27; &#x27;c&#x27;
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



# データの結合 `xr.concat`

複数の xarray オブジェクトを結合させたいこともよくあります。
例えば、色々なパラメータで得られた実験データ・シミュレーションデータを合わせて解析したいときとかでしょうか。


これは `np.concatenate`, `np.stack` に相当する関数です。
インデクシングの逆操作だと思うとわかりやすいと思います。

シンタックスは

```python
xr.concat([data1, data2, ..., ], dim='dim')
```

です。
最初の引数として、複数の DataArray （もしくは Dataset）を渡します。
2つ目の引数として、結合の方向を指定します。
結合方向によって、 concatenate 動作と stack 動作に分けられます。

### 既存の軸を指定する(concatenate動作)

`xr.concat`関数の基本的な操作は、既存の軸方向に複数の DataArray をつなげることができます。

例えば以下の2つの DataArray 




```python
da0 = da.isel(time=slice(0, 10))
da1 = da.isel(time=slice(10, None))
da0
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 10, lat: 25, lon: 53)&gt;
array([[[241.2    , 242.5    , ..., 235.5    , 238.59999],
        [243.79999, 244.5    , ..., 235.29999, 239.29999],
        ...,
        [295.9    , 296.19998, ..., 295.9    , 295.19998],
        [296.29   , 296.79   , ..., 296.79   , 296.6    ]],

       [[242.09999, 242.7    , ..., 233.59999, 235.79999],
        [243.59999, 244.09999, ..., 232.5    , 235.7    ],
        ...,
        [296.19998, 296.69998, ..., 295.5    , 295.1    ],
        [296.29   , 297.19998, ..., 296.4    , 296.6    ]],

       ...,

       [[244.79999, 244.39   , ..., 242.7    , 244.79999],
        [246.7    , 247.09999, ..., 237.79999, 240.2    ],
        ...,
        [297.79   , 297.19998, ..., 296.4    , 295.29   ],
        [297.9    , 297.69998, ..., 297.19998, 297.     ]],

       [[243.89   , 243.79999, ..., 240.29999, 242.59999],
        [245.5    , 245.79999, ..., 236.59999, 239.     ],
        ...,
        [297.6    , 297.     , ..., 295.6    , 295.     ],
        [298.1    , 297.69998, ..., 296.79   , 297.1    ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2013-01-03T06:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



を `time`方向につなげることで、もとのデータを復元できます。


```python
xr.concat([da0, da1], dim='time')
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lat: 25, lon: 53)&gt;
array([[[241.2    , 242.5    , 243.5    , ..., 232.79999, 235.5    ,
         238.59999],
        [243.79999, 244.5    , 244.7    , ..., 232.79999, 235.29999,
         239.29999],
        [250.     , 249.79999, 248.89   , ..., 233.2    , 236.39   ,
         241.7    ],
        ...,
        [296.6    , 296.19998, 296.4    , ..., 295.4    , 295.1    ,
         294.69998],
        [295.9    , 296.19998, 296.79   , ..., 295.9    , 295.9    ,
         295.19998],
        [296.29   , 296.79   , 297.1    , ..., 296.9    , 296.79   ,
         296.6    ]],

       [[242.09999, 242.7    , 243.09999, ..., 232.     , 233.59999,
         235.79999],
        [243.59999, 244.09999, 244.2    , ..., 231.     , 232.5    ,
         235.7    ],
        [253.2    , 252.89   , 252.09999, ..., 230.79999, 233.39   ,
         238.5    ],
...
        [293.69   , 293.88998, 295.38998, ..., 295.09   , 294.69   ,
         294.29   ],
        [296.29   , 297.19   , 297.59   , ..., 295.29   , 295.09   ,
         294.38998],
        [297.79   , 298.38998, 298.49   , ..., 295.69   , 295.49   ,
         295.19   ]],

       [[245.09   , 244.29   , 243.29   , ..., 241.68999, 241.48999,
         241.79   ],
        [249.89   , 249.29   , 248.39   , ..., 239.59   , 240.29   ,
         241.68999],
        [262.99   , 262.19   , 261.38998, ..., 239.89   , 242.59   ,
         246.29   ],
        ...,
        [293.79   , 293.69   , 295.09   , ..., 295.29   , 295.09   ,
         294.69   ],
        [296.09   , 296.88998, 297.19   , ..., 295.69   , 295.69   ,
         295.19   ],
        [297.69   , 298.09   , 298.09   , ..., 296.49   , 296.19   ,
         295.69   ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



ここで、 `time` 軸も同じように復元されていることにも気をつけてください。

### 新しい軸を指定する （stack 動作）
例えば以下のように2つのDataArrayがあるとします。



```python
da0 = da.isel(time=0)
da1 = da.isel(time=1)
da0
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (lat: 25, lon: 53)&gt;
array([[241.2    , 242.5    , 243.5    , ..., 232.79999, 235.5    , 238.59999],
       [243.79999, 244.5    , 244.7    , ..., 232.79999, 235.29999, 239.29999],
       [250.     , 249.79999, 248.89   , ..., 233.2    , 236.39   , 241.7    ],
       ...,
       [296.6    , 296.19998, 296.4    , ..., 295.4    , 295.1    , 294.69998],
       [295.9    , 296.19998, 296.79   , ..., 295.9    , 295.9    , 295.19998],
       [296.29   , 296.79   , 297.1    , ..., 296.9    , 296.79   , 296.6    ]],
      dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
    time     datetime64[ns] 2013-01-01
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



これを、`time`軸方向に結合するには以下のようにします。


```python
xr.concat([da0, da1], dim='time')
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2, lat: 25, lon: 53)&gt;
array([[[241.2    , 242.5    , 243.5    , ..., 232.79999, 235.5    ,
         238.59999],
        [243.79999, 244.5    , 244.7    , ..., 232.79999, 235.29999,
         239.29999],
        [250.     , 249.79999, 248.89   , ..., 233.2    , 236.39   ,
         241.7    ],
        ...,
        [296.6    , 296.19998, 296.4    , ..., 295.4    , 295.1    ,
         294.69998],
        [295.9    , 296.19998, 296.79   , ..., 295.9    , 295.9    ,
         295.19998],
        [296.29   , 296.79   , 297.1    , ..., 296.9    , 296.79   ,
         296.6    ]],

       [[242.09999, 242.7    , 243.09999, ..., 232.     , 233.59999,
         235.79999],
        [243.59999, 244.09999, 244.2    , ..., 231.     , 232.5    ,
         235.7    ],
        [253.2    , 252.89   , 252.09999, ..., 230.79999, 233.39   ,
         238.5    ],
        ...,
        [296.4    , 295.9    , 296.19998, ..., 295.4    , 295.1    ,
         294.79   ],
        [296.19998, 296.69998, 296.79   , ..., 295.6    , 295.5    ,
         295.1    ],
        [296.29   , 297.19998, 297.4    , ..., 296.4    , 296.4    ,
         296.6    ]]], dtype=float32)
Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 2013-01-01T06:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



このように、 `da0`, `da1` に `time`という scalar coordinate が存在している場合、その方向につなげることで
（`dim='time'` というように`dim`キーワードに名前を指定することで）、`time` 軸を新たな座標軸（dimension coordinate）にすることができます。



ほかにも、`dim`キーワードにDataArrayを与えることで、その名前の座標軸を新しく作ることも可能です。

# ファイルIO

xarray のもう一つのメリットは、DataArray や Dataset などのデータを無矛盾・自己完結的に保存できることです。
いくつかファイルフォーマットをサポートしていますが、最も使いやすいのは [netCDF](https://www.unidata.ucar.edu/software/netcdf/docs/netcdf_introduction.html) でしょう。

Wikipedia にも説明があります。

> NetCDFは、コンピュータの機種に依存しないバイナリ形式であり（機種非依存）、データを配列として読み書きすることができ（配列指向）、さらにデータに加えそのデータに関するに関する説明を格納できる（自己記述的）という特徴がある。NetCDFは、オープンな地理情報標準を策定する国際的なコンソーシアムである、Open Geospatial Consortiumにおける国際標準である。

現在、バージョン3と4がメンテナンスされていますが、バージョン4を使っておくほうが無難でしょう。
バージョン4はHDF5の標準に準拠していますので、通常の HDF5 リーダでも内容を読み込むことができます。

なお、正式な拡張子は `.nc` です。

## netCDF 形式でファイルを保存・読み込みをする

netCDF 形式で保存するためには、 パッケージ `netcdf4` が必要です。

```bash
pip install netcdf4
```
もしくは
```bash
conda install netcdf4
```
を実行してシステムにインストールしておきましょう。

xarray オブジェクトを netCDF 形式で保存するには、`.to_netcdf`メソッドを実行します。


```python
data.to_netcdf('test.nc')
```

    /home/keisukefujii/miniconda3/envs/xarray/lib/python3.7/site-packages/ipykernel_launcher.py:1: SerializationWarning: saving variable air with floating point data as an integer dtype without any _FillValue to use for NaNs
      """Entry point for launching an IPython kernel.


これで、現在のパスに `test.nc` が保存されたと思います。

逆に、保存されたファイルを読み込むためには、 `xr.load_dataset` を行います。


```python
xr.load_dataset('test.nc')
```




<pre>&lt;xarray.Dataset&gt;
Dimensions:           (lat: 25, lon: 53, time: 2920)
Coordinates:
  * lat               (lat) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0
  * lon               (lon) float32 200.0 202.5 205.0 ... 325.0 327.5 330.0
  * time              (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Data variables:
    air               (time, lat, lon) float32 241.2 242.5 ... 296.19 295.69
    mean_temperature  (lat, lon) float32 260.37564 260.1826 ... 297.30502
Attributes:
    Conventions:  COARDS
    title:        4x daily NMC reanalysis (1948)
    description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
    platform:     Model
    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...</pre>



このように、配列の部分だけでなく、軸の名前、座標の値、 `attributes` にあるその他の情報まで保存されていることがわかります。

### `scipy.io.netcdf` に関する注意

xarray の `.to_netcdf` は、デフォルトでは `netcdf4` ライブラリを使いますが、システムにインストールされていない場合、 `scipy.io.netcdf` を使います。 ただし `scipy.io.netcdf` でサポートされているのはバージョン3で、バージョン4では互換性がないので混乱の元になります。

システムに `netcdf4` をインストールするのを忘れないようにしましょう。
もしくは、 `.to_netcdf('test.nc', engine='netcdf4')` というように、陽に使うパッケージを指定することもできます。

## Out-of-memory 配列

netCDF 形式のファイルはランダム読み込みができます。そのため、ファイルのある部分のデータだけを使いたいというようなときは、全てのデータをまずハードディスクから読み込むのではなく、必要なときに読み込むようにすれば効率がよいでしょう。
特に、メモリに格納できないくらい大きなサイズのファイルを扱うときは、そういった処理が本質的に重要になってきます。
これを out-of-memory 処理と呼びます。

xarray における out-of-memory 処理は結構重要な機能と位置づけられていて、様々な実装が進んでいます。
その詳細はまた別の記事で紹介することにして、ここでは最も基本的な内容だけを述べようと思います。

読み込み時に `load_dataset` ではなく `open_dataset` をすることで、out-of-memory 処理が行われます。



```python
unloaded_data = xr.open_dataset('test.nc')
unloaded_data
```




<pre>&lt;xarray.Dataset&gt;
Dimensions:           (lat: 25, lon: 53, time: 2920)
Coordinates:
  * lat               (lat) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0
  * lon               (lon) float32 200.0 202.5 205.0 ... 325.0 327.5 330.0
  * time              (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Data variables:
    air               (time, lat, lon) float32 ...
    mean_temperature  (lat, lon) float32 ...
Attributes:
    Conventions:  COARDS
    title:        4x daily NMC reanalysis (1948)
    description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
    platform:     Model
    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...</pre>



`Data variables` のところを見てもらうと、例えば `air` のところに数字が入っておらず、単に `...` となっていることがわかります。

これはデータが読み込まれていないことを示しています。
ただし座標軸は常に読み込まれるので、上記のインデクシング作業はこのまま可能です。


```python
unloaded_data['air'].sel(lat=60, method='nearest')
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lon: 53)&gt;
[154760 values with dtype=float32]
Coordinates:
    lat      float32 60.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Attributes:
    long_name:     4xDaily Air temperature at sigma level 995
    units:         degK
    precision:     2
    GRIB_id:       11
    GRIB_name:     TMP
    var_desc:      Air temperature
    dataset:       NMC Reanalysis
    level_desc:    Surface
    statistic:     Individual Obs
    parent_stat:   Other
    actual_range:  [185.16 322.1 ]</pre>



変数を選んだり、インデクシングしたりしてもまだロードされません。
演算されたり、`.data` や `.values` が呼ばれるとメモリに初めてロードされることになります。

明示的にロードするには、`.compute()` メソッドを使っても良いと思います。


```python
unloaded_data.compute()
```




<pre>&lt;xarray.Dataset&gt;
Dimensions:           (lat: 25, lon: 53, time: 2920)
Coordinates:
  * lat               (lat) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0
  * lon               (lon) float32 200.0 202.5 205.0 ... 325.0 327.5 330.0
  * time              (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
Data variables:
    air               (time, lat, lon) float32 241.2 242.5 ... 296.19 295.69
    mean_temperature  (lat, lon) float32 260.37564 260.1826 ... 297.30502
Attributes:
    Conventions:  COARDS
    title:        4x daily NMC reanalysis (1948)
    description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
    platform:     Model
    references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...</pre>



### 注意

`open_dataset`, `open_dataarray` では、ファイルをロックし今後の読み込みに備えます。そのため、別のプログラムで同じ名前のファイルを触るとエラーが起こります。

例えば、時間のかかる計算の結果を `.nc` ファイルとしてその都度出力しているときを考えます。途中経過を見たいと思って `open_dataset` で開くと、その後そのファイルに書き込もうとすると失敗します。
失敗して止まるだけならましですが、既存のデータを壊してしまうことも多いようです。
この out-of-memory 処理は注意して使う必要があります。

以下のように、 `with` 文を使うと、確実にファイルを閉じてロックを解除してくれます。


```python
with xr.open_dataset('test.nc') as f:
  print(f)
```

    <xarray.Dataset>
    Dimensions:           (lat: 25, lon: 53, time: 2920)
    Coordinates:
      * lat               (lat) float32 75.0 72.5 70.0 67.5 ... 22.5 20.0 17.5 15.0
      * lon               (lon) float32 200.0 202.5 205.0 ... 325.0 327.5 330.0
      * time              (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00
    Data variables:
        air               (time, lat, lon) float32 ...
        mean_temperature  (lat, lon) float32 ...
    Attributes:
        Conventions:  COARDS
        title:        4x daily NMC reanalysis (1948)
        description:  Data is from NMC initialized reanalysis\n(4x/day).  These a...
        platform:     Model
        references:   http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...


# おすすめデータ構造

上記のように、netcdf 形式でデータを保存しておくと xarray でスムースにデータ解析を行うことができます。
実験やシミュレーションなどでデータを作成する時には、その後の解析のことを考えたデータ構造を作っておいて、それを含めた netcdf 形式で保存しておくとよいでしょう。

以下を抑えたデータにしておくことをお勧めします。

+ 座標軸データも保存する
+ 実験時刻やスキャンするパラメータは、 `attributes` ではなく `scalar coordinate` とする
+ `attributes` には変化させないもの（実験に使った機器の情報や、データに関する説明文など） を入れるとよいでしょう。

実験の場合、以下のような構造で保存すれば良いと思います。




```python
from datetime import datetime
```


```python
# これをカメラで撮影した画像データだとします。
raw_data = np.arange(512 * 2048).reshape(512, 2048)
# 画像のx軸・y軸の座標です。
x = np.arange(512)
y = np.arange(2048)

# 画像の各ピクセルに対する位置データ
r = 300.0
X = r * np.cos(x / 512)[:, np.newaxis] + r * np.sin(y / 2048)
Y = r * np.sin(x / 512)[:, np.newaxis] - r * np.cos(y / 2048)

# 実験で使ったパラメータ
p0 = 3.0
p1 = 'normal'

# 実験時刻
now = datetime.now()

data = xr.DataArray(
    raw_data, dims=['x', 'y'],
    coords={'x': x, 'y': y, 
            'X': (('x', 'y'), X),
            'Y': (('x', 'y'), Y),
            'p0': p0, 'p1': p1,
            'datetime': now},
    attrs={'camera type': 'flash 4', 
           'about': 'Example of the recommended data structure.'})
data
```




<pre>&lt;xarray.DataArray (x: 512, y: 2048)&gt;
array([[      0,       1,       2, ...,    2045,    2046,    2047],
       [   2048,    2049,    2050, ...,    4093,    4094,    4095],
       [   4096,    4097,    4098, ...,    6141,    6142,    6143],
       ...,
       [1042432, 1042433, 1042434, ..., 1044477, 1044478, 1044479],
       [1044480, 1044481, 1044482, ..., 1046525, 1046526, 1046527],
       [1046528, 1046529, 1046530, ..., 1048573, 1048574, 1048575]])
Coordinates:
  * x         (x) int64 0 1 2 3 4 5 6 7 8 ... 504 505 506 507 508 509 510 511
  * y         (y) int64 0 1 2 3 4 5 6 7 ... 2041 2042 2043 2044 2045 2046 2047
    X         (x, y) float64 300.0 300.1 300.3 300.4 ... 414.7 414.8 414.9 414.9
    Y         (x, y) float64 -300.0 -300.0 -300.0 -300.0 ... 89.66 89.79 89.91
    p0        float64 3.0
    p1        &lt;U6 &#x27;normal&#x27;
    datetime  datetime64[ns] 2020-11-13T03:31:56.219372
Attributes:
    camera type:  flash 4
    about:        Example of the recommended data structure.</pre>



実験中に変わりうるものを scalar coordinate にまとめておくことで、後にそれらに対する依存性を知りたい時にも解析が容易になります。

# その他の有用な関数・メソッド

## Reduction `.sum`, `.mean`, etc

`np.sum`, `np.mean` などの reduction 処理はメソッドとして提供されています。
例えば `np.sum` では、どの軸方向に足し合わせるかを指定できます。xarray では同様に、軸を軸名として指定できます。


```python
data = xr.tutorial.open_dataset('air_temperature')['air']

data.sum('lat')
```




<pre>&lt;xarray.DataArray &#x27;air&#x27; (time: 2920, lon: 53)&gt;
array([[6984.9497, 6991.6606, 6991.5303, ..., 6998.77  , 7007.8804,
        7016.5605],
       [6976.4307, 6988.45  , 6993.2407, ..., 6994.3906, 7006.7505,
        7019.941 ],
       [6975.2603, 6982.02  , 6988.77  , ..., 6992.0503, 7004.9404,
        7020.3506],
       ...,
       [6990.7505, 6998.3496, 7013.3496, ..., 6995.05  , 7008.6504,
        7019.4497],
       [6984.95  , 6991.6504, 7007.949 , ..., 6994.15  , 7008.55  ,
        7020.8506],
       [6981.75  , 6983.85  , 6997.0503, ..., 6985.6494, 6999.2495,
        7012.0493]], dtype=float32)
Coordinates:
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00</pre>



このようにある軸方向に足し合わせを行えば、その方向の座標軸がなくなります。
他の座標軸は残っています。

なお、`.sum` などを実行する場合、デフォルトでは `np.sum` ではなく `np.nansum` が呼ばれ、`Nan`をスキップした計算が行われます。 
`Nan` をスキップしてほしくない場合(`np.nansum`ではなく`np.sum`を使いたい場合）は、`skipna=False` を与える必要があります。

## 描画 `.plot`

`DataArray` には `plot` メソッドが備わっており、データの内容を簡単に可視化することができます。



```python
data.sel(time='2014-08-30 08:00', method='nearest').plot()
```




    <matplotlib.collections.QuadMesh at 0x7f9a98557310>




    
![output_102_1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/104970/e25e2ded-9062-1b5a-0a11-db8d0cd8e4d7.png)

    


論文・報告の最終的な図を作るためというよりは、簡易的なものと思っておいたほうがよいかもしれませんが、Jupyter や Ipython などでデータを可視化して理解しながら解析を進めるために使うには非常に便利です。