2015-09-24

Python Dask で並列 DataFrame 処理

Python pandas Dask

はじめに

先日のエントリで少し記載した Dask について、その使い方を書く。Dask を使うと、NumPy や pandas の API を利用して並列計算/分散処理を行うことができる。また、Dask は Out-Of-Core (データ量が多くメモリに乗らない場合) の処理も考慮した実装になっている。

sinhrks.hatenablog.com

上にも書いたが、Daskは NumPy や pandas を置き換えるものではない。数値計算のためのバックエンドとして NumPy や pandas を利用するため、むしろこれらのパッケージが必須である。

Dask は NumPy や pandas の API を完全にはサポートしていないため、並列 / Out-Of-Core 処理が必要な場面では Dask を、他では NumPy / pandas を使うのがよいと思う。pandasとDask のデータはそれぞれ簡単に相互変換できる。

補足とはいえ都度の変換は手間なので、pandas の処理実行時に Dask を利用するオプションをつけるという検討はされている。

インストール

pip もしくは conda で。

pip install dask

準備

まずは必要なパッケージを import する。

import numpy as np
import pandas as pd
import dask.dataframe as dd

np.__version__
# '1.9.3'

pd.__version__
# '0.16.2'

# バージョン表示のためにインポート。
import dask
dask.__version__
# '0.7.1'

`pandas` から `Dask` への変換

サンプルデータはすでにメモリ上にある pd.DataFrame とする。

df = pd.DataFrame({'X': np.arange(10), 
                   'Y': np.arange(10, 20),
                   'Z': np.arange(20, 30)},
                  index=list('abcdefghij'))
df
#    X   Y   Z
# a  0  10  20
# b  1  11  21
# c  2  12  22
# d  3  13  23
# e  4  14  24
# f  5  15  25
# g  6  16  26
# h  7  17  27
# i  8  18  28
# j  9  19  29

pandas のデータ構造から Dask に変換するには dd.from_pandas。2つめの引数でデータをいくつのパーティションに分割するかを指定している。結果は dask.dataframe.DataFrame ( dd.DataFrame ) となる。

divisions はデータがどこで分割されたかを示す。表示から、1つ目のパーティションには .index が "a" 〜 "e" までのデータが、2つ目のには "f" 〜 "j" までのデータが含まれていることがわかる。

重要 dd.DataFrame の処理全般について、行の順序は保証されない。各パーティションには divisions で示される .index を持つ行が含まれるが、パーティション内が常にソートされているとは限らない。

ddf = dd.from_pandas(df, 2)
ddf
# dd.DataFrame<from_pandas-b08addf72f0693a9fa1bb6c21d91f3d4, divisions=('a', 'f', 'j')>

# DataFrame の列名
ddf.columns
# ('X', 'Y', 'Z')

# DataFrame の index
ddf.index
# dd.Index<from_pandas-b08addf72f0693a9fa1bb6c21d91f3d4-index, divisions=('a', 'f', 'j')>

# DataFrame の divisions (パーティションの分割箇所)
ddf.divisions
# ('a', 'f', 'j')
 
# DataFrame のパーティション数
ddf.npartitions
# 2

dd.DataFrame からデータを取得する (計算処理を実行する) には .compute()。結果、元の pd.DataFrame が得られる。

ddf.compute()
#    X   Y   Z
# a  0  10  20
# b  1  11  21
# c  2  12  22
# d  3  13  23
# e  4  14  24
# f  5  15  25
# g  6  16  26
# h  7  17  27
# i  8  18  28
# j  9  19  29

内部処理

ここから、dd.DataFrame でどういった処理ができるのか、内部動作とあわせて記載する。といっても難しいことは全くやっていない。

まず、データ全体について 1 加算する処理を考える。これは各パーティションごとに 1 加算して連結するのと同じ。

ddf + 1
# dd.DataFrame<elemwise-5b9ae0407254158903343113fac41421, divisions=('a', 'f', 'j')>

(ddf + 1).compute()
# 略

f:id:sinhrks:20150924210655p:plain

次に、列ごとの合計をとる処理。これは、各パーティションごとに列の合計をとって連結し、もう一度合計をとる処理と同じ。

列の合計をとるため、結果は Series になる。ddf.sum() の時点では .compute() が呼ばれていないため実際の計算は行われていないが、Dask が結果の型や divisions を推測し正しく表示している。

ddf.sum()
# dd.Series<dataframe-sum--second-7ba12c9d58c17f61406b84b6c30d7a26, divisions=(None, None)>

ddf.sum().compute()
# X     45
# Y    145
# Z    245
# dtype: int64

f:id:sinhrks:20150924210729p:plain

Dask ではこのような形で、計算処理をパーティションごとに並列 / Out-Of-Core 実行できる形に読み替えている。これらの処理は内部的には Computational Graph ( Dask Graph ) として表現され、.compute() によって実行される。

各処理の Dask Graph は、.visualize() メソッドを利用して確認できる。Graph 上で縦につながっていない処理同士は並列で実行できる。

ddf.sum().visualize()

f:id:sinhrks:20150924212007p:plain

各列の平均をとる場合、内部的には各列の .sum()と各列の .count() をそれぞれ計算して除算。

ddf.mean().compute()
# X     4.5
# Y    14.5
# Z    24.5
# dtype: float64

ddf.mean().visualize()

f:id:sinhrks:20150924212112p:plain

DataFrame 同士の演算や、演算をチェインすることもできる。互いのパーティションが異なる場合はそれらが一致するよう調整が行われる。

((ddf - (ddf * 2)) == - ddf).visualize()

f:id:sinhrks:20150924212136p:plain

また、累積関数 ( cumxxx ) やウィンドウ関数 ( rolling_xxx ) なども利用できる。

ddf.cumsum().compute()
#     X    Y    Z
# a   0   10   20
# b   1   21   41
# c   3   33   63
# d   6   46   86
# e  10   60  110
# f  15   75  135
# g  21   91  161
# h  28  108  188
# i  36  126  216
# j  45  145  245

ddf.cumsum().visualize()

f:id:sinhrks:20150924213152p:plain

concat, join などの連結 / 結合もできる。通常の演算と同じく、dd.DataFrame 同士のパーティションは適当に調整される。

df2 = pd.DataFrame({'A': np.arange(5), 
                    'B': np.arange(10, 15)},
                   index=list('adefg'))
df2
#    A   B
# a  0  10
# d  1  11
# e  2  12
# f  3  13
# g  4  14

ddf2 = dd.from_pandas(df2, 2)
ddf2
# dd.DataFrame<from_pandas-667963fc37e22688843f02da80df5963, divisions=('a', 'f', 'g')>

ddf.join(ddf2).compute()
#    X   Y   Z   A   B
# a  0  10  20   0  10
# b  1  11  21 NaN NaN
# c  2  12  22 NaN NaN
# d  3  13  23   1  11
# e  4  14  24   2  12
# f  5  15  25   3  13
# g  6  16  26   4  14
# h  7  17  27 NaN NaN
# i  8  18  28 NaN NaN
# j  9  19  29 NaN NaN

ddf.join(ddf2).visualize()

f:id:sinhrks:20150924214828p:plain

サポートされている処理の一覧は以下のAPIドキュメントを。一部利用できない引数が明記されていないが、次バージョンにて改訂。

API — dask 0.7.1 documentation

9/26 追記 処理結果については、行の順序以外は pandas の処理と一致するはず。例外は quantile のような percentile をとる処理。これらは Out-Of-Core 処理のための近似アルゴリズムを使っており、正確な値とずれることがある。

実データでの利用例

こちらが良エントリ (英語)。

まとめ

Dask を利用して DataFrame を並列処理する方法を記載した。手順は、

dd.from_pandas を利用して pd.DataFrame を dd.DataFrame へ変換。
実行したいメソッド / 演算を dd.DataFrame に対して適用。
.compute() で計算を実行し、結果を取得する。計算処理は Dask にて自動的に並列化される。

最後、pandas 0.16.2 時点では並列処理による速度向上は大きくはない。これは Python の GIL (Global Interpreter Lock ) により並列実行できる処理が限定されているため。今月中にリリース予定の pandas 0.17.0 ではいくつかの処理で Cython から明示的に GIL 解放するよう実装を変更しており、並列化による速度向上は大きくなる。

参考 Pandas Releasing the GIL

2015-09-20

Python でパイプ演算子を使いたい

pandas Python

ネタ記事です。/ This is a joke post which makes no practical sense.

はじめに

Python pandas では主要な操作を以下のようにメソッドチェインの形で書くことができる。

# Python (pandas)
df.assign(x=df['y'] + df['z']).groupby('x').sum()

pandas v0.16.2 で DataFrame と Series に .pipe というメソッドが追加され、このチェインを外部の関数/メソッドに対して連結できるようになった。利用例は以下のリンクを。

補足 matplotlib でも v1.5.0 でラベルデータ対応が追加され、各関数が .pipe から利用できるようになる予定。

このメソッドチェインによる "処理の連結" を、R ( {magrittr} ) のようにパイプ演算子 %>% を使って統一的に書きたい。

# R (magrittr + dplyr)
df %>% mutate(x = y + z) %>% group_by(x) %>% summarize_each(funs(sum))

R と異なり、 Python では自作の演算子を定義することはできない。Python の演算子の中では右方向ビットシフト演算子 >> がパイプ演算子っぽく見えなくもない。これを使って R のような動作をさせたい。

演算子のオーバーライド

まずはサンプルデータを作成する。

import numpy as np
import pandas as pd
pd.__version__
# '0.16.2'

df = pd.DataFrame({'X': [1, 2, 3, 4, 5],
                   'Y': ['a', 'b', 'a', 'b', 'a']})
df
#    X  Y
# 0  1  a
# 1  2  b
# 2  3  a
# 3  4  b
# 4  5  a

普通に .pipe を使うとこのような感じになる。

func = lambda x: x
df.pipe(func)
#    X  Y
# 0  1  a
# 1  2  b
# 2  3  a
# 3  4  b
# 4  5  a

この処理を >> 演算子を使って書きたい。そのままでは、演算方法が未定義のため TypeError となる。

df >> func
# TypeError: unsupported operand type(s) for >>: 'DataFrame' and 'function'

Python では各演算子に対応する処理はマジックメソッドとして実装する。そのため、>> に対応する __rshift__ に .pipe と同じ処理をさせればよい。また、R のように右辺には関数を渡したいので、pandas のメソッドを適当にラップする関数 pipify を定義する。

pd.DataFrame.__rshift__ = pd.DataFrame.pipe

def pipify(f):
    def wrapped(*args, **kwargs):
        return lambda self: f(self, *args, **kwargs)
    return wrapped

head = pipify(pd.DataFrame.head)

できた。

df >> head(n=2)
#    X  Y
# 0  1  a
# 1  2  b

groupby も同じように書ける。集約関数は DataFrame や DataFrameGroupBy 双方から使うため以下のように定義する。

groupby = pipify(pd.DataFrame.groupby)

def sum(*args, **kwargs):
    return lambda self: self.sum(*args, **kwargs)

また、groupby した結果からも >> で演算を連結したい。そのためにDataFrameGroupBy.__rshift__ を定義する。

def _pipe(self, func):
    return func(self)

pd.core.groupby.DataFrameGroupBy.__rshift__ = _pipe

できた。

df >> groupby('Y') >> sum()
#    X
# Y   
# a  9
# b  6

df >> sum()
# X       15
# Y    ababa
# dtype: object

右辺側クラスによるオーバーライド

DataFrame が左辺にある場合、パイプ演算は上記のように定義できた。加えて、DataFrame が右辺にある場合も同じようなことがしたい。

以下のような記法で、左辺の辞書から DataFrame を作成したい。これも (当然) そのままではできない。

data = {'X': [1, 2, 3, 4, 5],
        'Y': ['a', 'b', 'a', 'b', 'a']} 

data >> pd.DataFrame()
# TypeError: unsupported operand type(s) for >>: 'dict' and 'DataFrame'

Python では、左辺側に右辺を処理できるマジックメソッドがない場合は、右辺側のマジックメソッド ( DataFarme.__rrshift__ ) が利用される。これを適当に定義すればよい。

def _get(a, b):
    return pd.DataFrame(b)

pd.DataFrame.__rrshift__ = _get

できた。

data = {'X': [1, 2, 3, 4, 5],
        'Y': ['a', 'b', 'a', 'b', 'a']} 

data >> pd.DataFrame()
#    X  Y
# 0  1  a
# 1  2  b
# 2  3  a
# 3  4  b
# 4  5  a

NumPy に対するオーバーライド

DataFrame は np.ndarray からも作成することができる。が、>> ではうまく動かない。

data = np.array([[1, 2], [3, 4]])

pd.DataFrame(data)
#    0  1
# 0  1  2
# 1  3  4

data >> pd.DataFrame()
# TypeError: ufunc 'right_shift' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

ndarray のマジックメソッドはスカラーや ndarray ライクな引数を汎用的に処理できるように定義されている。具体的には、引数が ndarray のサブクラス、もしくは Array Interface ( __array__ ) をもつとき、ndarray によって処理が行われる。

DataFrame は ndarray のサブクラスではないが Array Interface をもつ。右辺が不正なため通常の例外が送出されてしまい、右辺側のマジックメソッドに処理が移らない。

参考 Subclassing ndarray — NumPy v1.11 Manual

issubclass(pd.DataFrame, np.ndarray)
# False

df.__array__()
# array([[1, 'a'],
#        [2, 'b'],
#        [3, 'a'],
#        [4, 'b'],
#        [5, 'a']], dtype=object)

ここで DataFrame 側のマジックメソッドを利用したい場合、上記リンクに記載されている __array_priority__ プロパティを設定すればよい。DataFrame 側の優先度を上げると期待した処理が行われることがわかる。

pd.DataFrame.__array_priority__ = 1000

data >> pd.DataFrame()
#    0  1
# 0  1  2
# 1  3  4

まとめ

R 使ったほうがいい。Use R!

補足こんなパッケージもあります。

2016/10/18 追記

続き(?)です。

sinhrks.hatenablog.com

2015-07-26

Python xray で多次元データを pandas ライクに扱う

pandas Python 前処理

はじめに

pandas では 2 次元、表形式のデータ ( DataFrame ) を主な対象としているが、ときには 3 次元以上のデータを扱いたい場合がある。そういった場合以下のような方法がある。

MultiIndex を使い、2 次元のデータにマッピングする。
3 次元データ構造である Panel、4 次元の Panel4D、もしくは任意の次元のデータ構造 ( PanelND ) をファクトリ関数で定義して使う。
numpy.ndarray のまま扱う。

自分は MultiIndex を使うことが多いが、データを 2 次元にマップしなければならないため種類によっては直感的に扱いにくい。Panel や PanelND は DataFrame と比べると開発が活発でなく、特に Panel4D、PanelND は現時点で Experimental 扱いである。また、今後の扱いをどうするかも議論がある。numpy.ndarray ではデータのラベル付けができない。

`xray` とは

ラベル付きの多次元データを多次元のまま直感的に扱えるパッケージとして xray がある。作者は pandas 開発チーム仲間の shoyer だ。そのため、API は pandas にかなり近いものになっている。

xray は大きく以下ふたつのデータ構造を持つ。これらを使うと多次元データをより直感的に操作することができる。

xray.DataArray: numpy の多次元配列にラベルでのアクセスを追加したもの。データは任意の次元を持つことができる。
xray.Dataset: 複数の xray.DataArray をまとめるセット。

インストール

pip で。

$ pip install xray

データの準備

まず、必要なパッケージをインポートする。

import numpy as np
import pandas as pd

pd.__version__
# '0.16.2'

import xray
xray.__version__
# '0.5.2'

サンプルデータとして、気象庁から 2015年7月20日〜25日の東京、八王子、大島の最高気温のデータを使いたい。以下のサイトからダウンロードした。データは日付、場所 2 次元の配列となる。

出典：気象庁ホームページ (URL: http://www.data.jma.go.jp/gmd/risk/obsdl/index.php)

f:id:sinhrks:20150726225910p:plain

data2 = np.array([[34.2, 30.2, 33.5],
                  [36.0, 29.3, 34.9],
                  [35.3, 29.7, 32.8],
                  [30.1, 27.6, 30.4],
                  [33.6, 30.1, 33.9],
                  [34.1, 28.0, 33.1]])

この numpy.ndarray から xray.DataArray インスタンスを作成する。データへアクセスするためのラベル ( pandas でいう index のようなもの ) は coords キーワードで指定する。また、dims キーワードを使って各次元それぞれにも名前をつけることができる。ここでは以下のような DataArray を作成している。

データは 2 つの次元 date, location を持つ。
date 次元は dates で指定される 6 つの日付のラベルを持つ。
location 次元は locs で指定される 3 つの場所のラベルを持つ。

2 次元のため、pandas.DataFrame でいうと date が indexに、 location が columns に対応しているイメージ。

locs = [u'八王子', u'大島', u'東京']
dates = pd.date_range('2015-07-20', periods=6, freq='D')
da2 = xray.DataArray(data2, coords=[dates, locs], dims=['date', 'location'])
da2
# <xray.DataArray (date: 6, location: 3)>
# array([[ 34.2,  30.2,  33.5],
#        [ 36. ,  29.3,  34.9],
#        [ 35.3,  29.7,  32.8],
#        [ 30.1,  27.6,  30.4],
#        [ 33.6,  30.1,  33.9],
#        [ 34.1,  28. ,  33.1]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

# 各次元の名前
da2.dims
# ('date', 'location')

# 各次元の詳細
da2.coords
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

# date 次元のラベル
da2.coords['date']
# <xray.DataArray 'date' (date: 6)>
# array(['2015-07-20T09:00:00.000000000+0900',
#        '2015-07-21T09:00:00.000000000+0900',
#        '2015-07-22T09:00:00.000000000+0900',
#        '2015-07-23T09:00:00.000000000+0900',
#        '2015-07-24T09:00:00.000000000+0900',
#        '2015-07-25T09:00:00.000000000+0900'], dtype='datetime64[ns]')
# Coordinates:
#   * date     (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...

データの選択

xray でのデータ選択は pandas と類似の方法で行える。pandas でのデータ選択についてはこちらを。

Python pandas データ選択処理をちょっと詳しく <前編> - StatsFragments

ある次元から、特定のラベルをもつデータを選択したい場合、pandas と同じく .loc が使える。

# 7/25 のデータを選択
da2.loc[pd.Timestamp('2015-07-25')]
# <xray.DataArray (location: 3)>
# array([ 34.1,  28. ,  33.1])
# Coordinates:
#     date      datetime64[ns] 2015-07-25
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

また、対象の次元が datetime64 型の場合は、pandas と同じく日時文字列でも指定が可能 (部分文字列によるスライシング、詳細以下)。

pandas 日時まわりのリサンプリング/オフセット処理 - StatsFragments

# 7/25 のデータを選択
da2.loc['2015-07-25']
# <xray.DataArray (location: 3)>
# array([ 34.1,  28. ,  33.1])
# Coordinates:
#     date      datetime64[ns] 2015-07-25
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

また、ラベルではなく位置 ( n番目など ) によって選択する場合は __getitem__ を使う。pandas では .iloc に対応。

# 末尾 = 最新の日付のデータを取得
da2[-1]
# <xray.DataArray (location: 3)>
# array([ 34.1,  28. ,  33.1])
# Coordinates:
#     date      datetime64[ns] 2015-07-25
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

2 次元目以降をラベル / 位置によって指定する場合は、引数をカンマで区切り、選択に利用する次元の位置に対応する値を渡す。

# date は全選択、location が東京のデータを選択
da2.loc[:, u'東京']
# <xray.DataArray (date: 6)>
# array([ 33.5,  34.9,  32.8,  30.4,  33.9,  33.1])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location  <U3 u'\u6771\u4eac'

# date は全選択、location が 3 番目 = 東京のデータを選択
da2[:, 2]
# <xray.DataArray (date: 6)>
# array([ 33.5,  34.9,  32.8,  30.4,  33.9,  33.1])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location  <U3 u'\u6771\u4eac'

データ選択に利用する次元自体もラベルで指定したい場合は .sel を使う。.sel では次元の順序に関係なく値を指定できるため、高次元になった場合もシンプルだ。

# 東京のデータを選択
da2.sel(location=u'東京')
# <xray.DataArray (date: 6)>
# array([ 33.5,  34.9,  32.8,  30.4,  33.9,  33.1])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location  <U3 u'\u6771\u4eac'

# 7/22, 7/23 の 東京のデータを選択
da2.sel(location=u'東京', date=pd.DatetimeIndex(['2015-07-22', '2015-07-23']))
# <xray.DataArray (date: 2)>
# array([ 32.8,  30.4])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-22 2015-07-23
#     location  <U3 u'\u6771\u4eac'

もしくは .loc, __getitem__ に以下のような辞書を渡してもよい。

# 東京のデータを選択
da2.loc[dict(location=u'東京')]
# <xray.DataArray (date: 6)>
# array([ 33.5,  34.9,  32.8,  30.4,  33.9,  33.1])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location  <U3 u'\u6771\u4eac'

公式ドキュメントでは、xray.DataArray からのデータ選択を以下のような表で整理している。

次元の指定	値の指定	`DataArray` のメソッド
位置	位置	`da2[:, 0]`
位置	ラベル	`da2.loc[:, u'東京']`
ラベル	位置	`da2.isel(location=0)`, `da2[dict(space=0)]`
ラベル	ラベル	`da2.sel(location=u'東京')`, `da2.loc[dict(location=u'東京')]`

次元の追加

ここまでは 2 次元のデータを使っていたが、最低気温と平均気温のデータを追加して 3 次元のデータとする。

# 最低気温
data3 = np.array([[24.7, 25.1, 25.8],
                  [23.4, 24.3, 25.4],
                  [22.7, 23.8, 25.4],
                  [24.5, 24.4, 24.7],
                  [24.9, 24.6, 25.0],
                  [23.7, 24.8, 24.8]])
# 平均気温
data4 = np.array([[27.8, 26.5, 28.8],
                  [29.7, 26.3, 29.4],
                  [29.5, 26.5, 28.9],
                  [26.9, 25.3, 27.0],
                  [27.7, 26.4, 28.1],
                  [29.0, 26.1, 28.6]])
data = np.dstack([data2, data3, data4])

# 日時、場所、データの種類 の 3 次元
data.shape
# (6, 3, 3)

da3 = xray.DataArray(data, coords=[dates, locs, [u'最高', u'最低', u'平均']],
                     dims=['date', 'location', 'type'])
da3
# <xray.DataArray (date: 6, location: 3, type: 3)>
# array([[[ 34.2,  24.7,  27.8],
#         [ 30.2,  25.1,  26.5],
#         [ 33.5,  25.8,  28.8]],
# 
#      ... 
# 
#        [[ 34.1,  23.7,  29. ],
#         [ 28. ,  24.8,  26.1],
#         [ 33.1,  24.8,  28.6]]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...
#   * type      (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'

3 次元以上の場合もデータ選択のルールは一緒なのでわかりやすい。

# 7/24 のデータを選択
da3.loc['2015-07-24']
# <xray.DataArray (location: 3, type: 3)>
# array([[ 33.6,  24.9,  27.7],
#        [ 30.1,  24.6,  26.4],
#        [ 33.9,  25. ,  28.1]])
# Coordinates:
#     date      datetime64[ns] 2015-07-24
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...
#   * type      (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'

# 7/24 の 最高気温のデータを選択
da3.loc['2015-07-24', :, u'最高']
# <xray.DataArray (location: 3)>
# array([ 33.6,  30.1,  33.9])
# Coordinates:
#     date      datetime64[ns] 2015-07-24
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...
#     type      <U2 u'\u6700\u9ad8'

# 東京 の 最高気温のデータを選択
da3.sel(location=u'東京', type=u'最高')
# <xray.DataArray (date: 6)>
# array([ 33.5,  34.9,  32.8,  30.4,  33.9,  33.1])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location  <U3 u'\u6771\u4eac'
#     type      <U2 u'\u6700\u9ad8'

算術演算

pandas と同じく、xray.DataArray 同士での算術演算が可能。各日の最高気温と最低気温の差を求めると、

da3.sel(type=u'最高') - da3.sel(type=u'最低') 
# <xray.DataArray (date: 6, location: 3)>
# array([[  9.5,   5.1,   7.7],
#        [ 12.6,   5. ,   9.5],
#        [ 12.6,   5.9,   7.4],
#        [  5.6,   3.2,   5.7],
#        [  8.7,   5.5,   8.9],
#        [ 10.4,   3.2,   8.3]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

# 結果から 八王子 のデータだけを選択
(da3.sel(type=u'最高') - da3.sel(type=u'最低')).sel(location=u'八王子')
# <xray.DataArray (date: 6)>
# array([  9.5,  12.6,  12.6,   5.6,   8.7,  10.4])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location  <U3 u'\u516b\u738b\u5b50'

データのグループ化 / 集約

グループ化 / 集約も pandas とほぼ同じ形式でできる。location によってグループ化し、期間中の最高気温を出してみる。

da3.groupby('location')
# <xray.core.groupby.DataArrayGroupBy at 0x109a56f90>

da3.groupby('location').max()
# <xray.DataArray (location: 3)>
# array([ 36. ,  30.2,  34.9])
# Coordinates:
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

グループ化した結果は、イテレーションによって順番に処理することもできる。

for name, g in da3.groupby('location'):
    print(name)
    print(g)
# 八王子
# <xray.DataArray (date: 6, type: 3)>
# array([[ 34.2,  24.7,  27.8],
#        [ 36. ,  23.4,  29.7],
#        [ 35.3,  22.7,  29.5],
#        [ 30.1,  24.5,  26.9],
#        [ 33.6,  24.9,  27.7],
#        [ 34.1,  23.7,  29. ]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * type      (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'
#     location  <U4 u'\u516b\u738b\u5b50'
# 
# 以降略

データの結合 / 連結

上で作成した DataArray に、別のデータを追加したい。適当なデータを探したところ、ある成人男性のお住まいの気温データを見つけた。これを日別で集計して連結したい。

# GitHub からデータを取得
df = pd.read_csv('https://raw.githubusercontent.com/dichika/mydata/master/room.csv')
df['time'] = pd.to_datetime(df['time'])

# light は光量、temperatur は気温
df.head()
#                  time  light  temperature
# 0 2015-02-12 01:45:04  463.0       26.784
# 1 2015-02-12 02:00:03  473.0       26.630
# 2 2015-02-12 02:15:04  467.9       25.983
# 3 2015-02-12 02:30:04    0.0       25.453
# 4 2015-02-12 02:45:04    0.0       23.650

# 期間中にフィルタ
df = df[df['time'] >= pd.Timestamp('2015-07-20')]
# 日時でグループ化 / 集約
agg = df.groupby(pd.Grouper(key='time', freq='D'))['temperature'].agg(['max', 'min', 'mean'])
agg
#                max     min       mean
# time                                 
# 2015-07-20  28.374  26.450  27.484604
# 2015-07-21  33.790  26.800  30.132792
# 2015-07-22  34.180  27.070  30.134375
# 2015-07-23  28.779  27.412  28.290918

xray でデータを連結するためには、連結する次元 (ここでは location) 以外のデータの要素数を一致させる必要がある。上記のデータは数日遅れで公開されているため、直近を NaN でパディングして xray.DataArray を作成する。

agg.loc[pd.Timestamp('2015-07-24'), :] = np.nan
agg.loc[pd.Timestamp('2015-07-25'), :] = np.nan
agg
#                max     min       mean
# time                                 
# 2015-07-20  28.374  26.450  27.484604
# 2015-07-21  33.790  26.800  30.132792
# 2015-07-22  34.180  27.070  30.134375
# 2015-07-23  28.779  27.412  28.290918
# 2015-07-24     NaN     NaN        NaN
# 2015-07-25     NaN     NaN        NaN

d = xray.DataArray(agg.values.reshape(6, 1, 3), coords=[agg.index, [u'誰かの家'], [u'最高', u'最低', u'平均']],
                   dims=['date', 'location', 'type'])
d
# <xray.DataArray (date: 6, location: 1, type: 3)>
# array([[[ 28.374     ,  26.45      ,  27.48460417]],
#        [[ 33.79      ,  26.8       ,  30.13279167]],
#        [[ 34.18      ,  27.07      ,  30.134375  ]],
#        [[ 28.779     ,  27.412     ,  28.29091765]],
#        [[         nan,          nan,          nan]],
#        [[         nan,          nan,          nan]]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location  (location) <U4 u'\u8ab0\u304b\u306e\u5bb6'
#   * type      (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'

データの連結は xray.concat で可能。新しい場所のデータを追加 (連結) したいので、対象の次元として location を指定する。

da3 = xray.concat([da3, d], dim='location')
da3.sel(location=u'誰かの家')
# <xray.DataArray (date: 6, type: 3)>
# array([[ 28.374     ,  26.45      ,  27.48460417],
#        [ 33.79      ,  26.8       ,  30.13279167],
#        [ 34.18      ,  27.07      ,  30.134375  ],
#        [ 28.779     ,  27.412     ,  28.29091765],
#        [         nan,          nan,          nan],
#        [         nan,          nan,          nan]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * type      (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'
#     location  <U4 u'\u8ab0\u304b\u306e\u5bb6'

ほか、merge で結合もできる。

`Dataset` の利用

xray では、DataArray クラスを複数まとめて Dataset クラスとして扱うことができる。元となる DataArray は同じ次元でなくてもよい。

Dataset を作成するため、降水量、湿度それぞれ 2 次元の DataArray を用意する。

# 降水量
precip = np.array([[0, np.nan, 0],
                   [0, np.nan, np.nan],
                   [0, 0, np.nan],
                   [6.5, 1, 4.5],
                   [30.0, 0, 7.0],
                   [0, np.nan, np.nan]])
precip = xray.DataArray(precip, coords=[dates, locs], dims=['date', 'location'])
precip
# <xray.DataArray (date: 6, location: 3)>
# array([[  0. ,   nan,   0. ],
#        [  0. ,   nan,   nan],
#        [  0. ,   0. ,   nan],
#        [  6.5,   1. ,   4.5],
#        [ 30. ,   0. ,   7. ],
#        [  0. ,   nan,   nan]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

# 湿度
humid = np.array([[np.nan, 88, 75],
                  [np.nan, 85, 65],
                  [np.nan, 87, 61],
                  [np.nan, 91, 80],
                  [np.nan, 86, 83],
                  [np.nan, 88, 80]])
humid = xray.DataArray(humid, coords=[dates, locs], dims=['date', 'location'])
humid
# <xray.DataArray (date: 6, location: 3)>
# array([[ nan,  88.,  75.],
#        [ nan,  85.,  65.],
#        [ nan,  87.,  61.],
#        [ nan,  91.,  80.],
#        [ nan,  86.,  83.],
#        [ nan,  88.,  80.]])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location  (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...

これらと気温データをあわせて xray.Dataset を作成する。Dataset が持つ次元は Dimenstions, Coordinates に表示される。 Dataset に含まれる DataArray は Data variables 中に表示され、それぞれどの次元を含んでいるかがわかる。元データと同じく、気温は 3 次元、ほかは 2 次元のデータとなっている。

ds = xray.Dataset({'temperature': da3,
                   'precipitation': precip,
                   'humidity': humid})
ds
# <xray.Dataset>
# Dimensions:        (date: 6, location: 4, type: 3)
# Coordinates:
#   * date           (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#   * location       (location) object u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...
#   * type           (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'
# Data variables:
#     precipitation  (date, location) float64 0.0 nan 0.0 nan 0.0 nan nan nan ...
#     temperature    (date, location, type) float64 34.2 24.7 27.8 30.2 25.1 ...
#     humidity       (date, location) float64 nan 88.0 75.0 nan nan 85.0 65.0 ...

ds.data_vars
# Data variables:
#     precipitation  (date, location) float64 0.0 nan 0.0 nan 0.0 nan nan nan ...
#     temperature    (date, location, type) float64 34.2 24.7 27.8 30.2 25.1 ...
#     humidity       (date, location) float64 nan 88.0 75.0 nan nan 85.0 65.0 ...

Dataset からのデータ選択は、Dataset に含まれるすべての DataArray に対して行われる。DataArray とは異なり、次元は必ずラベルで指定する必要がある。

次元の指定	値の指定	`Dataset` のメソッド
位置	位置	なし
位置	ラベル	なし
ラベル	位置	`ds.isel(location=2)`, `ds[dict(location=2)]`
ラベル	ラベル	`ds.sel(location=u'東京')`, `ds.loc[dict(location=u'東京')]`

# 東京 のデータを選択
ds.sel(location=u'東京')
# <xray.Dataset>
# Dimensions:        (date: 6, type: 3)
# Coordinates:
#   * date           (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location       object u'\u6771\u4eac'
#   * type           (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'
# Data variables:
#     precipitation  (date) float64 0.0 nan nan 4.5 7.0 nan
#     temperature    (date, type) float64 33.5 25.8 28.8 34.9 25.4 29.4 32.8 ...
#     humidity       (date) float64 75.0 65.0 61.0 80.0 83.0 80.0

# 東京 の 平均気温 を選択
ds.sel(location=u'東京', type=u'平均')['temperature']
# <xray.DataArray 'temperature' (date: 6)>
# array([ 28.8,  29.4,  28.9,  27. ,  28.1,  28.6])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     type      <U2 u'\u5e73\u5747'
#     location  object u'\u6771\u4eac'

# 東京 の 湿度 を選択
ds.sel(location=u'東京')['humidity']
# <xray.DataArray 'humidity' (date: 6)>
# array([ 75.,  65.,  61.,  80.,  83.,  80.])
# Coordinates:
#   * date      (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location  object u'\u6771\u4eac'

# 対応する DataArray がない場合は NaN となる
ds.sel(location=u'誰かの家')
# <xray.Dataset>
# Dimensions:        (date: 6, type: 3)
# Coordinates:
#   * date           (date) datetime64[ns] 2015-07-20 2015-07-21 2015-07-22 ...
#     location       object u'\u8ab0\u304b\u306e\u5bb6'
#   * type           (type) <U2 u'\u6700\u9ad8' u'\u6700\u4f4e' u'\u5e73\u5747'
# Data variables:
#     precipitation  (date) float64 nan nan nan nan nan nan
#     temperature    (date, type) float64 28.37 26.45 27.48 33.79 26.8 30.13 ...
#     humidity       (date) float64 nan nan nan nan nan nan

そのほかの操作も DataArray と同じようにできる。

# location ごとに期間中の最大値を計算
ds.groupby('location').max()
# <xray.Dataset>
# Dimensions:        (location: 3)
# Coordinates:
#   * location       (location) <U3 u'\u516b\u738b\u5b50' u'\u5927\u5cf6' ...
# Data variables:
#     precipitation  (location) float64 30.0 1.0 7.0
#     temperature    (location) float64 36.0 30.2 34.9
#     humidity       (location) float64 nan 91.0 83.0

`pandas` のデータ形式への変換

DataArray, Dataset はそれぞれ pandas のデータに変換できる。元データが 3 次元以上の場合、変換後の pandas のデータは MultiIndex を持つことになる。

また、上の例では次元が異なる DataArray から Dataset を作成した。このとき、存在しない次元 ( ここでは type ) のデータはすべて同じ値でパディングされる。

ds.to_dataframe()

f:id:sinhrks:20150726231436p:plain

パディングしたくない場合は、個々の DataArray ごとに DataFrame に変換すればよい。

ds['humidity'].to_dataframe()

f:id:sinhrks:20150726231444p:plain

まとめ

xray を使えば多次元のラベル付きデータを多次元のまま、pandas に近い方法で扱うことができる。

2015-07-18

Python geopandas + Bokeh で地理情報をプロットしたい

pandas Python 可視化

数日前、pandas を利用して地理情報をプロットするという非常によいエントリが翻訳されていた。

postd.cc

上のエントリ、前処理が手間に見えるが pd.read_html や .str アクセサを使えばもっと簡単に書けると思う、、、がそれは本題でない。

pandas で地理情報を扱う場合、geopandas という拡張パッケージを利用すると便利なため、その使い方を書きたい。また、処理を Python で完結させるため、QGIS ではなく Bokeh でプロットしたい。

`geopandas` のインストール

pip で。

$ pip install geopandas geopy

このエントリでは依存パッケージである shapely、geopy の機能も利用する。shapely は自動的にインストールされるはずだが、geopy については上のように別途インストールが必要。

地理情報の読み込み

以降の操作は IPython Notebook で行う。まずは必要なパッケージをインポートする。

%matplotlib inline
import numpy as np
import pandas as pd
import geopandas as gpd

例として、東京都の地図を読み込んで利用したい。geopandas は、依存パッケージである fiona がサポートする地理情報ファイル形式を読み込むことができる。GitHub に地球地図日本のデータから作成した GeoJSON ファイルがあったため、こちらを利用させていただく。

dataofjapan/land · GitHub
- "出典: 地球地図日本(国土地理院技術資料D・1-No.576)" を以下記載のスクリプトで加工して利用。

ファイルの読み込みは、geopandas.read_file で行う。データは pandas.DataFrame を継承した geopandas.GeoDataFrame インスタンスとして読み込まれる。

df = gpd.read_file('tokyo.geojson')

type(df)
# geopandas.geodataframe.GeoDataFrame

isinstance(df, pd.DataFrame)
# True

df.head()
#   area_en area_ja    code                                           geometry  
# 0  Tokubu     都区部  131211  POLYGON ((139.821051 35.815077, 139.821684 35....   
# 1  Tokubu     都区部  131059  POLYGON ((139.760933 35.732206, 139.761002 35....   
# 2  Tokubu     都区部  131016  POLYGON ((139.770135 35.705352, 139.770172 35....   
# 3  Tokubu     都区部  131067  POLYGON ((139.809714 35.728135, 139.809705 35....   
# 4  Tokubu     都区部  131091  (POLYGON ((139.719199 35.641847, 139.719346 35...

GeoDataFrame.plot を使うと、データを matplotlib を利用してプロットできる。

df.plot()

f:id:sinhrks:20150718210549p:plain

島嶼部があると細部がわからないため、都区部だけをフィルタしてプロットする。

df[df['area_en'] == 'Tokubu'].plot()

f:id:sinhrks:20150718210600p:plain

地理情報の処理

GeoDataFrame の各行はそれぞれひとつの領域 (ジオメトリオブジェクト) に対応している。各領域のジオメトリオブジェクトの実体は geometry キー中に shapely の Polygon もしくは MultiPolygon インスタンスとして保存されている。

Polygon には、領域各点の経度/緯度からなる tuple がリストとして保存されている。以下の処理で領域各点の経度/緯度をそれぞれ別のリストとして取り出すことができる。この処理は、続けて記載する Bokeh でのプロット時に必要になる。

# 1行目 (1領域目) を取得
s = df.iloc[0, :]
s
# area_en                                                Tokubu
# area_ja                                                   都区部
# code                                                   131211
# geometry    POLYGON ((139.821051 35.815077, 139.821684 35....
# ward_en                                             Adachi Ku
# ward_ja                                                   足立区
# Name: 0, dtype: object

type(s)
# <class 'pandas.core.series.Series'>

# 領域のジオメトリオブジェクトの型を表示
type(s['geometry'])
# <class 'shapely.geometry.polygon.Polygon'>

# 領域の経度/緯度を取得
list(s['geometry'].exterior.coords[:5])
# [(139.821051, 35.815077), (139.821684, 35.814887), (139.822599, 35.814509),
#  (139.8229070000001, 35.81437), (139.822975, 35.814339)]

# 経度、緯度別のリストに分割
lons, lats = zip(*list(s['geometry'].exterior.coords))

# 経度を表示
lons[:5]
# (139.821051, 139.821684, 139.822599, 139.8229070000001, 139.822975, 139.822977)

地理情報のプロット

このデータを Bokeh でプロットしたい。Bokeh についてはこちらを。

Chainer で Deep Learning: Bokeh で Live Monitoring したい - StatsFragments

Bokeh で地理情報をプロットするサンプルとして、ギャラリーにコロプレス図を描く例がある。Bokeh では描画する各領域の x 座標 (経度)、y 座標 (緯度) をリストで指定する。各領域は複数座標からなるため、"各領域の x/y 座標のリスト" のリスト ( 例: [[境界1の点1のx座標, 境界1の点2のx座標, ...], [領域2の点1のx座標, ...], ...] ) を引数として用意する必要がある。

import shapely

# "各領域の x 座標 (経度) のリスト" のリスト
countries_xs = []

# "各領域の y 座標 (緯度) のリスト" のリスト
countries_ys = []

for i, row in df.iterrows():
    if row['area_en'] != 'Tokubu':
        # 都区部のみ描画
        continue

    polygons = row['geometry']
    if isinstance(polygons, shapely.geometry.multipolygon.MultiPolygon):
        # MultiPolygon に含まれる Polygon を取り出す
        polygons = [x for x in polygons]
    elif isinstance(polygons, shapely.geometry.polygon.Polygon):
        polygons = [polygons]
    else:
        raise ValueError
        
    for p in polygons: 
        # 境界の各点を経度、緯度に分割
        lons, lats = zip(*list(p.exterior.coords))
        countries_xs.append(lons)
        countries_ys.append(lats)

作成したデータを Bokeh に渡す。出力先は IPython Notebook にした。

from bokeh import plotting

# 出力先を指定
plotting.output_notebook()

# figure を作成
p = plotting.figure(toolbar_location="left", plot_width=900, plot_height=700)

# 各領域を描画 (塗り分けはしない)
p.patches(countries_xs, countries_ys, fill_color='white',
          line_color="black", line_width=0.5)
plotting.show(p)

f:id:sinhrks:20150718210612p:plain

Geocoding の利用

geopandas は geopy を利用して Geocoding (住所から緯度経度への変換) を行うための API geopandas.geocode.geocode を持っている。が、geopandas v0.1.1, geopy v1.10.0 では、これらのインターフェースに不整合がありそのままでは利用できない。詳細と回避策は以下 Stack Overflow を。

python - Geocoding error with geopandas and geopy - Stack Overflow

from geopandas.geocode import geocode

# geopandas.geocode の API をそのまま使うとエラー
geocode(['千代田区千代田1-1'])
# AttributeError: 'module' object has no attribute 'MapQuest'

# 回避策として geopy の以下プロパティを明示的に None にする
import geopy
geopy.geocoders.MapQuest = None
geocode(['千代田区千代田1-1'])
#                                              address                        geometry
# 0  1-1 Chiyoda, Chiyoda-ku, Tōkyō-to 100-0001, Japan  POINT (139.7539454 35.6838012)

Geocoding を利用したデータの追加

上で作成した地図に、道路交通センサスにある各地点の交通量をプロットしたい。もっと細かいデータが見られるとよいのだが、公開情報ではこれくらいしかなさそうだ。以下のリンクから東京都の箇所別基本表 kasyo13.xls をダウンロードし、利用させていただく。

平成22年度道路交通センサス
- 出典：「平成22年度　全国道路・街路交通情勢調査（道路交通センサス）一般交通量調査　集計表」（国土交通省）

pdf = pd.read_excel('kasyo13.xls', header=8, index_col=False)
pdf.head()
#           (km)  Unnamed: 1  Unnamed: 2   Unnamed: 3  Unnamed: 4   Unnamed: 5  
# 0  13110100010           1        1010       東名高速道路           1  14110100090   
# 1  13111010010           1        1101  中央自動車道富士吉田線           4  13200400170   
# 2  13111010020           1        1101  中央自動車道富士吉田線           6          NaN   
# 3  13111010030           1        1101  中央自動車道富士吉田線           6          NaN   
# 4  13111010040           1        1101  中央自動車道富士吉田線           6          NaN

# ヘッダがセル結合されており正しく読めないため、連番で上書き
pdf.columns = range(len(pdf.columns))

# 計測地点の住所が入力されている列を表示
pdf.iloc[:, 28].head()
# 0             東京～東名川崎
# 1    高井戸ＩＣ～世田谷区・東京都境間
# 2    高井戸ＩＣ～世田谷区・東京都境間
# 3     世田谷区・東京都境～調布ＩＣ間
# 4    高井戸ＩＣ～世田谷区・東京都境間
# Name: 28, dtype: object

# NaN, 正しい住所でないデータをフィルタ
pdf = pdf.dropna(subset=[28])
pdf = pdf[pdf[28].str.startswith('東京都')]

# 計測地点住所と 24時間交通量の列のみにフィルタ
pdf = pdf.iloc[:, [28, 36]]
pdf.columns = ['計測地点', '交通量合計']
pdf.head()
#             計測地点   交通量合計
# 22  東京都中央区八丁堀３丁目  107885
# 23  東京都中央区八丁堀１丁目   98482
# 24   東京都中央区新富１丁目   89089
# 25   東京都中央区銀座１丁目   94143
# 26   東京都中央区築地１丁目   93844

上のデータから、計測地点 カラムの住所をもとに、Geocoding で緯度経度を取得する。

codes = geocode(pdf['計測地点'])

type(codes)
# <class 'geopandas.geodataframe.GeoDataFrame'>

codes.head()
#                                               address                        geometry
# 22  3 Chome Hatchōbori, Chūō-ku, Tōkyō-to 104-0032...  POINT (139.7756944 35.6755445) 
# 23  1 Chome Hatchōbori, Chūō-ku, Tōkyō-to 104-0032...  POINT (139.7771192 35.6774165)  
# 24  1 Chome Shintomi, Chūō-ku, Tōkyō-to 104-0041, ...  POINT (139.7725002 35.6729586) 
# 25  Ginza Itchome Station, 1 Chome-6 Ginza, Chūō-k...     POINT (139.767045 35.67435) 
# 26  1 Chome Tsukiji, Chūō-ku, Tōkyō-to 104-0045, J...  POINT (139.7716809 35.6702089)

取得したデータは上記のような GeoDataFrame となる。geometry カラムに経度/緯度を示す Point インスタンスが含まれる。この座標と交通量を利用してバブルチャートが描きたい。Bokeh のギャラリーにバブルチャートをプロットする例があるため、それに倣うと、

pdf = pdf.join(codes)

# 各点の x 座標
points_x = []

# 各点の y 座標
points_y = []

# 各点の円の半径
points_v = []

for i, row in pdf.iterrows():
    points_x.append(row['geometry'].x)
    points_y.append(row['geometry'].y)
    points_v.append(row['交通量合計'])
# 円の半径を適当に調整
points_v = np.sqrt(np.array(points_v)) / 50000

# Bokeh での描画
plotting.output_notebook()
p = plotting.figure(toolbar_location="left", plot_width=900, plot_height=700)
p.patches(countries_xs, countries_ys, fill_color='white',
          line_color="black", line_width=0.5)
# バブルチャートを描画
p.scatter(points_x, points_y, radius=points_v, fill_alpha=0.1)
plotting.show(p)

f:id:sinhrks:20150718210626p:plain

まとめ

geopandas を使って、地理情報に対する以下のような処理を行う方法を記載した。

各種地理情報フォーマットの読み込み/書き込み
各領域 (ジオメトリオブジェクト) からの座標 (緯度/経度) の取得
matplotlib, Bokeh でのプロット
Geocoding

公式ドキュメントには、上以外にジオメトリオブジェクトを操作 / 変換する例も記載されている。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
出版社/メーカー: オライリージャパン
発売日: 2013/12/26
メディア: 大型本
この商品を含むブログ (12件) を見る

2015-07-11

Python pandas パフォーマンス維持のための 3 つの TIPS

pandas Python

pandas でそこそこ大きいデータを扱う場合、その処理速度が気になってくる。公式ドキュメントではパフォーマンス向上のために Cython や Numba を使う方法を記載している。

Enhancing Performance — pandas 0.16.2 documentation

が、軽く試したいだけなのにわざわざ Cythonや Numba を使うのは手間だし、かといってあまりに遅いのも嫌だ。そんなとき、pandas 本来のパフォーマンスをできるだけ維持するためのポイントを整理したい。

pandas に限らず、パフォーマンス改善の際にはボトルネックの箇所によってとるべき対策は異なる。pandas では速度向上/エッジケース処理のためにデータの型や条件によって内部で処理を細かく分けており、常にこうすれば速くなる！という方法を出すのは難しい。以下はこの前提のうえで、内部実装からみてまあほとんどの場合あてはまるだろう、、、という内容をまとめる。

環境構築

環境は EC2 の c4.xlarge インスタンス上に作成する。パフォーマンスを重視する場合、環境構築の時点で以下を行っておく。

numpy を各種数値計算ライブラリ BLAS, LAPACK, ATLAS とリンクさせる。
numexpr, bottleneck をインストールする。

EC2 の Amazon Linux であれば以下のコマンドでできる。

$ sudo yum -y install blas-devel lapack-devel atlas-devel 
$ sudo python -m pip install numpy
$ sudo python -m pip install numexpr bottleneck

環境構築が正しくできたかを確認する。以降の処理は IPython Notebook 上で行う。

import numpy as np
import numpy.distutils.system_info as sysinfo

sysinfo.get_info('blas')
# {'libraries': ['blas'], 'library_dirs': ['/usr/lib64'], 'language': 'f77'}

sysinfo.get_info('lapack')
# {'libraries': ['lapack'], 'library_dirs': ['/usr/lib64'], 'language': 'f77'}

sysinfo.get_info('atlas')
# {'define_macros': [('ATLAS_INFO', '"\\"3.8.4\\""')],
#  'include_dirs': ['/usr/include'],
#  'language': 'f77',
#  'libraries': ['lapack', 'f77blas', 'cblas', 'atlas'],
#  'library_dirs': ['/usr/lib64/atlas']}

import pandas as pd
pd.show_versions()
# ...
# pandas: 0.16.2
# ...
# bottleneck: 1.0.0
# ...
# numexpr: 2.4.3
# ...

サンプルデータの作成

サンプルデータ作成には pandas がテスト用に用意しているメソッドを使う。もっとも、似たようなデータが作れれば方法はなんでもいい。

pd.util.testing.rands_array: 指定した文字数 / 要素数のランダムな文字列の np.ndarray を作成。
pd.util.testing.choice: np.ndarray から size 個をサンプリング。

# 2 文字 / 3 要素の array を作成
pd.util.testing.rands_array(2, 3)
# array(['Xn', 'ZC', 'zj'], dtype=object)

# array から 5 回サンプリング
pd.util.testing.choice(np.array(['a', 'b']), size=5)
# array(['a', 'b', 'b', 'b', 'a'], dtype='|S1')

サンプルデータとして 100 万レコード / 3 列のデータを用意した。

np.random.seed(0)

# レコード数
N = 1000000

chars1 = pd.util.testing.rands_array(5, 100)
chars2 = pd.util.testing.rands_array(5, 10000)
df = pd.DataFrame({'x': np.random.randn(N),
                   'y': pd.util.testing.choice(chars1, size=N),
                   'z': pd.util.testing.choice(chars2, size=N)})
df.shape
# (1000000, 3)

df.head()

f:id:sinhrks:20150711222824p:plain

y 列は 100 通りの値を、z 列は 10000 通りの値をとる。

chars1[:10]
# array(['SV1ad', '7dNjt', 'vYKxg', 'yym6b', 'MNxUy', 'rLzni', 'juZqZ',
#        'fpVas', 'JyXZD', 'ttoNG'], dtype=object)

len(np.unique(chars1))
# 100

len(np.unique(chars2))
# 10000

1. 行に対するループ / `DataFrame.apply` は使わない

pandas.DataFrame は列ごとに異なる型を持つことができる。DataFrame は内部的に同じ型の列をまとめて np.ndarray として保持している。列ごとに連続したデータを持つことになるため、そもそも行に対するループには向かない。また、DataFrame.iterrows でのループの際には異なる型を持つ列の値を Series として取り出すため、そのインスタンス化にも時間がかかる。

また行ごと / 列ごとに関数を適用するメソッドに DataFrame.apply があるが、このメソッドでは Python の関数を繰り返し呼び出すためのコストがかかる。apply は利便性を重視したメソッドのため、パフォーマンスを気にする場合は避けたほうがよい。

参考 Python pandas データのイテレーションと関数適用、pipe - StatsFragments

上記ふたつの処理の組み合わせである、各行への関数適用 DataFrame.apply(axis=1) について処理時間を %timeit で計測する。まずは単純に y 列と z 列の値を文字列結合する場合。

def f1(s):
    return s['y'] + s['z']

%timeit df.apply(f1, axis=1)
# 1 loops, best of 3: 13.7 s per loop

この処理は、Series として列を取り出しベクトルとして行うほうが格段に速い。

%timeit df['y'] + df['z']
# 10 loops, best of 3: 74.9 ms per loop

ただ、関数中に条件分岐を含む場合など、そのままベクトル化できない場合もある。x 列の値によって結合の順序を変える例を考える。

def f2_1(s):
    if s['x'] > 0:
        return s['y'] + s['z']
    else:
        return s['z'] + s['y']

%timeit df.apply(f2_1, axis=1)
# 1 loops, best of 3: 16 s per loop

こういった場合は np.vectorize で関数をベクトル化し、引数として各列を Series (もしくは np.ndarray ) として渡したほうが apply よりは速い。

def f2_2(x, y, z):
    if x > 0:
        return y + z
    else:
        return z + y

%timeit pd.Series(np.vectorize(f2_2)(df['x'], df['y'], df['z']), index=df.index)
# 1 loops, best of 3: 334 ms per loop

補足 numpy のドキュメントに記載されているとおり、np.vectorize もパフォーマンスを最重視した方法ではない。さらに速くしたい場合は個別にベクトル化した関数を用意する。参考として、Cython で書いた場合の処理時間は以下となる。np.vectorize の 4-5 倍は速いようだ。

%load_ext cython

%%cython

import numpy as np
from numpy cimport ndarray

def f2_3(ndarray[double, ndim=1] x,
         ndarray[object, ndim=1] y,
         ndarray[object, ndim=1] z):
    
    cdef:
        int i, length = len(x)
        double xval
        object yval, zval
        ndarray[object, ndim=1] result = np.empty(length, dtype=object)
        
    for i in range(length):
        xval = x[i]
        yval = y[i]
        zval = z[i]
        if xval > 0:
            result[i] = yval + zval
        else:
            result[i] = zval + yval
            
    return result

%timeit pd.Series(f2_3(df['x'].values, df['y'].values, df['z'].values), index=df.index)
# 10 loops, best of 3: 67.5 ms per loop

補足 Cython については最近以下の書籍が出ていた。翻訳はわからないが、原著はわかりやすく要点がまとまっていたのでおすすめ。

Cython ―Cとの融合によるPythonの高速化

作者: Kurt W. Smith,中田秀基,長尾高弘
出版社/メーカー: オライリージャパン
発売日: 2015/06/19
メディア: 大型本
この商品を含むブログ (3件) を見る

2. `object` 型は使わない

pandas には、グループ化 ( DataFrame.groupby ) や結合 ( DataFrame.merge ) など、データの値をキーにして行われる処理がある。こういった操作を行う場合は object 型を避けたほうが速くなることが多い。

pandas では文字列を含む列は object 型になることに注意する。numpy の文字列型とならない理由は簡単にいうと欠損値 ( NaN ) の処理のため。例えば文字列で保存された商品コードをキーにしてグループ化 / 結合する場合など、そのままでは object 型として処理されてしまう。

df.dtypes
# x    float64
# y     object
# z     object
# dtype: object

object 型の列をキーとしてグループ化 / 集約したときの処理時間は以下。

%timeit df.groupby('y').mean()
# 10 loops, best of 3: 52.3 ms per loop

処理速度をあげるためには、キーとなる値をカテゴリ型 ( pd.Categorical ) に変換しておくとよい。これは R でいう factor にあたる型。

df['y'] = df['y'].astype('category')
%timeit df.groupby('y').mean()
# 100 loops, best of 3: 6.89 ms per loop

カテゴリ型では、カテゴリのラベル (pd.Categorical.categories ) とラベルの位置 (pd.Categorical.codes ) を分けて扱う。内部処理は int 型で保存されたラベルの位置について行われるため、object に対する処理と比較すると速くなる。

# カテゴリ型の内部データを表示
df['y'].values
# [ewJ6t, JaFfE, ttoNG, F2Sy9, OopuJ, ..., dS9oG, juZqZ, 5gvFn, 9brWJ, 9fxRG]
# Length: 1000000
# Categories (100, object): [0dmK0, 1DdJN, 1IZE1, 1m5Qu, ..., x7c5I, ydsVd, yym6b, zYGBj]

# カテゴリのラベル
df['y'].cat.categories
# Index([u'0dmK0', u'1DdJN', u'1IZE1', u'1m5Qu', u'219tH', u'2aMtU', u'403al',
#        ...
#        u'yym6b', u'zYGBj'],
#       dtype='object')

# カテゴリのラベルの位置
df['y'].cat.codes.head()
# 0    66
# 1    30
# 2    90
# 3    23
# 4    38
# dtype: int8

列がとりうる値の数が増えた場合も、カテゴリ型にしたほうが速い。

# object 型
%timeit df.groupby('z').mean()
# 10 loops, best of 3: 78.4 ms per loop

# カテゴリ型
df['z'] = df['z'].astype('category')
%timeit df.groupby('z').mean()
# 100 loops, best of 3: 17.9 ms per loop

列の型 dtype 以外に処理に影響を与えうるのは、欠損値 ( NaN ) の有無 (無いほうが速い)、object 型内での型の混在 (文字列と数値が混ざっている) など。

3. ユニークでない / ソートされていない `index` は使わない

最後。index についても、上のとおり object 型は避ける。また、特に理由がない場合値をユニークにし、かつソートしておく。上記のサンプルデータ作成時点では、特に index を指定していないため int 型の index が昇順で振られている。

df.index
# Int64Index([     0,      1,      2,      3,      4,      5,      6,      7,
#                  8,      9, 
#             ...
#             999990, 999991, 999992, 999993, 999994, 999995, 999996, 999997,
#             999998, 999999],
#            dtype='int64', length=1000000)

この index で結合 ( DataFrame.join ) した場合の処理時間は以下 (自身との結合のため意味のない処理だが)。

%timeit df.join(df, rsuffix='right_')
# 10 loops, best of 3: 70.6 ms per loop

このとき、 DataFrameに対して前処理 (スライス / サンプリングなど) を行うと、処理に応じて index も並びかわる。例えばサンプリング ( DataFrame.sample ) を行った場合、

df2 = df.sample(n=len(df))
df2.shape
# (1000000, 3)

df2.index
# Int64Index([372149, 955220, 320961, 244572, 254656,  74192, 143279, 246307,
#             764579,  96091, 
#             ...
#             827030, 681861, 492455, 894210, 758153, 327280, 245717, 952466,
#              26440, 532620],
#            dtype='int64', length=1000000)

このような DataFrame で結合を行うと、上と比較して処理時間が長くなっていることがわかる。

%timeit df2.join(df, rsuffix='right_')
# 1 loops, best of 3: 185 ms per loop

index の値がユニークかどうか / ソートされているかは以下のプロパティで確認できる。このプロパティによって一部の内部処理が分かれる。

# ユニークかどうか
df.index.is_unique
# True

# 昇順でソートされているか
df.index.is_monotonic_increasing
# True

df2.index.is_unique
# True

df2.index.is_monotonic_increasing
# False

まとめ

pandas にその本来のパフォーマンスを発揮させるためには、以下 3 点の処理を避けるとよい。

行に対するループ / DataFrame.apply は使わない
object 型は使わない
ユニークでない / ソートされていない index は使わない

上記を行った上で処理速度に不満がある場合は Cython や Numba で高速化する。もしくは、Dask を使って並列化を考える。

9/24 追記 Dask についてはこちらを。

sinhrks.hatenablog.com

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
出版社/メーカー: オライリージャパン
発売日: 2013/12/26
メディア: 大型本
この商品を含むブログ (10件) を見る

2015-06-18

Python pandas データのイテレーションと関数適用、pipe

pandas Python

pandas ではデータを列や表形式のデータ構造として扱うが、これらのデータから順番に値を取得 (イテレーション) して何か操作をしたい / また何らかの関数を適用したい、ということがよくある。このエントリでは以下の 3 つについて整理したい。

イテレーション
関数適用
pipe (0.16.2 で追加)

それぞれ、Series、DataFrame、GroupBy (DataFrame.groupbyしたデータ) で可能な操作が異なるため、順に記載する。

まずは必要なパッケージを import する。

import numpy as np
import pandas as pd

イテレーション

`Series`

Series は以下 2つのイテレーション用メソッドを持つ。各メソッドの挙動は以下のようになる。

__iter__: Series の値 ( values ) のみをイテレーション
iteritems: Series の index と values からなる tuple をイテレーション

図で表すとこんな感じ。矢印が処理の方向、枠内が 1 処理単位。

f:id:sinhrks:20150618213404p:plain

s = pd.Series([1, 2, 3], index=['a', 'b', 'c'])

for v in s:
    print(v)
# 1
# 2
# 3

for i, v in s.iteritems():
    print(i)
    print(v)
    print('')
# a
# 1
# 
# b
# 2
# 
# c
# 3

`DataFrame`

DataFrame は以下 3つのイテレーション用メソッドを持つ。同様に挙動を示す。

__iter__: DataFrame の列名 ( columns ) のみをイテレーション
iteritems: DataFrame の列名と列の値 ( Series ) からなる tuple をイテレーション
iterrows: DataFrame の行名と行の値 ( Series ) からなる tuple をイテレーション

f:id:sinhrks:20150618222231p:plain

df = pd.DataFrame({'A': [1, 2, 3], 'B': [4, 5, 6]}, index=['a', 'b', 'c'])
df
#    A  B
# a  1  4
# b  2  5
# c  3  6

for col in df:
    print(col)
# A
# B

for key, column in df.iteritems():
    print(key)
    print(column)
    print('')
# A
# a    1
# b    2
# c    3
# Name: A, dtype: int64
# 
# B
# a    4
# b    5
# c    6
# Name: B, dtype: int64

for key, row in df.iterrows():
    print(key)
    print(row)
    print('')
# a
# A    1
# B    4
# Name: a, dtype: int64
# 
# b
# A    2
# B    5
# Name: b, dtype: int64
# 
# c
# A    3
# B    6
# Name: c, dtype: int64

`GroupBy`

GroupBy は以下のイテレーション用メソッドを持つ。

__iter__: GroupBy のグループ名とグループ ( DataFrame もしくは Series ) からなる tuple をイテレーション

f:id:sinhrks:20150618213434p:plain

df = pd.DataFrame({'group': ['g1', 'g2', 'g1', 'g2'],
                   'A': [1, 2, 3, 4], 'B': [5, 6, 7, 8]}, 
                  columns=['group', 'A', 'B'])
df
#   group  A  B
# 0    g1  1  5
# 1    g2  2  6
# 2    g1  3  7
# 3    g2  4  8

grouped = df.groupby('group')

for name, group in grouped:
    print(name)
    print(group)
    print('')
# g1
#   group  A  B
# 0    g1  1  5
# 2    g1  3  7
# 
# g2
#   group  A  B
# 1    g2  2  6
# 3    g2  4  8

関数適用

`Series`

Series の各値に対して関数を適用する方法は以下の 2 つ。挙動はほぼ一緒だが、関数適用する場合は apply を使ったほうが意図が明確だと思う

Series.apply: Series の各値に対して関数を適用。
Series.map: Series の各値を、引数を用いてマッピングする。引数として、dict や Series も取れる。

f:id:sinhrks:20150618213448p:plain

s = pd.Series([1, 2, 3], index=['a', 'b', 'c'])

s.apply(lambda x: x * 2)
# a    2
# b    4
# c    6
# dtype: int64

# apply の引数には、Series の値そのものが渡っている
s.apply(type)
# a    <type 'numpy.int64'>
# b    <type 'numpy.int64'>
# c    <type 'numpy.int64'>
# dtype: object

# 関数が複数の返り値を返す場合、結果は tuple の Series になる
s.apply(lambda x: (x, x * 2))
# a    (1, 2)
# b    (2, 4)
# c    (3, 6)
# dtype: object

# 結果を DataFrame にしたい場合は、返り値を Series として返す
s.apply(lambda x: pd.Series([x, x * 2], index=['col1', 'col2']))
#    col1  col2
# a     1     2
# b     2     4
# c     3     6

# map の挙動は apply とほぼ同じ (map では結果を DataFrame にすることはできない)
s.map(lambda x: x * 2)
# a    2
# b    4
# c    6
# dtype: int64

s.map(type)
# a    <type 'numpy.int64'>
# b    <type 'numpy.int64'>
# c    <type 'numpy.int64'>
# dtype: object

# map は 関数以外に、 mapping 用の dict や Series を引数として取れる
s.map(pd.Series(['x', 'y', 'z'], index=[1, 2, 3]))
# a    x
# b    y
# c    z
# dtype: object

`DataFrame`

DataFrame に対して関数を適用する方法は以下の 2 つ。

DataFrame.apply: DataFrame の各列もしくは各行に対して関数を適用。行 / 列の指定は axis キーワードで行う。
DataFrame.applymap: DataFrame の各値に対して関数を適用。

f:id:sinhrks:20150618213500p:plain

df = pd.DataFrame({'A': [1, 2, 3], 'B': [4, 5, 6]}, index=['a', 'b', 'c'])
df
#    A  B
# a  1  4
# b  2  5
# c  3  6

# 各列に対して関数適用
df.apply(lambda x: np.sum(x))
A     6
B    15
dtype: int64

# 各行に対して関数適用
df.apply(lambda x: np.sum(x), axis=1)
a    5
b    7
c    9
dtype: int64

# 各値に対して関数適用
df.applymap(lambda x: x * 2)
#    A   B
# a  2   8
# b  4  10
# c  6  12

# apply で適用される関数には、各列もしくは各行が Series として渡される
df.apply(type)
# A    <class 'pandas.core.series.Series'>
# B    <class 'pandas.core.series.Series'>
# dtype: object

# applymap で適用される関数には、値そのものが引数として渡される
df.applymap(type)
#                       A                     B
# a  <type 'numpy.int64'>  <type 'numpy.int64'>
# b  <type 'numpy.int64'>  <type 'numpy.int64'>
# c  <type 'numpy.int64'>  <type 'numpy.int64'>

`GroupBy`

GroupBy については、GroupBy.apply で各グループに関数を適用できる。

f:id:sinhrks:20150618213517p:plain

df = pd.DataFrame({'group': ['g1', 'g2', 'g1', 'g2'],
                   'A': [1, 2, 3, 4], 'B': [5, 6, 7, 8]}, 
                  columns=['group', 'A', 'B'])
df
#   group  A  B
# 0    g1  1  5
# 1    g2  2  6
# 2    g1  3  7
# 3    g2  4  8

grouped = df.groupby('group')

grouped.apply(np.mean)
#        A  B
# group      
# g1     2  6
# g2     3  7

補足処理最適化のため、対象となるグループの数 == 関数適用の実行回数とはならないので注意。関数中で破壊的な処理を行うと意図しない結果になりうる。

# 適用される関数
def f(x):
    print('called')
    return x

# グループ数は 2
grouped.ngroups
# 2

# f の実行は 3 回
grouped.apply(f)
# called
# called
# called

`pipe`

先日リリースされた v0.16.2 にて pipe メソッドが追加された。これは R の {magrittr} というパッケージからインスパイアされたもので、データへの連続した操作をメソッドチェイン (複数のメソッドの連続した呼び出し) で記述することを可能にする。

Series.pipe、DataFrame.pipe それぞれ、x.pipe(f, *args, **kwargs) は f(x, *args, **kwargs) と同じ。つまり、データ全体に対する関数適用になる。

f:id:sinhrks:20150618213530p:plain

補足 GroupBy.pipe は v0.16.2 時点では存在しない。

# 渡される型は呼び出し元のインスタンス
s.pipe(type)
# pandas.core.series.Series

df.pipe(type)
# pandas.core.frame.DataFrame

np.random.seed(1)
df = pd.DataFrame(np.random.randn(10, 10))

# DataFrame を引数として heatmap を描く関数を定義
def heatmap(df):
    import matplotlib.pyplot as plt
    fig, ax = plt.subplots()
    return ax.pcolor(df.values, cmap='Greens')

# heatmap(df) と同じ。
df.pipe(heatmap)

f:id:sinhrks:20150618221027p:plain

まとめ

イテレーション、関数適用、pipe について整理した。特に関数適用はデータの前処理時に頻出するため、パターンを覚えておくと便利。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
出版社/メーカー: オライリージャパン
発売日: 2013/12/26
メディア: 大型本
この商品を含むブログ (9件) を見る

2015-06-14

Python pandas + folium + Jupyter でリーフレット / コロプレス図を描きたい

Python 可視化 pandas

引き続き、 R の可視化を Python に持ってくるシリーズ。R には以下のようなパッケージがあり、地図上へのリーフレット配置やコロプレス図の描画がカンタンにできる。それぞれの概要はリンク先を。

{leaflet}: リーフレット配置
- Leaflet for R - Introduction (英語)
{choroplethr}: コロプレス図の描画
- choroplethrで大阪市のコロプレス図を描く - Technically, technophobic.

これを Python でやりたい。調べてみると folium というパッケージが上記両方をサポートしているようなので使ってみる。

<a href="https://github.com/python-visualization/folium">python-visualization/folium</a>github.com

インストール

pip で。

pip install folium

準備

以降の操作は Jupyter Notebook から行う。まずはパッケージをロードする。

import numpy as np
import pandas as pd
import folium

folium はプロット結果を html としてエクスポートすることを想定して作成されているようだ。そのため、結果を Jupyter 上に埋め込みたい場合は以下のような関数を定義する必要がある。

from IPython.display import HTML

def inline_map(m):
    # 中間生成される json が必要なプロットがあるため、一度 html として書き出し
    m.create_map(path='tmp.html')
    iframe = '<iframe srcdoc=\"{srcdoc}\" style=\"width: 100%; height: 400px; border: none\"></iframe>'
    return HTML(iframe.format(srcdoc=m.HTML.replace('\"', '&quot;')))

リーフレット

手順は以下のようになる。

folium.Map で地図を描画する範囲を緯度経度/ズームにより指定
Map.simple_marker で緯度経度を指定してリーフレットを配置 (複数配置する場合は繰り返し)
Map.create_map で地図を含む html を生成 ( Jupyter 上に描画する場合は上で定義した関数 inline_map を呼ぶ )

m = folium.Map(location=[33.763, -84.392], zoom_start=17)
m.simple_marker([33.763006, -84.392912], popup='World of Coca-Cola')
inline_map(m)

f:id:sinhrks:20150614211105p:plain

補足 Jupyter 上ではスクロール / 拡大縮小できる。

既定では OpenStreetMap が利用されるが、Mapbox や maps.stamen を使うこともできる。また、リーフレットのマーカーとしては以下のものが利用できる。

Simple Markers: シンプルなマーカー (上のもの)
Circle Markers: 円形のマーカー
Polygon Markers: 多角形のマーカー
Lat/Lng Popups: 緯度経度を表示するマーカー
Click-for-Marker: クリックで配置可能なマーカー
Vincent Popups: Vincent によるプロットを埋め込めるマーカー

それぞれの描画サンプルは README で確認することができる。

コロプレス図

コロプレス図を描くには以下２つのデータソースが必要である。

GeoJSON もしくは TopoJSON 形式のファイル
コロプレス図を色分けするための値を含む pandas の DataFrame

サンプルとして国別のマクドナルドの店舗数をプロットする。

`GeoJSON` ファイルの準備

国別にプロットするため、国別の GeoJSON ファイルがほしい。以下リポジトリの countries.geo.json をローカルに保存して使う。

johan/world.geo.json · GitHub

中身は以下のように ISO 3166-1 alpha-3 の国コードを id としたデータとなっている。

{"type":"FeatureCollection",
 "features":[{"type":"Feature","id":"AFG",
              "properties":{"name":"Afghanistan"},
              "geometry":{"type":"Polygon","coordinates":[[[61.210817,35.650072],[62.230651,35.270664],...

`DataFrame` の準備

DataFrame は上で準備した GeoJSON と紐づけるためのキー (一般には GeoJSON の id ) と値の 2 列をもつ必要がある。国別のマクドナルドの店舗数データを pd.read_html で読み込む。

url = "https://en.wikipedia.org/wiki/List_of_countries_with_McDonald%27s_restaurants"
df = pd.read_html(url, header=0, index_col=0)[0]
# 列名を変更
df.columns = ['Country', 'Date', 'First outlet location', 'Number', 'Source', 'Note', 'CEO']
df[['Country', 'Number']].head()
#                Country Number
# #                            
# 1        United States  14267
# 2               Canada   1427
# 3          Puerto Rico    108
# 4  U.S. Virgin Islands      6
# 5           Costa Rica     54

こちらのデータには国名が入っているため、GeoJSON と紐づけるためには ISO 国コードに変換する必要がある。変換には pycountry を使う。インストールしていない方は pip で。

import pycountry
pycountry.countries.get(name='Japan').alpha3
# 'JPN'

# データ中の国名が pycountry のものと違う場合の mapper
countries = {'United Kingdom': 'United Kingdom',
             'Russia': 'Russian Federation',
             'South Korea': 'Korea, Republic of',
             'Taiwan': 'Taiwan, Province of China',
             'Vietnam': 'Viet Nam', }
def f(x):
    for k, v in pd.compat.iteritems(countries):
        if k in x:
            x = v
    try:
        return pycountry.countries.get(name=x).alpha3
    except KeyError:
        return np.nan

# 国コードへの変換
df.loc[:, 'Code'] = df['Country'].apply(f)
# 国コードに変換できなかったデータは除外
df = df.dropna(subset=['Code'])

# 欠損値を 0 でパディング
df.loc[:, 'Number'] = df['Number'].fillna('0')
# 数値に変換できない文字列があるため、余計な文字を削除
df.loc[:, 'Number'] = df['Number'].str.replace('[+,]', '')
# 数値 (float) 型に変換
df.loc[:, 'Number'] = df['Number'].astype(float)
df[['Code', 'Country', 'Number']].sort('Number', ascending=False).head()
#    Code        Country  Number
# #                             
# 1   USA  United States   14267
# 8   JPN          Japan    3164
# 49  CHN          China    2000
# 11  DEU        Germany    1468
# 2   CAN         Canada    1427

これで、国コード / プロットする値をもつ DataFrame が準備できた。アメリカすごいな、、、。

コロプレス図の描画

Map.geo_json でコロプレス図を描画するレイヤーを Map に追加できる。ここで使っている引数の意味は以下。

geo_path: GeoJSON ファイルのパス
data: 色分けのための値をもつ DataFrame
columns: data 中 GeoJSON と紐づけるキー / 値を含む列名
key_on: GeoJSON 側で紐付けに使うキー
threshold_scale: 色分けをする際の閾値
fill_color 色分けに使う color-brewer の名前
reset: 既存のレイヤがある場合に削除する

m = folium.Map(location=[10, 35], zoom_start=1.5)

geojson = r'countries.geo.json'
m.geo_json(geo_path=geojson, data=df,
           columns=['Code', 'Number'],
           key_on='feature.id',
           threshold_scale=[1, 100, 500, 1000, 2000, 4000],
           fill_color='BuPu', reset=True)
inline_map(m)

f:id:sinhrks:20150614211127p:plain

まとめ

folium を使えばリーフレット / コロプレス図の描画がカンタンにできる。Jupyter 上で Javascript を使用してのデータ可視化、結構使えるのでは。

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
出版社/メーカー: オライリージャパン
発売日: 2013/12/26
メディア: 大型本
この商品を含むブログ (9件) を見る

2015-06-13

Python pandas のデータを Highcharts/Highstock + Jupyter でプロットしたい

pandas Python 可視化

R を使っている方はご存知だと思うが、R には {htmlwidgets} というパッケージがあり、R 上のデータを任意の Javascript ライブラリを使ってプロットすることが比較的カンタンにできる。{htmlwidgets} って何？という方にはこちらの説明がわかりやすい。

RPubs - htmlwidgetsでJavascriptの可視化をRに

同じことを Python + pandas を使ってやりたい。サンプルとして利用する Javascript ライブラリは上の資料と同じく Highcharts、Highstock にする。

<a href="http://www.highcharts.com/">Interactive JavaScript charts for your webpage | Highcharts</a>www.highcharts.com

補足 pandas-highcharts という Python パッケージもあるが、このエントリでは任意の Javascript ライブラリで使えるであろう方法を記載する。

`Highcharts` でのプロット

以降の操作は Jupyter Notebook 上で行う。まずは必要パッケージをロードする。

import numpy as np
import pandas as pd
from IPython.display import HTML

続けて、Highcharts、Highstock を読み込む。これは %%html magic を使うのが楽。

%%html
    <script src="http://code.highcharts.com/highcharts.js"></script>
    <script src="http://code.highcharts.com/stock/highstock.js"></script>
    <script src="http://code.highcharts.com/modules/exporting.js"></script>

Highcharts でプロットする準備ができたため、%%html magic を使ってサンプル Your first chart に記載のサンプルをプロットしてみる。これでプロット時のスクリプト / データ構造が確認できる。

%%html
    <div id="container" style="width:100%; height:400px;"></div>
    <script>
        plot = function () { 
            $('#container').highcharts({
                chart: {
                    type: 'bar'
                },
                title: {
                    text: 'Fruit Consumption'
                },
                xAxis: {
                    categories: ['Apples', 'Bananas', 'Oranges']
                },
                yAxis: {
                    title: {
                        text: 'Fruit eaten'
                    }
                },
                series: [{
                    name: 'Jane',
                    data: [1, 0, 4]
                }, {
                    name: 'John',
                    data: [5, 7, 3]
                }]
            });
        };
        plot();
    </script>

f:id:sinhrks:20150613152238p:plain

補足 Jupyter から実行すればアニメーションする。

pandas のデータをプロットしたい場合は、上のスクリプトと同じように .highcharts() の引数にあわせた形式でデータを渡し、HTML としてレンダリングしてやればうまくいきそうだ。pandas で元データとなる DataFrame を定義して、上のスクリプトの形式に変換していく。

df = pd.DataFrame({'Jane': [1, 0, 4], 'John': [5, 7, 3]},
                  index=['Apples', 'Bananas', 'Oranges'])
df
#          Jane  John
# Apples      1     5
# Bananas     0     7
# Oranges     4     3

スクリプトは文字列結合で作ってもよいが、引数となる json 形式に対応する辞書型のデータをつくってからjson.dumps したほうが簡単だろう。DataFrame.to_json でデータを直接変換できると楽なのだが、フォーマットが違うため無理そうだ。

# NG!
df.to_json()
# '{"Jane":{"Apples":1,"Bananas":0,"Oranges":4},
#   "John":{"Apples":5,"Bananas":7,"Oranges":3}}'

そのため、個々の要素ごとに変換を考えていく。うまいこと辞書ができたら json.dumps する。

[{'name': c, 'data': col.tolist()} for c, col in df.iteritems()]
# [{'data': [1, 0, 4], 'name': 'Jane'}, {'data': [5, 7, 3], 'name': 'John'}]

chartdict = {'chart': {'type': 'bar'},
             'title': {'text': 'Fruit Consumption'},
             'xAxis': {'categories': df.index.tolist()}, 
             'yAxis': {'title': {'text': 'Fruit eaten'}},
             'series': [{'name': c, 'data': col.tolist()} for c, col in df.iteritems()]
            }

import json
json.dumps(chartdict)
# '{"series": [{"data": [1, 0, 4], "name": "Jane"}, {"data": [5, 7, 3], "name": "John"}],
#   "yAxis": {"title": {"text": "Fruit eaten"}}, "chart": {"type": "bar"},
#   "xAxis": {"categories": ["Apples", "Bananas", "Oranges"]},
#   "title": {"text": "Fruit Consumption"}}'

あとは必要な HTML / Javascript のテンプレートを文字列として作って format すればよい。

template = """
<script src="http://code.highcharts.com/highcharts.js"></script>
<script src="http://code.highcharts.com/modules/exporting.js"></script>
<div id="{chart}" style="width:100%; height:400px;"></div>
<script type="text/javascript">
plot = function () {{
    $("#{chart}").highcharts({data});
}};
plot();
</script>
"""
HTML(template.format(chart='container2', data=json.dumps(chartdict)))
# 略

`Highstock` でのプロット

同様に、Highstock へプロットすることもできる。サンプル Single line series を pandas のデータを使ってプロットしてみる。

補足株価の取得には以下のパッケージを使う。

Python pandas 日本語環境向けのちょっとしたパッケージ作った - StatsFragments

import japandas as jpd
toyota = jpd.DataReader(7203, 'yahoojp', start='2015-01-01')
toyota.head()
#               始値    高値    安値    終値       出来高  調整後終値*
# 日付                                                  
# 2015-01-05  7565  7575  7416  7507   9515300    7507
# 2015-01-06  7322  7391  7300  7300  12387900    7300
# 2015-01-07  7256  7485  7255  7407  11465400    7407
# 2015-01-08  7500  7556  7495  7554  10054500    7554
# 2015-01-09  7630  7666  7561  7609  10425400    7609

Highstock に渡すデータの形式は以下のサイトがわかりやすい。引数としては [UNIX時間(ミリ秒), 始値, 高値, 安値, 終値] を入れ子のリストにして渡せばよいようだ。

highstockの使い方初心者向け

時刻は UNIX時間で渡す必要があるので、少し操作が必要だ。上で取得した DataFrame は日時型の Index である DatetimeIndex を持っている。DatetimeIndex はUNIXエポックを基準とした現在時刻をナノ秒で保存しているため、int 型に変換して 1000000 で割ればUNIX時間(ミリ秒)となる。したがって、Highstock へ渡すデータは以下のようにして作れる。

toyota.index.astype(int) / 1000000
# array([1420416000000, 1420502400000, 1420588800000, 1420675200000,
#        1420761600000, 1421107200000, 1421193600000, 1421280000000,
#        ....
#        1433721600000, 1433808000000, 1433894400000, 1433980800000,
#        1434067200000])

toyota['time'] = toyota.index.astype(int) / 1000000
toyota[['time', u'始値', u'高値', u'安値', u'終値']].values.tolist()
# [[1420416000000, 7565, 7575, 7416, 7507],
#  [1420502400000, 7322, 7391, 7300, 7300],
#  ...
#  [1433980800000, 8250, 8326, 8241, 8322],
#  [1434067200000, 8387, 8394, 8329, 8394]]

描画する。アニメーションも含め、うまく動いているようだ。

chartdict = {'rangeSelector': {'selected': 1},
             'title': {'text' : 'Stock Price'},
             'series': [{'name' : u'トヨタ',
                         'data': toyota[['time', u'始値', u'高値', u'安値', u'終値']].values.tolist(),
                         'tooltip': {'valueDecimals': 2}}]}

template = """
<div id="{chart}" style="width:100%; height:400px;"></div>
<script type="text/javascript">
plot = function () {{
    $('#{chart}').highcharts('StockChart', {data});
}};
plot();
</script>
HTML(template.format(chart='container3', data=json.dumps(chartdict)))

f:id:sinhrks:20150613145918p:plain

まとめ

pandas のデータを任意の Javascript ライブラリでプロットする際には、

当該の Javascript ライブラリが利用するデータ構造を確認する
そのデータ構造にあうように pandas のデータを変換し、辞書型を作る
json.dumps でスクリプトに渡す

とやればだいたいできると思う。

補足ここで今話題の spyre 上でもプロットしたいな？と思ったのだが、spyre で普通に HTML としてレンダリングするとうまく動かない。spyre が内部で使っている cherrypy も DOM を書き換えているのが原因かと思っているのだが、自分にはそのあたりの知識がまったくないのでよくわからない。できた方がいれば教えてください。

<a href="http://sinhrks.hatenablog.com/entry/2015/06/13/003255">Python spyre によるデータ分析結果のWebアプリ化 - StatsFragments</a>sinhrks.hatenablog.com

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
出版社/メーカー: オライリージャパン
発売日: 2013/12/26
メディア: 大型本
この商品を含むブログ (9件) を見る

2015-06-13

Python spyre によるデータ分析結果のWebアプリ化

Python 可視化 pandas

R を使っている方はご存知だと思うが、R には {Shiny} というパッケージがあり、データ分析の結果をインタラクティブな Web アプリとして共有することができる。{Shiny} って何？という方にはこちらの説明がわかりやすい。

<a href="http://qiita.com/hoxo_m/items/c8365117f3444fb51df4">RStudio Shiny チュートリアルレッスン1 ようこそ Shiny へ - Qiita</a>qiita.com

Python でも {Shiny} のようなお手軽可視化フレームワークがあるといいよね、とたびたび言われていたのだが、spyre というなんかそれっぽいパッケージがあったので触ってみたい。

<a href="https://github.com/adamhajari/spyre">adamhajari/spyre</a>github.com

インストール

pip で。

pip install dataspyre

使い方

現時点で ドキュメンテーションはない ので、README と examples ディレクトリを見る。サンプルとして株価を取得してプロットするWebアプリを作ってみたい。spyre で Webアプリを作る手順は以下の3つ。

spyre.server.App を継承したクラスを作る。
描画をコントロール/指示するクラス変数を指定する。
描画を行うメソッド getData, getPlot を書く。メソッド名は描画内容 (output_type) ごとに固定。

最初は examples のものを書き換えながら作るのが楽。各クラス変数/メソッドは相互に関連するため、プログラム全体を示した上で必要と思われる箇所にコメントを入れた。

補足株価の取得には以下のパッケージを使う。

Python pandas 日本語環境向けのちょっとしたパッケージ作った - StatsFragments

#!/usr/bin/env python
# -*- coding: utf-8 -*-
from spyre import server

import pandas as pd
pd.options.display.mpl_style = 'default'

# あらかじめデータを取得しておく
# 終値のみを取得し、一つのDataFrameに結合
import japandas as jpd
toyota = jpd.DataReader(7203, 'yahoojp', start='2015-01-01')[[u'終値']]
toyota.columns = [u'トヨタ']
honda = jpd.DataReader(7267, 'yahoojp', start='2015-01-01')[[u'終値']]
honda.columns = [u'ホンダ']
df = toyota.join(honda)


# spyre.server.App を継承したクラスを作る
class StockExample(server.App):
    title = u"株価のプロット"

    # 左側のペインに表示する UI 要素を辞書のリストで指定
    # ここではドロップダウン一つだけを表示
    inputs = [{"input_type":'dropdown', 
               # ドロップダウン自体の表示ラベル
               "label": 'Frequency',
               # ドロップダウンの選択項目を指定
               # label はドロップダウン項目の表示ラベル
               # value は各項目が選択された時にプログラム中で利用される値
               "options" : [ {"label": "月次", "value":"M"},
                             {"label": "週次", "value":"W"},
                             {"label": "日次", "value":"B"}],
               # 各 UI 要素の入力は各描画メソッド (getData, getPlot) に
               # 辞書型の引数 params として渡される
               # その辞書から値を取り出す際のキー名
               "variable_name": 'freq',
               "action_id": "update_data" }]

    # 画面を更新する設定
    controls = [{"control_type" : "hidden",
                 "label" : "update",
                 "control_id" : "update_data"}]

    # 描画するタブの表示ラベルを文字列のリストで指定
    tabs = [u"トヨタ", u"ホンダ", u"データ"]

    # tabs で指定したそれぞれのタブに描画する内容を辞書のリストで指定
    outputs = [{"output_type" : "plot", # matplotlib のプロットを描画する
                "output_id" : "toyota", # 描画するタブに固有の id
                "control_id" : "update_data",
                "tab" : u"トヨタ",       # 描画するタブの表示ラベル (tabs に含まれるもの)
                "on_page_load" : True },

               {"output_type" : "plot",
                "output_id" : "honda",
                "control_id" : "update_data",
                "tab" : u"ホンダ",
                "on_page_load" : True },

               {"output_type" : "table", # DataFrameを描画する
                "output_id" : "table_id",
                "control_id" : "update_data",
                "tab" : u"データ",
                "on_page_load" : True }]

    def getData(self, params):
        """
        output_type="table" を指定したタブを描画する際に呼ばれるメソッド
        DataFrameを返すこと

        params は UI 要素の入力 + いくつかのメタデータを含む辞書
        UI 要素の入力は inputs で指定した variable_name をキーとして行う
        """
        # ドロップダウンの値を取得
        # 値にはユーザの選択によって、options -> value で指定された M, W, B いずれかが入る
        freq = params['freq']

        # freq でグループ化し平均をとる
        tmp = df.groupby(pd.TimeGrouper(freq)).mean()
        return tmp

    def getPlot(self, params):
        """
        output_type="plot" を指定したタブを描画する際に呼ばれるメソッド
        matplotlib.Figureを返すこと
        """
        tmp = self.getData(params)

        # 同じ output_type で複数のタブを描画したい場合は、 params に含まれる
        # output_id で分岐させる
        # output_id は タブの表示ラベルではなく、outputs 中で指定した output_id
        if params['output_id'] == 'toyota':
            ax = tmp[[u'トヨタ']].plot(legend=False)
            return ax.get_figure()
        elif params['output_id'] == 'honda':
            ax = tmp[[u'ホンダ']].plot(legend=False)
            return ax.get_figure()
        else:
            raise ValueError


app = StockExample()
# port 9093 で Webサーバ + アプリを起動
app.launch(port=9093)

実行後、ブラウザで http://127.0.0.1:9093 を開くと以下のような画面が表示される。inputs で指定した項目が左側のメニューに、 tabs で指定した項目がタブとして表示されている。

f:id:sinhrks:20150613001354p:plain

ドロップダウンやタブの選択を切り替えると、動的にグラフやデータが更新されることがわかる。

f:id:sinhrks:20150613001400p:plain

DataFrame の描画はちょっとおかしく、index である日時自体が描画されずにその名前 ("日付") だけが表示されている。これは後日 issue あげて直そうかと思う。

f:id:sinhrks:20150613001406p:plain

06/13追記 index を描画しないのは spyre の仕様で、そのとき名前が残ってしまうのは pandas のバグです。v0.17.0で修正予定。

BUG: df.to_html(index=False) renders index.name by sinhrks · Pull Request #10346 · pydata/pandas · GitHub

まとめ

お手軽可視化フレームワーク spyre を試してみた。現行の {Shiny} と比べると機能/UI とも十分ではないが、いくつかのデータ / プロットをインタラクティブに共有する用途であれば十分使えそうだ。

使いごこちの印象として、はじめて {Shiny} を触ったあの日の気持ちを思い出したような気がする、、。

2015-05-18

pandas 日時まわりのリサンプリング/オフセット処理

pandas Python

こちらの続き。

Python pandas で日時関連のデータ操作をカンタンに - StatsFragments

今回のサンプルデータには自分の歩数のデータを使いたい。インスパイヤ元は以下のサイトだ。

<a href="http://d.hatena.ne.jp/dichika/20150208/p1">下半身のデータをオープンにした - 盆栽日記</a>d.hatena.ne.jp

データの読み込み

歩数データは iPhone の Health アプリから Export できる。形式は XML なので、そのままでは pandas で読み込めない。一度 XML から必要な属性を辞書のリストとして取り出した後、pandas に読み込ませる。

import pandas as pd
from xml.etree import ElementTree
tree = ElementTree.parse('export.xml')
root = tree.getroot()
# 属性の辞書のリストを作る
data = [e.attrib for e in root.findall('.//Record')]
data[0]
# {'endDate': '20150511220900+0900',
#  'recordCount': '104',
#  'source': u'xxx',
#  'startDate': '20150511210900+0900',
#  'type': 'HKQuantityTypeIdentifierStepCount',
#  'unit': 'count',
#  'value': '531'}

df = pd.DataFrame(data)
# 必要な列のみにフィルタ
df = df[['startDate', 'type', 'unit', 'value']]
df.head()
#              startDate                               type   unit value
# 0  20150511210900+0900  HKQuantityTypeIdentifierStepCount  count   531
# 1  20150512000900+0900  HKQuantityTypeIdentifierStepCount  count     6
# 2  20150512060900+0900  HKQuantityTypeIdentifierStepCount  count   629
# 3  20150512070900+0900  HKQuantityTypeIdentifierStepCount  count   312
# 4  20150512080900+0900  HKQuantityTypeIdentifierStepCount  count  1483

Series.value_counts() を使うと、ある列で出現する要素の個数を集計できる。"type" 列を集計すると、Export したデータには歩数以外のデータも含まれていることがわかる。歩数以外は不要なので削除する。

df['type'].value_counts()
# HKQuantityTypeIdentifierStepCount                 375
# HKQuantityTypeIdentifierDistanceWalkingRunning    374
# HKQuantityTypeIdentifierFlightsClimbed            216
# dtype: int64

# 歩数のみにフィルタ
df = df[df['type'] == 'HKQuantityTypeIdentifierStepCount']

また、XML で読み込んだデータはすべて文字列型になっているため、適宜型変換を行う。

# 歩数データは数値に
df['value'] = df['value'].astype(float)

# 日付は index として設定し、日付型 (DatetimeIndex) に変換
df = df.set_index('startDate')

# この時点で index は文字列型になっている
df.index 
# Index([u'20150511210900+0900', u'20150512000900+0900', u'20150512060900+0900',
#        u'20150512070900+0900', u'20150512080900+0900', u'20150512090900+0900',
#        ...
#        u'20150518165000+0900', u'20150518165900+0900', u'20150518171600+0900',
#        u'20150518181200+0900'],
#       dtype='object', name=u'startDate', length=965)

# to_datetime で日時型に変換 / タイムゾーンを表すオフセットを適宜設定
pd.to_datetime(df.index, utc=True).tz_convert('Asia/Tokyo')
# DatetimeIndex(['2015-05-11 21:09:00+09:00', '2015-05-12 00:09:00+09:00',
#                '2015-05-12 06:09:00+09:00', '2015-05-12 07:09:00+09:00',
#                ...
#                '2015-05-18 16:50:00+09:00', '2015-05-18 16:59:00+09:00',
#                '2015-05-18 17:16:00+09:00', '2015-05-18 18:12:00+09:00'],
#               dtype='datetime64[ns]', length=965, freq=None, tz='Asia/Tokyo')

# index を上書き
df.index = pd.to_datetime(df.index, utc=True).tz_convert('Asia/Tokyo')

# タイムゾーンは使わないので削除
df.index = df.index.tz_localize(None)

# 日付の昇順にソート
df = df.sort_index()

df.head()
#                                                   type   unit     value
# 2014-11-25 21:09:00  HKQuantityTypeIdentifierStepCount  count   1396.00
# 2014-11-26 21:09:00  HKQuantityTypeIdentifierStepCount  count   7020.37
# 2014-11-27 21:09:00  HKQuantityTypeIdentifierStepCount  count   6413.63
# 2014-11-28 21:09:00  HKQuantityTypeIdentifierStepCount  count   8396.00
# 2014-11-29 21:09:00  HKQuantityTypeIdentifierStepCount  count  12411.00

今の機種に買い替えたのが昨年11末なので、買い替え以降すべてのデータが入っているようだ。

部分文字列によるスライシング

pandas での一般的なデータ選択については以下の記事にまとめた。

Python pandas データ選択処理をちょっと詳しく <前編> - StatsFragments

加えて、データが日付型の Index ( DatetimeIndex ) を持つとき、日付-like な文字列で __getitem__ すると、その期間にあてはまるデータを抽出してくれる。例えば 2015-05-17 日分のデータを抽出したければ以下のようにする。この方法を使うと好きな期間のデータを簡単に確認することができる。詳細は公式ドキュメントを参照。

df['2015-05-17']
#                                                   type   unit     value
# 2015-05-17 08:09:00  HKQuantityTypeIdentifierStepCount  count   18.0000
# 2015-05-17 11:09:00  HKQuantityTypeIdentifierStepCount  count  180.0000
# 2015-05-17 12:09:00  HKQuantityTypeIdentifierStepCount  count  934.0000
# 2015-05-17 13:09:00  HKQuantityTypeIdentifierStepCount  count  469.0000
# ...                                                ...    ...       ...
# 2015-05-17 21:16:00  HKQuantityTypeIdentifierStepCount  count   88.5087
# 2015-05-17 21:17:00  HKQuantityTypeIdentifierStepCount  count   23.1214
# 2015-05-17 21:18:00  HKQuantityTypeIdentifierStepCount  count   99.3283
# 2015-05-17 21:19:00  HKQuantityTypeIdentifierStepCount  count   41.6284
# 
# [14 rows x 3 columns]

リサンプリング

データは日時で処理したいので、まずは 1日ごとの合計を算出したい。こういうときは DataFrame.resample。リサンプリングする期間や集約関数は引数として指定できる。

s = df.resample('D', how='sum')
#                    value
# 2014-11-25   1396.000000
# 2014-11-26   7020.370000
# 2014-11-27   6413.630000
# 2014-11-28   8396.000000
# ...                  ...
# 2015-05-15   7561.998000
# 2015-05-16  10615.000000
# 2015-05-17   9208.000000
# 2015-05-18   8085.994898
# 
# [175 rows x 1 columns]

日次に集約した結果をプロットしてみる。

s.plot()

f:id:sinhrks:20150518222638p:plain

あまり意識はしていないのだがそこそこ歩いているようだ。5月初旬はちょっと出かけていたため歩数に異常値が出ている。

日付の標準化

日時データの適当な期間での集約は DataFrame.resample でできた。が、時には日時の補正のみを行い、集約はしたくない場合がある。

DatetimeIndex を日付ごとにまとめるのに一番簡単なのは .normalize。

pd.Timestamp('2015-05-01 23:59').normalize()
# Timestamp('2015-05-01 00:00:00')
pd.Timestamp('2015-05-02 00:00').normalize()
# Timestamp('2015-05-02 00:00:00')

他の日時関連のメソッドと同じく、Timestamp, DatetimeIndex 両方で使える。

df.tail().index
# DatetimeIndex(['2015-05-18 18:55:00', '2015-05-18 18:56:00',
#                '2015-05-18 18:57:00', '2015-05-18 18:58:00',
#                '2015-05-18 18:59:00'],
#               dtype='datetime64[ns]', freq=None, tz=None)

df.tail().index.normalize()
# DatetimeIndex(['2015-05-18', '2015-05-18',
#                '2015-05-18', '2015-05-18',
#                '2015-05-18'],
#               dtype='datetime64[ns]', freq=None, tz=None)

当然、集約すれば日付での .resample と同じ結果になる。

s = df.groupby(df.index.normalize())[['value']].sum()
s 
# 略

オフセット

また、より柔軟に日時補正を行うためにオフセットが提供されている ( 前回記事 )。オフセットを利用してデータを補正/集約する例を書きたい。

オフセットを使わずに .resample で月次平均をとると以下のような結果になる。

s.resample('M', how='mean')
#                    value
# 2014-11-30   7012.500000
# 2014-12-31   8435.741935
# 2015-01-31   9134.032258
# 2015-02-28   9323.821429
# 2015-03-31   9326.356129
# 2015-04-30   9938.533333
# 2015-05-31  10973.055439

オフセットを使って同じ処理をするには、日時を月末に補正するオフセット MonthEnd を使う。オフセットの一覧は公式ドキュメントを参照。

m = offsets.MonthEnd()
m
# <MonthEnd>

オフセットにはいくつか共通のメソッドがあるため、順に記載する。

まず、ある日付が当該のオフセット上に存在するか調べるには .onOffset。MonthEnd.onOffset の場合は、引数が月末の日付であるとき True になる。

m.onOffset(pd.Timestamp('2015-04-29'))
# False
m.onOffset(pd.Timestamp('2015-04-30'))
# True
m.onOffset(pd.Timestamp('2015-05-01'))
# False

オフセットを加減算することにより、日時を次の/前のオフセットへと移動できる。また、オフセットの加算と同一の処理として　.apply がある。

# 日付を常に 次のオフセットへ移動させる。
pd.Timestamp('2015-04-29') + m
# Timestamp('2015-04-30 00:00:00')
pd.Timestamp('2015-04-30') + m
# Timestamp('2015-05-31 00:00:00')
pd.Timestamp('2015-05-01') + m
# Timestamp('2015-05-31 00:00:00')

# 日付を常に 前のオフセットへ移動させる。
pd.Timestamp('2015-04-29') - m
# Timestamp('2015-03-31 00:00:00')
pd.Timestamp('2015-04-30') - m
# Timestamp('2015-03-31 00:00:00')
pd.Timestamp('2015-05-01') - m
# Timestamp('2015-04-30 00:00:00')

# 日付を常に 次のオフセットへ移動させる。
m.apply(pd.Timestamp('2015-04-29'))
# Timestamp('2015-04-30 00:00:00')
m.apply(pd.Timestamp('2015-04-30'))
# Timestamp('2015-05-31 00:00:00')
m.apply(pd.Timestamp('2015-05-01'))
# Timestamp('2015-05-31 00:00:00')

オフセットに乗っている日時は動かしたくなければ .rollforward / .rollback というメソッドを使う。

# 日付がオフセットに乗っていない場合 次のオフセットへ移動させる。
m.rollforward(pd.Timestamp('2015-04-29'))
# Timestamp('2015-04-30 00:00:00')
m.rollforward(pd.Timestamp('2015-04-30'))
# Timestamp('2015-04-30 00:00:00')
m.rollforward(pd.Timestamp('2015-05-01'))
# Timestamp('2015-05-31 00:00:00')

# 日付がオフセットに乗っていない場合 前のオフセットへ移動させる。
m.rollback(pd.Timestamp('2015-04-29'))
# Timestamp('2015-03-31 00:00:00')
m.rollback(pd.Timestamp('2015-04-30'))
# Timestamp('2015-04-30 00:00:00')
m.rollback(pd.Timestamp('2015-05-01'))
# Timestamp('2015-04-30 00:00:00')

よって、月次平均の算出をオフセットのメソッドを使って行う場合は以下のようになる。

# 例示のため一部データをスライス
s['2014-11']
#                value
# 2014-11-25   1396.00
# 2014-11-26   7020.37
# 2014-11-27   6413.63
# 2014-11-28   8396.00
# 2014-11-29  12411.00
# 2014-11-30   6438.00

# 日時を月末に補正
s['2014-11'].index.map(m.rollforward)
# array([Timestamp('2014-11-30 00:00:00'), Timestamp('2014-11-30 00:00:00'),
#        Timestamp('2014-11-30 00:00:00'), Timestamp('2014-11-30 00:00:00'),
#        Timestamp('2014-11-30 00:00:00'), Timestamp('2014-11-30 00:00:00')], dtype=object)

s.groupby(s.index.map(m.rollforward)).mean()
#                    value
# 2014-11-30   7012.500000
# 2014-12-31   8435.741935
# 2015-01-31   9134.032258
# 2015-02-28   9323.821429
# 2015-03-31   9326.356129
# 2015-04-30   9938.533333
# 2015-05-31  10973.055439

補足オフセットは .resample の引数として渡すこともでき、同じ結果になる。

s.resample(m, how='mean')
# 略

オフセットを活用した集計

つづけて、カレンダー上の平日 / 休日での差異をみたい。アクセサ ( こちらの記事参照 ) を使って曜日で集計してから処理してもよいが、 BusinessDay オフセットを使えば簡単。BusinessDay.onOffset を使えば平日が True / 休日を False として集計できる。以下の結果をみると、休日の方が歩いている傾向があるようだ。

bday = offsets.BusinessDay()
s.groupby(bday.onOffset).mean()
#               value
# False  10860.216800
# True    8716.657583

BusinessDay オフセットは土日のみを休日として扱うが、任意の祝日も含めて休日として扱いたい場合は CustomBusinessDay オフセットが使える。拙作のパッケージで日本の祝日のカレンダーを定義しているので、それを使って、

import japandas as jpd
calendar = jpd.JapaneseHolidayCalendar()
cday = pd.offsets.CDay(calendar=calendar)

s.groupby(cday.onOffset).mean()
#               value
# False  10799.033448
# True    8600.419640

これまでの内容を使って、月ごとに平日 / 土日 / 祝日の歩数平均をクロス集計したい。クロス集計を行うには pd.pivot_table 。

avg = pd.pivot_table(s, index=s.index.map(m.rollforward),
                     columns=[s.index.map(bday.onOffset), s.index.map(cday.onOffset)],
                     values='value', aggfunc='mean')
avg
#                    False    True              
#                    False    False        True 
# 2014-11-30   9424.500000      NaN  5806.500000
# 2014-12-31  10382.275000  11704.0  7579.354545
# 2015-01-31   9534.193333  10162.5  8851.113000
# 2015-02-28  10943.125000   9446.0  8635.578947
# 2015-03-31  11148.877778      NaN  8580.779091
# 2015-04-30  11349.625000   7110.0  9535.666667
# 2015-05-31  12769.000000  11582.7  9572.544211

左から、

土日: BusinessDay CustomBusinessDay ともに False
祝日: BusinessDay は True, CustomBusinessDay は False
平日: BusinessDay CustomBusinessDay ともに True

となる。ということで上の結果からも 5月の休日は結構歩いたなってことがわかった。

まとめ

日付を index として持つデータに対するリサンプリング/オフセット処理をまとめた。これらは以下のような使い分けをすればよいと思う。

比較的単純な集約はリサンプリング
リサンプリングではできないより細かい補正やクロス集計をしたい場合はオフセット

Pythonによるデータ分析入門 ―NumPy、pandasを使ったデータ処理

作者: Wes McKinney,小林儀匡,鈴木宏尚,瀬戸山雅人,滝口開資,野上大介
出版社/メーカー: オライリージャパン
発売日: 2013/12/26
メディア: 大型本
この商品を含むブログ (9件) を見る

はじめに

インストール

準備

pandas から Dask への変換

内部処理

実データでの利用例

まとめ

はじめに

演算子のオーバーライド

右辺側クラスによるオーバーライド

NumPy に対するオーバーライド

まとめ

はじめに

xray とは

インストール

データの準備

次元の追加

まとめ

geopandas のインストール

地理情報の読み込み

地理情報の処理

地理情報のプロット

Geocoding の利用

Geocoding を利用したデータの追加

まとめ

環境構築

サンプルデータの作成

1. 行に対するループ / DataFrame.apply は 使わない

2. object 型は使わない

3. ユニークでない / ソートされていない index は使わない

まとめ

Series

DataFrame

GroupBy

関数適用

Series

DataFrame

GroupBy

pipe

まとめ

インストール

準備

リーフレット

コロプレス図

GeoJSON ファイルの準備

DataFrame の準備

コロプレス図の描画

まとめ

Highcharts でのプロット

Highstock でのプロット

まとめ

インストール

使い方

まとめ

データの読み込み

部分文字列によるスライシング

リサンプリング

日付の標準化

オフセット

オフセットを活用した集計

まとめ

`pandas` から `Dask` への変換

`xray` とは

`geopandas` のインストール

1. 行に対するループ / `DataFrame.apply` は使わない

2. `object` 型は使わない

3. ユニークでない / ソートされていない `index` は使わない

`Series`

`DataFrame`

`GroupBy`

`Series`

`DataFrame`

`GroupBy`

`pipe`

`GeoJSON` ファイルの準備

`DataFrame` の準備

`Highcharts` でのプロット

`Highstock` でのプロット