xarrayの使い方　－　備忘録

#xarrayの使い方

xarrayはPython用のライブラリ。多次元配列を扱いやすいように設計されている。
numpyでも多次元配列を効率的に扱えるが、xarrayでは、各軸が何を示しているのかをデータ自体に含めることができる（indexing）という点が異なる。
また、配列に対してindexingできるといえば、pandasだが、 pandasが得意とするとことは2次元配列であり、多次元配列は苦手な模様。 多次元配列に拡張するならばxarrayが便利のよう。
私自身はもっぱら画像を扱っており、普段はNumpy使いのため、このノートブックでは、xarrayの基本をさらっとおさえて、numpyとの互換性やnumpyをベースにしたxarrayの使い心地について確認する。 多次元配列と言いながら、2次元配列しか動かしていないということに最後に気づきましたが・・・。

## 参考ページ
本家
[xarray: N-D labeled arrays and datasets in Python](http://xarray.pydata.org/en/stable/index.html)

@fujiisoup さんの記事より
[多次元データ解析ライブラリ　xarray](https://qiita.com/fujiisoup/items/0e9b2e74902551aa2855)
こちらの方がだいぶ丁寧ですし、紹介している機能も豊富です。
## 使い方の確認
###インポート
~~~python
import numpy as np
import pandas as pd
import xarray as xr
~~~
ちなみに私はanaconda3を使っています。xarrayはデフォルトでは入っていませんでしたが、`conda install xarray`で素直にインストールできました。
### xarray型データ(DataArray)の生成

xarray型（？）のデータはDataArrayで作る。DataSetクラスもあるのですが、本記事ではスルーします。

~~~python
data = xr.DataArray(np.random.rand(2,3))
data
<xarray.DataArray (dim_0: 2, dim_1: 3)>
array([[0.806013, 0.534489, 0.460619],
       [0.291888, 0.422426, 0.904364]])
Dimensions without coordinates: dim_0, dim_1
~~~

ただのndarrayを入力すると、index情報はデフォルト値となる。 次に、軸が何を表すかを明示する。

~~~python
data = xr.DataArray(np.random.rand(2,3), dims=['time1', 'time2'])
data

<xarray.DataArray (time1: 2, time2: 3)>
array([[0.343279, 0.718976, 0.912658],
       [0.228575, 0.602243, 0.450417]])
Dimensions without coordinates: time1, time2
~~~

各軸がtime1, time2を表すことを明示できた。 さらにtime1, time2がそれぞれどのような値であるかを定義し、入力する。

~~~python
time1, time2 = np.arange(2), np.arange(3)
data = xr.DataArray(np.random.rand(2,3), dims=['time1', 'time2'], coords=[('time1', time1), ('time2', time2)])
data

<xarray.DataArray (time1: 2, time2: 3)>
array([[0.801929, 0.36534 , 0.008653],
       [0.037607, 0.598957, 0.092912]])
Coordinates:
  * time1    (time1) int32 0 1
  * time2    (time2) int32 0 1 2
~~~

各軸の意味を明示した多次元配列が作成できた。自分でClassを作ってもいいけど、xarrayであればこのように簡単に作れる。

### 操作
このように作成したxarray配列は,ndarray配列のように扱えるのだろうか。

~~~python
data[0,0]

<xarray.DataArray ()>
array(0.801929)
Coordinates:
    time1    int32 0
    time2    int32 0
~~~

この書き方はxarrayでは次のように、よりデータの意味を明示しながら書くこともできる。

~~~python
data.sel(time1=0, time2=0)

<xarray.DataArray ()>
array(0.801929)
Coordinates:
    time1    int32 0
    time2    int32 0
~~~

その他の確認

~~~python
data[0,:]

<xarray.DataArray (time2: 3)>
array([0.801929, 0.36534 , 0.008653])
Coordinates:
    time1    int32 0
  * time2    (time2) int32 0 1 2

data[0,::-1]

<xarray.DataArray (time2: 3)>
array([0.008653, 0.36534 , 0.801929])
Coordinates:
    time1    int32 0
  * time2    (time2) int32 2 1 0

data * 2.

<xarray.DataArray (time1: 2, time2: 3)>
array([[1.603857, 0.73068 , 0.017307],
       [0.075215, 1.197913, 0.185823]])
Coordinates:
  * time1    (time1) int32 0 1
  * time2    (time2) int32 0 1 2

data + 1.

<xarray.DataArray (time1: 2, time2: 3)>
array([[1.801929, 1.36534 , 1.008653],
       [1.037607, 1.598957, 1.092912]])
Coordinates:
  * time1    (time1) int32 0 1
  * time2    (time2) int32 0 1 2
~~~

スライスや四則演算は普通にできるみたい。

~~~python
np.sum(data, axis=1)

<xarray.DataArray (time1: 2)>
array([1.175922, 0.729476])
Coordinates:
  * time1    (time1) int32 0 1
~~~

numpy関数も適用できるみたい。出力は同じくxarrayのDataArray型となるようだ。
ちなみにこれはxarrayでは次のように直感的に書ける。

~~~python
data.sum(dim='time2')

<xarray.DataArray (time1: 2)>
array([1.175922, 0.729476])
Coordinates:
  * time1    (time1) int32 0 1
~~~
scipy関数も一応試してみる。

~~~python
from scipy.signal import medfilt2d

medfilt2d(data, kernel_size=3)

array([[0.        , 0.30608039, 0.        ],
       [0.        , 0.30608039, 0.        ]])
~~~
scipy関数も普通に適用できるみたい。この場合は、xarray配列が出力されるわけではないのか。

### 行列演算
多次元配列を扱いやすくするということで、行列計算で威力を発揮する。numpyのブロードキャスト周りが、より便利になる。おさらいとしてnumpyでの場合を見てみる。

~~~python
np_data1, np_data2 = np.arange(6).reshape(2,3), np.arange(2)
np_data1, np_data2

(array([[0, 1, 2],
        [3, 4, 5]]), array([0, 1]))

np_data1 * np_data2

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-83-6d392ce997d0> in <module>()
----> 1 np_data1 * np_data2

ValueError: operands could not be broadcast together with shapes (2,3) (2,) 

np_data1.shape, np_data2.shape

((2, 3), (2,))
~~~

次元に整合性がないため、行列演算もされないし、ブロードキャストもされない。 例えばブロードキャストするためには、data2の「形」をdata1に合わせてしっかり作ってあげなくてはいけない。

~~~python
np_data1 * np_data2.reshape(2,1)  # ブロードキャスト

array([[0, 0, 0],
       [3, 4, 5]])

np_data1.T @ np_data2.reshape(2,1) # 行列積

array([[3],
       [4],
       [5]])
~~~

それではxarray型では？

~~~python
xr_data1 = xr.DataArray(np.arange(6).reshape(2,3), dims=['time1', 'time2'])

xr_data2 = xr.DataArray(np.arange(2), dims=['time1'])

xr_data1, xr_data2

(<xarray.DataArray (time1: 2, time2: 3)>
 array([[0, 1, 2],
        [3, 4, 5]])
 Dimensions without coordinates: time1, time2, <xarray.DataArray (time1: 2)>
 array([0, 1])
 Dimensions without coordinates: time1)

xr_data1 * xr_data2

<xarray.DataArray (time1: 2, time2: 3)>
array([[0, 0, 0],
       [3, 4, 5]])
Dimensions without coordinates: time1, time2
~~~
data2の次元はtime1であることを教えている。そのため、data1の1次元目と対応するということがわかるため、ブロードキャストが可能となる。ポイントは、次元の名前とその要素数が同じであること。ブロードキャスト時にいちいち次元（行列のサイズ）を意識しなくてもいいため、便利ではある。
ただし行列の積演算をする`@`演算子（？）については、利用できないようだ。要注意。

data1を転置しても、各軸の持つ意味は変わらないため、ブロードキャストには影響しない。

~~~python
xr_data1.T * xr_data2

<xarray.DataArray (time2: 3, time1: 2)>
array([[0, 3],
       [0, 4],
       [0, 5]])
Dimensions without coordinates: time2, time1
~~~

### 内挿
内挿が簡単にできる。これは便利だと思う。

~~~python
data = xr.DataArray(np.arange(12).reshape(4,3), coords=[('y', np.arange(4)), ('x', np.arange(3)+10)], dims=['y', 'x'])
data

<xarray.DataArray (y: 4, x: 3)>
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11]])
Coordinates:
  * y        (y) int32 0 1 2 3
  * x        (x) int32 10 11 12

data.sel(y=0)

<xarray.DataArray (x: 3)>
array([0, 1, 2])
Coordinates:
    y        int32 0
  * x        (x) int32 10 11 12

data.interp(y=0.5)

<xarray.DataArray (x: 3)>
array([1.5, 2.5, 3.5])
Coordinates:
  * x        (x) int32 10 11 12
    y        float64 0.5

data.interp(y=np.arange(3)+0.5, x=np.arange(2)+10.5)

<xarray.DataArray (y: 3, x: 2)>
array([[2., 3.],
       [5., 6.],
       [8., 9.]])
Coordinates:
  * y        (y) float64 0.5 1.5 2.5
  * x        (x) float64 10.5 11.5
~~~

### アクセスの速度は？
最後に配列の各要素へのアクセス速度が気になったので、簡単なケースで走らせてみました。
比較は3ケースで、

- numpy arrayへのインデックス指定
- xarray dataarrayへのインデックス指定
- xarray datearrayへのロケーションによる指定

~~~python
def access_speed_measurement():
    import datetime
    tmp = 0.
    x_size, y_size = 100, 100
    np_data = np.random.rand(y_size,x_size)
    start_time = datetime.datetime.now()

    for j in range(y_size):
        for i in range(x_size):
            tmp = np_data[j, i]
    end_time = datetime.datetime.now()
    print('Elapsed time for accesing numpy ndarray: {0} [sec]'.format((end_time - start_time).seconds))

    xr_data = xr.DataArray(np_data, coords=[('y', np.arange(y_size)), ('x', np.arange(x_size))], dims=['y', 'x'])
    start_time = datetime.datetime.now()
    for j in range(y_size):
        for i in range(x_size):
            tmp = xr_data[j, i]
    end_time = datetime.datetime.now()
    print('Elapsed time for accesing xarray DateArray by indices: {0} [sec]'.format((end_time - start_time).seconds))

    start_time = datetime.datetime.now()
    for j in range(y_size):
        for i in range(x_size):
            tmp = xr_data.sel(y=j, x=i)
    end_time = datetime.datetime.now()
    print('Elapsed time for accesing xarray DateArray by locations: {0} [sec]'.format((end_time - start_time).seconds))

access_speed_measurement()

Elapsed time for accesing numpy ndarray: 0 [sec]
Elapsed time for accesing xarray DateArray by indices: 2 [sec]
Elapsed time for accesing xarray DateArray by locations: 3 [sec]
~~~
1番目と2番目で大きな差がついてほしくはなかったのですが、そのような結果にならなくて残念。ただし、上記では純粋なアクセス速度の比較にはなっていなくて、tmpへの値の入力時間がの差の可能性もある。純粋なアクセスを比較するためにはどうすればいいんだろう…。
2番目と3番目については、軸の値がインデックスでは何番目に相当するかっていうワンクッションが入ると思うので、そこを考えると予想通りではあります。
## 最後に
ということでxarrayライブラリについて、私が使うケースを想定して、確認したいことをつらつら書いていきました。
まだまだ使いどころはありそうですので、これからメインで十分使っていけそうです。
​