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

@fujiisoup さんの記事より
多次元データ解析ライブラリ　xarray
こちらの方がだいぶ丁寧ですし、紹介している機能も豊富です。

使い方の確認

インポート

import numpy as np
import pandas as pd
import xarray as xr

ちなみに私はanaconda3を使っています。xarrayはデフォルトでは入っていませんでしたが、conda install xarrayで素直にインストールできました。

xarray型データ(DataArray)の生成

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

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情報はデフォルト値となる。 次に、軸が何を表すかを明示する。

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がそれぞれどのような値であるかを定義し、入力する。

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配列のように扱えるのだろうか。

data[0,0]

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

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

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

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

その他の確認

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

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

np.sum(data, axis=1)

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

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

data.sum(dim='time2')

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

scipy関数も一応試してみる。

from scipy.signal import medfilt2d

medfilt2d(data, kernel_size=3)

array([[0.        , 0.30608039, 0.        ],
       [0.        , 0.30608039, 0.        ]])

scipy関数も普通に適用できるみたい。この場合は、xarray配列が出力されるわけではないのか。

行列演算

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

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に合わせてしっかり作ってあげなくてはいけない。

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型では？

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を転置しても、各軸の持つ意味は変わらないため、ブロードキャストには影響しない。

xr_data1.T * xr_data2

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

内挿

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

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へのロケーションによる指定

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ライブラリについて、私が使うケースを想定して、確認したいことをつらつら書いていきました。
まだまだ使いどころはありそうですので、これからメインで十分使っていけそうです。
​