MySQLにあるセンシングデータをヒートマップで可視化した話

前任者が残していった**レガシーコード**（C言語で2000行弱）をPythonで書き直してみたら、50行ほどになったでござる、という話です。

ここで言うレガシーコードとは、昔のプログラミング言語で書かれたコードという意味では無く、**変更が容易ではないコード（＝メンテナンスコストが高くつく）**という意味合いで用いています。[テストコード](https://qiita.com/tags/%E3%83%86%E3%82%B9%E3%83%88%E3%82%B3%E3%83%BC%E3%83%89)の無いコードをレガシーコードと呼ぶこともあるので注意。

# コードの処理概要
各IoTデバイスから5秒おきに流れてくるRAWデータ（MySQLに格納）を月別に集計し、ヒートマップを作成するプログラム。
ざっくり言ってそれだけ。

## 環境

- Ubuntu 18.04 (Amazon EC2)
- MySQL 5.7.33
- Python 3.7.3
 - pandas 1.1.5
 - seaborn 0.11.1

# 実装
## SQL
pandasでPIVOT化する前提で、SQLもシンプルに書き直しました。
MySQLにはDATE_TRUNC関数が無く、かつエポック秒で格納されているので、FROM_UNIXTIME関数で日時以下を切り捨てています。
DateTime型のカラムであれば、DATE_FORMAT関数で月別（または日別や時間別）に集計することになると思います。

```sql:SQL
SELECT
    device_id,                                      -- IoTデバイスID
    from_unixtime(rec_time, "%Y/%m") as dt,         -- UNIX時刻をyyyy/mm形式に変換
    sum(indicator_value) as iv_total,               -- 指標の集計
    day(last_day(from_unixtime(rec_time))) as lday  -- 当該月の日数
FROM sensor_table
GROUP BY device_id, dt, lday
ORDER BY device_id, dt;
```

```:結果
+-----------+---------+--------------------+------+
| device_id | dt      | iv_total           | lday |
+-----------+---------+--------------------+------+
| GB1736A   | 2020/11 | 13711366.145503998 |   30 |
| GB1736A   | 2020/12 | 14891288.108707428 |   31 |
| GB1736A   | 2021/01 |  3325523.908924103 |   31 |
| GB1736A   | 2021/02 | 10302060.682685852 |   28 |
| HS3912A   | 2020/11 |  5302481.319400787 |   30 |
| HS3912A   | 2020/12 | 2877198.6013946533 |   31 |
| HS3912A   | 2021/01 | 42210.279989242554 |   31 |
     :           :               :             :
```

１日当たりの合計を求めるため各月の日数（lday）も出しています。単純に30で割ると２月が少なくなってしまうからです。これはSQLじゃなくてプログラムで出しても良かったと思います。

年月の行を列に変換したいところですがMySQLにはピボット関数が無いため、力技でやろうとするとCASE文の羅列になってしまいます。
なので、Pythonの`pandas`モジュールにやらせます。

## pandasでPIVOT出力
`pandas.io.sql`でRDBとDataFrameを接続します。

```python:Python
import mysql.connector
import pandas as pd
import pandas.io.sql as psql

with mysql.connector.connect(
    host='localhost', port=3306, user='sensor_user', password='XXXXXXXX', database='sensor_db'
) as conn:
    sql = '''
        SELECT
            device_id,                                      -- IoTデバイスID
            from_unixtime(rec_time, "%Y/%m") as dt,         -- UNIX時刻をyyyy/mm形式に変換
            sum(indicator_value) as iv_total,               -- 指標の集計
            day(last_day(from_unixtime(rec_time))) as lday  -- 当該月の日数
        FROM sensor_table
        GROUP BY device_id, dt, lday
        ORDER BY device_id, dt;
    '''
    df = psql.read_sql(sql, conn)  # SQLの結果からDataFrameを作成

df['rt'] = df['iv_total'] / (24 * 60 * 12 * df['lday']) * 100  # 1日当たりの集計
df.drop(['iv_total', 'lday'], axis='columns', inplace=True)    # 余分な列を削除
df['rt'] = df['rt'].astype('int')                              # 整数に変換

print(df.pivot_table(index='device_id', columns='dt', fill_value=0))  # PIVOT化
```

まるまる１ヶ月間レコードが１件も追加されない、ということが起こり得たので、PIVOT化の際に`fill_value=0`で欠損値を穴埋めしています。
データベース接続パラメータは、`configparser`モジュール等で外だしした方が良いでしょう。

```:結果
               rt                                                  ...
dt        2018/04 2018/05 2018/06 2018/07 2018/08 2018/09 2018/10  ... 2020/08 2020/09 2020/10 2020/11 2020/12 2021/01 2021/02
device_id                                                          ...
GB2513A      5505    5487    5414    5478    5507    5539    5439  ...       0       0       0    2920    4797    1241       0
GB2657A      1695    3224       0    4329    1037    3236       0  ...       0    4676    5096       0       0       0       0
GB2660A      8213    8212    8197    8147    8140    8154    8200  ...    7986    8059    8059    8014    2389    7579    3408
GB2861A      8037    8006    8515    8592    8566    8561    8425  ...    8195    8272    8260    8347    8370    8331    2799
GB3250A      2511    3814       0       0     422    3257    4661  ...       0       0       0       0       0       0       0
GB3447A      5345    5200    5525    5340    5502    5608    5558  ...    5355    5458    5429    4862    5531    5559    2198
GB3665A      8441    8564    8525    8438    8486    8407    8399  ...    8467    8383    8419    8492    8516    8568    3435
GB4280A      5471    5570    5572    5336    5391    5504    5524  ...    5354    5409    5393    5411    5523    5535    2298
GB5971A      7970    7956    7960    7903    7966    7937    7968  ...    5419    4692    3659    3386    4743    6456    2326   
```

縦にID、横に年月が並ぶ表ができました。

## seabornでヒートマップ出力
[seaborn](https://seaborn.pydata.org)はPythonのデータ可視化ライブラリのひとつで、matplotlibと比べるとエレガントかつ少ないコード量で描画できるのでデータサイエンティストにも人気です。

先頭に次のコードを追加し、ライブラリをインポートします。
seabornはmatplotlibのラッパーモジュールなので、matplotlibも必要です。

```python:Python
import seaborn as sns
import matplotlib.pyplot as plt
```

末尾に次のコードを追加し、PIVOT化したDataFrameをseabornのheatmapメソッドに渡します。

```python:Python
plt.figure(figsize=(30, 15))  # 描画領域確保
sns.heatmap(df.pivot_table(index='device_id', columns='dt', fill_value=0), cmap='OrRd', annot=True, fmt='4d', linewidths=.5)
plt.title('Haetmap')        # タイトル
plt.xlabel('month')         # X軸のラベル
plt.ylabel('device id')     # Y軸のラベル
plt.savefig('heatmap.png')  # 画像保存
```

![heatmap2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/236222/121a743c-6162-f324-24fd-57a5ac10c286.png)
はい、これだけでヒートマップができました！　※データは架空のものです。
color map は[このサイト](https://pod.hatenablog.com/entry/2018/09/20/212527)を参考に決めました。

以上で目的は達成ですが、seabornには他にも便利なメソッドが沢山用意されているので一部を紹介しましょう。

### データ件数の棒グラフ - countplot
データの件数（頻度）を集計し、ヒストグラムとして出力します。基本中の基本ですね。

```python
sns.countplot(data=df, x='dt')
plt.xticks(rotation=45)  # X軸を45度傾ける
```

![countplot_x2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/236222/e1e451c4-a642-ecdf-e916-a15b766f13a4.png)
Y軸を指定すると横になります。

```python
sns.countplot(data=df, y='dt')
```

![countplot_y2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/236222/f6d810e9-8d68-bdd6-be8d-c677f79c9120.png)

### ヒストグラム - histplot
distplotメソッドが *deprecated* になったようで、本メソッドが推奨になりました。
重なる滑らかな曲線（カーネル密度推定）は`kde=False`で表示されなくなります。

```python
sns.histplot(df['rt'], bins=25, kde=True)
```
![histplot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/236222/36312ab6-c28f-ef40-0e4b-bea6b9a02d48.png)

### 散布図 - scatterplot
X軸とY軸の2つの項目で量を計測し、分布を表現するものです。

```python
sns.scatterplot(data=df.query('device_id.str.startswith("NM")'), x='dt', y='rt', hue='device_id')
plt.xticks(rotation=45)  # X軸を45度傾ける
```

2つの項目が原因と結果を表していれば、X軸に原因、Y軸に結果を与えることで相関関係が分かります。
下図では、指標と年月の間にそのような関係性は無いため、例題としては良くありません。
![scatterplot2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/236222/8766928c-cd1b-40a6-46c9-98aa9731021d.png)

# さいごに
**軽量プログラミング言語（Lightweight Language）**の中でもPythonは特に扱い易く、著名なライブラリを用いてシンプルにリファクタリングできました。いわゆる**技術負債の返済**が終わり、今後はプログラマの調達が容易になって低コスト＆スピード感を持って対応していけそうです。（C言語経験者も集め易いけど比較的ベテランさんが多いためか高単価になりがちな気がする）

なお、このIoTデバイスはNMEA0183を出力します。

```
$GPRMC,092047.000,V,0251.7300,N,00659.1050,E,,,101218,,
```

|カラム番号|内容|
|:-:|:--|
|1|$GPRMC
|2|協定世界時(UTC)での時刻。hhmmss.ss
|3|ステータス データ有効:A データ無効:V
|4|緯度 AABB.CCDD (0度の場合は0BB.CC.DD)
|5|北緯:N 南緯:S
|6|経度 AABB.CCDD (0度の場合は0BB.CC.DD)
|7|東経:E 西経:S
|8|対地速度 (knot)
|9|進行方向 (度)
|10|協定世界時(UTC)での日付。ddmmyy
|11|磁気偏差量 (度)
|12|磁気偏差方向 西偏:W 東偏:E

GPS位置情報を使い、GoogleMapやOpenStreetMapの地図上にヒートマップを重ねて表示できそうでワクワクします。
イメージとしては[こんな感じ](https://github.com/njanakiev/openstreetmap-heatmap)でやりたい。
いずれ取り組む機会があったら記事にしたいと思います。