前任者が残していったレガシーコード（C言語で2000行弱）をPythonで書き直してみたら、50行ほどになったでござる、という話です。
ここで言うレガシーコードとは、昔のプログラミング言語で書かれたコードという意味では無く、変更が容易ではないコード（＝メンテナンスコストが高くつく）という意味合いで用いています。テストコードの無いコードをレガシーコードと呼ぶこともあるので注意。
コードの処理概要
各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関数で月別（または日別や時間別）に集計することになると思います。
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を接続します。
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はPythonのデータ可視化ライブラリのひとつで、matplotlibと比べるとエレガントかつ少ないコード量で描画できるのでデータサイエンティストにも人気です。
先頭に次のコードを追加し、ライブラリをインポートします。
seabornはmatplotlibのラッパーモジュールなので、matplotlibも必要です。
import seaborn as sns
import matplotlib.pyplot as plt
末尾に次のコードを追加し、PIVOT化したDataFrameをseabornのheatmapメソッドに渡します。
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') # 画像保存
はい、これだけでヒートマップができました！ ※データは架空のものです。
color map はこのサイトを参考に決めました。
以上で目的は達成ですが、seabornには他にも便利なメソッドが沢山用意されているので一部を紹介しましょう。
データ件数の棒グラフ - countplot
データの件数（頻度）を集計し、ヒストグラムとして出力します。基本中の基本ですね。
sns.countplot(data=df, x='dt')
plt.xticks(rotation=45) # X軸を45度傾ける
Y軸を指定すると横になります。
sns.countplot(data=df, y='dt')
ヒストグラム - histplot
distplotメソッドが deprecated になったようで、本メソッドが推奨になりました。
重なる滑らかな曲線（カーネル密度推定）はkde=Falseで表示されなくなります。
sns.histplot(df['rt'], bins=25, kde=True)
散布図 - scatterplot
X軸とY軸の2つの項目で量を計測し、分布を表現するものです。
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軸に結果を与えることで相関関係が分かります。
下図では、指標と年月の間にそのような関係性は無いため、例題としては良くありません。
さいごに
軽量プログラミング言語（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の地図上にヒートマップを重ねて表示できそうでワクワクします。
イメージとしてはこんな感じでやりたい。
いずれ取り組む機会があったら記事にしたいと思います。