本PBLではコインランドリーの稼働履歴データと気象データを使い、気象条件がコインランドリー利用者数に与える影響を分析します。
またその結果をもとに、気象条件によって起こる稼働のピークを平滑化する施策を立案します。
分析には、ウォッシュプラス富士見1丁目店(千葉県浦安市)のデータを利用します。
import os
import glob
import numpy as np
import pandas as pd
# from datetime import datetime as dt
# from datetime import timedelta
pd.set_option('display.max_columns',None)
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
### matplotlibの日本語対応 ###
#---日本語設定をしていない場合は、以下のコメントアウトを外して実行して下さい。
# → matplotlibのバージョンによっては、以下の設定をしても日本語が表示されない場合があります。
# from matplotlib import rcParams
# rcParams['font.family'] = 'sans-serif'
# rcParams['font.sans-serif'] = ['Hiragino Maru Gothic Pro', 'Yu Gothic', 'Meirio']
washplus
フォルダにある、稼働履歴データを読み込みます。
稼働履歴データのファイルは通常のCSV形式なので、Pandasのread_csv
関数で簡単に読み込めますが、ファイル数が複数あります。
ウォーミングアップとしてまずファイルを1つ、どれでも良いので読み込んでみましょう。
ポイントとして文字コードがshift-jis
のため、read_csv関数のオプションとして指定する必要があります。
読み込んだら最初の5行を表示させてみましょう。
#--稼働履歴データを1つ読み込む
tmp = pd.read_csv("washplus/machine_operations_0002_all_20171201_20180228.csv", encoding="shift-jis")
tmp.head()
開始日時 | 終了日時 | 取出し日時 | 号機 | 機器種別 | コース番号 | コース区分 | 顧客ID | OS | 会員グループ名 | 料金 | チケット利用枚数 | 支払料金 | アプリ決済種別 | 決済備考 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018/02/01 11:47:57 | 2018/02/01 12:47:57 | 2018/02/01 14:51:27 | 2 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 600 | NaN | NaN |
1 | 2018/02/01 11:52:37 | 2018/02/01 12:22:47 | 2018/02/01 12:22:57 | 3 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 300 | NaN | NaN |
2 | 2018/02/01 12:03:57 | 2018/02/01 12:33:57 | 2018/02/01 12:48:17 | 1 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 300 | NaN | NaN |
3 | 2018/02/01 12:06:47 | 2018/02/01 12:26:57 | 2018/02/01 12:27:07 | 5 | 洗濯乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 200 | NaN | NaN |
4 | 2018/02/01 12:06:52 | 2018/02/01 12:16:52 | 2018/02/01 12:28:22 | 6 | 洗濯乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 100 | NaN | NaN |
次に全ファイルをまとめて読み込みます。
ここではglob
関数を使って全ファイル名を取得したのち、forループで順番に読み込んで、concat
関数で縦につなげる処理を行います。
こうすることで、全データを1つの変数に格納します。
#---全ファイルをまとめて読み込む
# ファイル名のリストを取得
dirname = "washplus"
files = glob.glob( os.path.join(dirname,"*.csv") )
print("読み込むファイル")
# 日本語を含むCSVファイルなので、文字コードShift-JISを指定して読み込む
laundry_org = pd.DataFrame()
for file in files:
print(file)
tmp = pd.read_csv(file, encoding="shift-jis")
laundry_org = pd.concat([laundry_org,tmp])
# 読み込んだデータを確認
laundry_org.head()
読み込むファイル washplus/machine_operations_0002_all_20200301_20200531.csv washplus/machine_operations_0002_all_20191201_20200229.csv washplus/machine_operations_0002_all_20180301_20180531.csv washplus/machine_operations_0002_all_20190301_20190531.csv washplus/machine_operations_0002_all_20190601_20190831.csv washplus/machine_operations_0002_all_20200901_20200930.csv washplus/machine_operations_0002_all_20190901_20191130.csv washplus/machine_operations_0002_all_20200601_20200831.csv washplus/machine_operations_0002_all_20171201_20180228.csv washplus/machine_operations_0002_all_20181201_20190228.csv washplus/machine_operations_0002_all_20180901_20181130.csv washplus/machine_operations_0002_all_20180601_20180831.csv
開始日時 | 終了日時 | 取出し日時 | 号機 | 機器種別 | コース番号 | コース区分 | 顧客ID | OS | 会員グループ名 | 料金 | チケット利用枚数 | 支払料金 | アプリ決済種別 | 決済備考 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2020/03/01 01:21:41 | 2020/03/01 01:51:51 | 2020/03/01 01:52:31 | 2 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 300 | NaN | NaN |
1 | 2020/03/01 01:22:11 | 2020/03/01 01:32:11 | 2020/03/01 01:51:41 | 4 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 100 | NaN | NaN |
2 | 2020/03/01 02:19:56 | 2020/03/01 03:21:21 | 2020/03/01 03:25:11 | 5 | 洗濯乾燥機 | 1 | NaN | NaN | NaN | NaN | NaN | NaN | 1100 | NaN | NaN |
3 | 2020/03/01 04:19:21 | 2020/03/01 04:49:21 | 2020/03/01 04:50:01 | 1 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 300 | NaN | NaN |
4 | 2020/03/01 05:04:27 | 2020/03/01 06:06:16 | 2020/03/01 06:12:56 | 7 | 洗濯乾燥機 | 1 | NaN | NaN | NaN | NaN | NaN | NaN | 1100 | NaN | NaN |
読み込んだファイル名を見ると、日付順になっていないかもしれません。
その場合、変数laundry_orgの中身も日付順になっていませんので、日付の古い順に並び替える処理をします。
ここでは、Pandasのto_datetime
関数で開始日時をTimestamp型に変換します。
Timestamp型にすると、年・月・日・時・曜日などの情報に分解でき、また日付の足し算(e.g. 1時間後を計算)などができて便利です。
Timestamp型に変換したら、sort_values
関数で開始日時の古い順(降順)にソートします。
#---開始日時でソート
# 開始日時をTimestamp型に変換して日付順にソートする
laundry_org["datetime"] = pd.to_datetime(laundry_org["開始日時"])
laundry_org = laundry_org.sort_values("datetime").reset_index(drop=True)
laundry_org.head()
開始日時 | 終了日時 | 取出し日時 | 号機 | 機器種別 | コース番号 | コース区分 | 顧客ID | OS | 会員グループ名 | 料金 | チケット利用枚数 | 支払料金 | アプリ決済種別 | 決済備考 | datetime | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018/02/01 11:47:57 | 2018/02/01 12:47:57 | 2018/02/01 14:51:27 | 2 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 600 | NaN | NaN | 2018-02-01 11:47:57 |
1 | 2018/02/01 11:52:37 | 2018/02/01 12:22:47 | 2018/02/01 12:22:57 | 3 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 300 | NaN | NaN | 2018-02-01 11:52:37 |
2 | 2018/02/01 12:03:57 | 2018/02/01 12:33:57 | 2018/02/01 12:48:17 | 1 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 300 | NaN | NaN | 2018-02-01 12:03:57 |
3 | 2018/02/01 12:06:47 | 2018/02/01 12:26:57 | 2018/02/01 12:27:07 | 5 | 洗濯乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 200 | NaN | NaN | 2018-02-01 12:06:47 |
4 | 2018/02/01 12:06:52 | 2018/02/01 12:16:52 | 2018/02/01 12:28:22 | 6 | 洗濯乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 100 | NaN | NaN | 2018-02-01 12:06:52 |
# 末尾も見てみましょう
laundry_org.tail()
開始日時 | 終了日時 | 取出し日時 | 号機 | 機器種別 | コース番号 | コース区分 | 顧客ID | OS | 会員グループ名 | 料金 | チケット利用枚数 | 支払料金 | アプリ決済種別 | 決済備考 | datetime | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
53047 | 2020/09/30 22:53:08 | 2020/09/30 23:13:13 | 2020/09/30 23:31:03 | 6 | 洗濯乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 200 | NaN | NaN | 2020-09-30 22:53:08 |
53048 | 2020/09/30 22:56:03 | 2020/09/30 23:16:03 | 2020/09/30 23:16:03 | 4 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 200 | NaN | NaN | 2020-09-30 22:56:03 |
53049 | 2020/09/30 23:18:23 | 2020/09/30 23:38:23 | 2020/09/30 23:38:33 | 4 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 200 | NaN | NaN | 2020-09-30 23:18:23 |
53050 | 2020/09/30 23:33:28 | 2020/09/30 23:43:38 | 2020/09/30 23:44:23 | 5 | 洗濯乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 100 | NaN | NaN | 2020-09-30 23:33:28 |
53051 | 2020/09/30 23:33:58 | 2020/09/30 23:44:02 | 2020/09/30 23:44:18 | 6 | 洗濯乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 100 | NaN | NaN | 2020-09-30 23:33:58 |
2018年2月〜2020年9月の稼働履歴データを読み込みました。
さっそく幾つかの切り口で集計をしてみましょう。
カラム名を見ると、機器種別・号機・コース番号などで集計できそうです。また稼働回数だけでなく、売上の計算もできそうです。
|号機|機械種別| |--|--| |1|乾燥機(二段)| |2|乾燥機(二段)| |3|乾燥機(二段)| |4|乾燥機(二段)| |5|洗濯乾燥機(中型)| |6|洗濯乾燥機(中型)| |7|洗濯乾燥機(中型)| |8|洗濯乾燥機(大型)| |9|洗濯乾燥機(大型)| |10|洗濯乾燥機(特大)| | |コース番号|コース名| |--|--| |1|洗濯乾燥(標準)| |2|洗濯乾燥(少なめ)| |3|洗濯| |4|乾燥| |
各カラムについては付属資料「washplusデータ基本情報.xlsx
」も参照して下さい。
ここでは気象との関連が大きそうなコース番号を中心に見ていきます。
まず稼働回数から集計します。
value_counts
関数を使うと、指定したカラム(e.g. 機器種別)の構成要素(e.g. 洗濯乾燥機, 乾燥機)について出現回数をカウントできます。
# 機器種別で稼働回数を集計
laundry_org["機器種別"].value_counts()
洗濯乾燥機 28023 乾燥機 25029 Name: 機器種別, dtype: int64
# コース番号で稼働回数を集計
laundry_org["コース番号"].value_counts()
4 33487 1 10697 2 5666 3 3202 Name: コース番号, dtype: int64
乾燥コースの利用頻度が最も多く、次いで洗濯乾燥(標準)コースが利用されています。
次に支払料金カラムを使って、売上を計算します。
groupby
関数を使うと、指定したカラム(e.g. 機器種別)の構成要素(e.g. 洗濯乾燥機, 乾燥機)でグルーピングし、集計したいカラム(e.g. 支払料金)についてグルーピングした構成要素ごとの合計(sum
関数)を計算することができます。
# 機器種別ごとに支払料金を集計
laundry_org.groupby(["機器種別"])[["支払料金"]].sum()
支払料金 | |
---|---|
機器種別 | |
乾燥機 | 7098800 |
洗濯乾燥機 | 24877000 |
# コース番号ごとに支払料金を集計
laundry_org.groupby(["コース番号"])[["支払料金"]].sum()
支払料金 | |
---|---|
コース番号 | |
1 | 14648300 |
2 | 5839900 |
3 | 2241500 |
4 | 9246100 |
# 機器種別とコース番号ごとに支払料金を集計
laundry_org.groupby(["機器種別","コース番号"])[["支払料金"]].sum()
支払料金 | ||
---|---|---|
機器種別 | コース番号 | |
乾燥機 | 4 | 7098800 |
洗濯乾燥機 | 1 | 14648300 |
2 | 5839900 | |
3 | 2241500 | |
4 | 2147300 |
稼働回数とは逆に、洗濯乾燥(標準)コースの方が売上は大きくなっています。
稼働回数は乾燥コースより少なくても、単価が大きい洗濯乾燥(標準)コースの方が売上は大きくなるようです。
ざっくりと集計してみましたので、今度は日付情報を使って月ごとや週ごとに利用者数(稼働回数)の特徴を見ていきます。
事前準備として、日付情報を年・月・日・時・曜日に分けて格納します。
上で開始日時をTimestamp型に変換してあるので、年・月・日・時・曜日の情報は以下のようにして簡単に取り出せます。
また月ごとに集計する都合上、ちょうど2年分のデータに絞ります。今回は2018年9月〜2020年8月を使います。
#--Timestamp型で日時の処理
# dateはYYYY-MM-DDのフォーマットに日付を変換しています
# weekdayは0〜6の数字で、それぞれ月曜〜日曜を意味します
laundry_org["date"] = laundry_org["datetime"].dt.strftime("%Y-%m-%d")
laundry_org["year"] = laundry_org["datetime"].dt.year
laundry_org["month"] = laundry_org["datetime"].dt.month
laundry_org["day"] = laundry_org["datetime"].dt.day
laundry_org["hour"] = laundry_org["datetime"].dt.hour
laundry_org["weekday"] = laundry_org["datetime"].dt.weekday
# 指定した日付の範囲を取り出す
laundry_org = laundry_org.query('"2018-09-01" <= date < "2020-09-01"')
laundry_org.head()
開始日時 | 終了日時 | 取出し日時 | 号機 | 機器種別 | コース番号 | コース区分 | 顧客ID | OS | 会員グループ名 | 料金 | チケット利用枚数 | 支払料金 | アプリ決済種別 | 決済備考 | datetime | date | year | month | day | hour | weekday | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
10357 | 2018/09/01 06:59:36 | 2018/09/01 07:53:06 | 2018/09/01 07:57:11 | 5 | 洗濯乾燥機 | 2 | NaN | NaN | NaN | NaN | NaN | NaN | 900 | NaN | NaN | 2018-09-01 06:59:36 | 2018-09-01 | 2018 | 9 | 1 | 6 | 5 |
10358 | 2018/09/01 07:17:16 | 2018/09/01 07:47:16 | 2018/09/01 07:58:36 | 1 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 300 | NaN | NaN | 2018-09-01 07:17:16 | 2018-09-01 | 2018 | 9 | 1 | 7 | 5 |
10359 | 2018/09/01 07:24:36 | 2018/09/01 07:44:36 | 2018/09/01 07:45:06 | 4 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 200 | NaN | NaN | 2018-09-01 07:24:36 | 2018-09-01 | 2018 | 9 | 1 | 7 | 5 |
10360 | 2018/09/01 07:41:06 | 2018/09/01 08:33:11 | 2018/09/01 08:37:56 | 6 | 洗濯乾燥機 | 2 | NaN | NaN | NaN | NaN | NaN | NaN | 900 | NaN | NaN | 2018-09-01 07:41:06 | 2018-09-01 | 2018 | 9 | 1 | 7 | 5 |
10361 | 2018/09/01 08:05:16 | 2018/09/01 08:45:16 | 2018/09/01 08:52:26 | 4 | 乾燥機 | 4 | NaN | NaN | NaN | NaN | NaN | NaN | 400 | NaN | NaN | 2018-09-01 08:05:16 | 2018-09-01 | 2018 | 9 | 1 | 8 | 5 |
今後、何回か積み上げ棒グラフで可視化するため、まず関数を定義しておきます。
# 積み上げ棒グラフを描画する関数
# 引数
# - data:クロス集計したデータ(データフレーム)
# - xlabl:グラフのx軸のラベル
def bar_plot(data, xlabel="Index"):
fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(1,1,1)
x = data.index
for i, t in enumerate(data.columns):
y = data.iloc[:, i:].sum(axis=1)
ax.bar(x, y, label="コース"+str(t))
ax.set_xticks( x )
ax.set_xlabel(xlabel, fontsize=15)
ax.set_ylabel('稼働回数', fontsize=15)
ax.legend()
plt.show()
次に、月ごと・コースごとにクロス集計し、各コースの稼働回数が月(季節)によりどう変化するのか見てみます。
crosstab
関数を使うとExcelのピボットテーブルのような集計(出現回数のカウント)ができます。
以下のセルでは、クロス集計結果を上で定義したbar_plot
関数に与えてグラフを作成し、クロス集計結果そのものも出力させています。
#---月ごと・コースごとにクロス集計
# crosstab関数を使う
laundry_monthly = pd.crosstab(laundry_org["month"], laundry_org["コース番号"])
bar_plot(laundry_monthly, "月")
# クロス集計結果を表示
laundry_monthly
コース番号 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
month | ||||
1 | 630 | 290 | 161 | 1885 |
2 | 561 | 289 | 174 | 1626 |
3 | 504 | 317 | 189 | 1985 |
4 | 588 | 285 | 193 | 1643 |
5 | 789 | 346 | 221 | 1725 |
6 | 625 | 327 | 192 | 2112 |
7 | 780 | 466 | 179 | 4033 |
8 | 747 | 435 | 281 | 1549 |
9 | 683 | 430 | 183 | 2425 |
10 | 705 | 394 | 176 | 2385 |
11 | 733 | 311 | 159 | 2269 |
12 | 808 | 380 | 200 | 2596 |
乾燥コースは7月が突出して多いことがわかります。ついで12月が多くなっています。
洗濯乾燥(標準)コースは乾燥コースに比べると月ごとの差は小さいものの、12月・7月の順に多くなっています。
これらの月で稼働回数が多いのは何故でしょうか?
7月は梅雨、12月は年末であることが関係していそうですね。
続いて曜日ごと・コースごとにクロス集計し、曜日による利用者数の変化を見てみましょう。
#---曜日ごと・コースごとにクロス集計
# crosstab関数を使う
laundry_weekday = pd.crosstab(laundry_org["weekday"], laundry_org["コース番号"])
laundry_weekday.index = ["月","火","水","木","金","土","日"]
bar_plot(laundry_weekday, "曜日")
# クロス集計結果を表示
laundry_weekday
コース番号 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
月 | 1118 | 524 | 314 | 3557 |
火 | 890 | 522 | 276 | 3166 |
水 | 1045 | 564 | 295 | 3331 |
木 | 875 | 527 | 321 | 3089 |
金 | 1028 | 538 | 284 | 3232 |
土 | 1477 | 719 | 364 | 4620 |
日 | 1720 | 876 | 454 | 5238 |
各コースとも土日、特に日曜日の稼働回数(利用者数)が多いことがわかります。
今度は時間ごと・コースごとにクロス集計し、時間帯による利用者数の変化を見てみましょう。
#---時間ごと・コースごとにクロス集計
# crosstab関数を使う
laundry_hourly = pd.crosstab(laundry_org["hour"], laundry_org["コース番号"])
bar_plot(laundry_hourly, "時刻")
# クロス集計結果を表示
laundry_hourly
コース番号 | 1 | 2 | 3 | 4 |
---|---|---|---|---|
hour | ||||
0 | 226 | 153 | 85 | 701 |
1 | 190 | 75 | 71 | 574 |
2 | 112 | 67 | 43 | 324 |
3 | 93 | 48 | 34 | 236 |
4 | 144 | 61 | 41 | 252 |
5 | 191 | 71 | 52 | 581 |
6 | 166 | 85 | 89 | 781 |
7 | 182 | 103 | 82 | 1263 |
8 | 355 | 149 | 91 | 1387 |
9 | 523 | 206 | 141 | 1454 |
10 | 635 | 233 | 181 | 1620 |
11 | 572 | 213 | 146 | 1577 |
12 | 476 | 207 | 129 | 1304 |
13 | 469 | 237 | 92 | 1425 |
14 | 445 | 188 | 96 | 1379 |
15 | 437 | 204 | 108 | 1468 |
16 | 418 | 227 | 105 | 1451 |
17 | 399 | 281 | 100 | 1448 |
18 | 414 | 267 | 117 | 1409 |
19 | 385 | 265 | 117 | 1257 |
20 | 365 | 283 | 97 | 1215 |
21 | 306 | 239 | 78 | 1193 |
22 | 355 | 217 | 98 | 986 |
23 | 295 | 191 | 115 | 948 |
全体として午前中に増加する稼働回数は10時ごろにピークに達し、夕方ごろまで多い状態が続いたのち、19時ごろから減少する傾向があります。
コースごとに見ても、洗濯乾燥(標準)コースと、乾燥コースは全体と同じような傾向があります。
一方で洗濯乾燥(少なめ)コースに関しては、午前10時ごろよりも夕方の方が稼働回数が多くなっています。
ここまでの分析結果を簡単にまとめます。
#---日付ごと・コースごとにクロス集計
# crosstab関数を使う
laundry_daily = pd.crosstab(laundry_org["date"], laundry_org["コース番号"])
laundry_daily.columns = [f"コース{i}" for i in laundry_daily.columns]
laundry_daily = laundry_daily.reset_index()
laundry_daily["date"] = pd.to_datetime(laundry_daily["date"])
laundry_daily.head()
date | コース1 | コース2 | コース3 | コース4 | |
---|---|---|---|---|---|
0 | 2018-09-01 | 7 | 10 | 3 | 24 |
1 | 2018-09-02 | 14 | 14 | 3 | 77 |
2 | 2018-09-03 | 6 | 2 | 1 | 79 |
3 | 2018-09-04 | 12 | 6 | 2 | 42 |
4 | 2018-09-05 | 11 | 11 | 1 | 30 |
次に気象データを読み込みます。気象データはweather
フォルダの中に格納されています。
今回はアメダス江戸川臨海の観測データを、気象庁ホームページより事前にダウンロードしておいたものを利用します。
ダウンロードしたファイルはCSV形式ですが、ヘッダが最初の行(0行目)にないなど、特殊なフォーマットとなっています。
そのため、うまく読み込むためには工夫が必要です。
ここでは以下のように、アメダスCSVデータを読み込むための関数を定義して使用します。
#---アメダスCSVファイルを読み込む関数
# 引数
# - file : ファイル名
# - all_params : ファイルに含まれるカラムの名前リスト
# - params : 最終的に欲しいカラムの名前リスト
# - start_index : データの始まりの行(最初の行は0行目と数える)
def Read_AMeDAS_CSV(file, all_params, params, start_index=6):
# カラム名を指定してファイルを読み込む
# データがスタートする行を指定する
data = pd.read_csv(file, encoding="shift-jis", names=all_params, skiprows=range(start_index))
# 品質番号が8or5以外だったら、品質が不十分としてその気象要素の値を欠損値にする
for col in params:
if col + "_q" in all_params:
data[col] = data[col].mask(data[col+"_q"] < 5, np.nan)
# 必要なカラムだけ残して返す
return data[params]
ダウンロードしたアメダスCSVにはカラム名として使える適切な行がないため、あらかじめ用意したカラム名リストを渡して読み込みます。
それが引数all_paramsです。
事前に用意したアメダスCSVデータは、品質情報と均質番号を含んだ状態で入手しました。
今回は品質情報が8
または5
の場合は良しとしますが、5未満
の場合は欠損値として扱うことにします。
また均質番号は期間通して1
であることを確認済みとし、今回は使用しません。
品質番号はカラム名を気象要素名_q
、均質番号は気象要素名_i
とすることで識別できるようにカラム名を付けます。
読み込んだデータのうち、均質情報などは読み込み時にしか使いませんので、関数からreturnする際に落とすことにします。
必要なカラム名のみのリストを渡してその処理をします。それが引数paramsです。
それではアメダス江戸川臨海の日別観測値を読み込みましょう。
#---江戸川臨海の日別観測値
# 全カラム名リスト
all_columns1 = ["年月日","平均気温","平均気温_q","平均気温_i",
"最高気温","最高気温_q","最高気温_i",
"最低気温","最低気温_q","最低気温_i",
"降水量","降水量_q","降水量_i",
"日照時間","日照時間_q","日照時間_i",
"最大風速","最大風速_q","最大風速の風向","最大風速の風向_q","最大風速_i",
"平均風速","平均風速_q","平均風速_i",
"最多風向","最多風向_q","最多風向_i",
"最大瞬間風速","最大瞬間風速_q","最大瞬間風速の風向","最大瞬間風速の風向_q","最大瞬間風速_i",
"最大1時間降水量","最大1時間降水量_q","最大1時間降水量_i"]
# 必要なカラム名リスト
select_columns1 = ["年月日","平均気温","最高気温","最低気温","降水量","日照時間",
"最大風速","最大風速の風向","平均風速","最多風向",
"最大瞬間風速","最大瞬間風速の風向","最大1時間降水量"]
# 読み込み
file1 = "weather/edogawarinkai_daily_20180101_20200930.csv"
obs1 = Read_AMeDAS_CSV(file1, all_columns1, select_columns1)
obs1.head()
年月日 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018/1/1 | 6.7 | 13.3 | 0.8 | 0.0 | 9.3 | 4.9 | 北北西 | 2.5 | 北北西 | 10.6 | 北北西 | 0.0 |
1 | 2018/1/2 | 7.1 | 10.6 | 1.7 | 0.0 | 9.4 | 6.7 | 北西 | 3.2 | 北北西 | 12.0 | 北西 | 0.0 |
2 | 2018/1/3 | 5.3 | 8.7 | 1.6 | 0.0 | 9.5 | 9.1 | 北西 | 5.1 | 北西 | 17.1 | 西北西 | 0.0 |
3 | 2018/1/4 | 4.8 | 9.5 | -1.0 | 0.0 | 9.4 | 5.3 | 北 | 2.9 | 北 | 9.7 | 北 | 0.0 |
4 | 2018/1/5 | 3.5 | 6.1 | 0.6 | 0.0 | 0.0 | 4.1 | 北北西 | 2.2 | 北 | 7.9 | 北北西 | 0.0 |
アメダス江戸川臨海は4要素(気温・風向風速・降水量・日照時間)の観測はありますが、相対湿度の観測がありません。
しかしながら、コインランドリーと相対湿度には何らかの相関があることも考えられ、できれば相対湿度も使いたいです。
そこで東京管区気象台の観測データから相対湿度だけダウンロードしました。これを利用することにします。
#---東京管区気象台の観測値
# 全カラム名リスト
all_columns2 = ["年月日","平均湿度","平均湿度_q","平均湿度_i",
"最小相対湿度","最小相対湿度_q","最小相対湿度_i"]
# 必要なカラム名リスト
select_columns2 = ["年月日","平均湿度","最小相対湿度"]
# 読み込み
file2 = "weather/tokyo_daily_20180101_20200930.csv"
obs2 = Read_AMeDAS_CSV(file2, all_columns2, select_columns2)
obs2.head()
年月日 | 平均湿度 | 最小相対湿度 | |
---|---|---|---|
0 | 2018/1/1 | 55.0 | 27.0 |
1 | 2018/1/2 | 42.0 | 24.0 |
2 | 2018/1/3 | 43.0 | 21.0 |
3 | 2018/1/4 | 41.0 | 25.0 |
4 | 2018/1/5 | 60.0 | 38.0 |
アメダス江戸川臨海と東京管区気象台のデータを結合します。
merge
関数を使い、年月日カラムを軸にして、東京管区気象台のデータをアメダス江戸川臨海に結合します。
#---合成
obs_daily = pd.merge(obs1, obs2, on="年月日", how="left")
obs_daily.head()
年月日 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018/1/1 | 6.7 | 13.3 | 0.8 | 0.0 | 9.3 | 4.9 | 北北西 | 2.5 | 北北西 | 10.6 | 北北西 | 0.0 | 55.0 | 27.0 |
1 | 2018/1/2 | 7.1 | 10.6 | 1.7 | 0.0 | 9.4 | 6.7 | 北西 | 3.2 | 北北西 | 12.0 | 北西 | 0.0 | 42.0 | 24.0 |
2 | 2018/1/3 | 5.3 | 8.7 | 1.6 | 0.0 | 9.5 | 9.1 | 北西 | 5.1 | 北西 | 17.1 | 西北西 | 0.0 | 43.0 | 21.0 |
3 | 2018/1/4 | 4.8 | 9.5 | -1.0 | 0.0 | 9.4 | 5.3 | 北 | 2.9 | 北 | 9.7 | 北 | 0.0 | 41.0 | 25.0 |
4 | 2018/1/5 | 3.5 | 6.1 | 0.6 | 0.0 | 0.0 | 4.1 | 北北西 | 2.2 | 北 | 7.9 | 北北西 | 0.0 | 60.0 | 38.0 |
これで日別の観測値を読み込むことができました。しかし、これだけでは晴れ・曇り・雨の天気カテゴリーが不明です。
東京管区気象台の天気の観測を使う方法もありますが、今回はせっかくアメダス江戸川臨海が近くにあるので、日照時間や降水量の観測から天気カテゴリーを推定する方法をとります。
(なお天気カテゴリーは、推計気象分布から入手する方法もあります。)
以下のようなロジックで天気を推定します。
(注意)
この方法で天気を推定した場合、例えば「早朝まで雨が降っていたが日中は晴れた」ケースであっても日天気が「雨」と判定されてしまう
ではまず観測日(年月日)をTimestamp型にして月を抽出し、月でグルーピングし、月ごとの日照時間の最大値を求めます。
obs_daily["date"] = pd.to_datetime(obs_daily["年月日"])
obs_daily.drop(columns="年月日", inplace=True)
obs_daily["month"] = obs_daily["date"].dt.month
max_sun_hour = obs_daily.groupby("month")[["日照時間"]].max().reset_index()
max_sun_hour.columns = ["month","日照時間max"]
max_sun_hour
month | 日照時間max | |
---|---|---|
0 | 1 | 10.0 |
1 | 2 | 10.7 |
2 | 3 | 11.8 |
3 | 4 | 12.8 |
4 | 5 | 13.5 |
5 | 6 | 13.8 |
6 | 7 | 13.7 |
7 | 8 | 12.9 |
8 | 9 | 11.8 |
9 | 10 | 11.1 |
10 | 11 | 10.3 |
11 | 12 | 9.6 |
概ね違和感のない値が求まりましたので、monthカラムを軸として、変数max_sun_hourを変数obs_dailyにマージします。
そうすると、変数obs_dailyに日照時間maxカラムが結合されます。
次に日照時間と日照時間maxの比を計算し、日照時間rateというカラムを作ります。
obs_daily = pd.merge(obs_daily, max_sun_hour, on="month", how="left")
obs_daily["日照時間rate"] = obs_daily["日照時間"] / obs_daily["日照時間max"]
obs_daily.head()
平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | date | month | 日照時間max | 日照時間rate | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 6.7 | 13.3 | 0.8 | 0.0 | 9.3 | 4.9 | 北北西 | 2.5 | 北北西 | 10.6 | 北北西 | 0.0 | 55.0 | 27.0 | 2018-01-01 | 1 | 10.0 | 0.93 |
1 | 7.1 | 10.6 | 1.7 | 0.0 | 9.4 | 6.7 | 北西 | 3.2 | 北北西 | 12.0 | 北西 | 0.0 | 42.0 | 24.0 | 2018-01-02 | 1 | 10.0 | 0.94 |
2 | 5.3 | 8.7 | 1.6 | 0.0 | 9.5 | 9.1 | 北西 | 5.1 | 北西 | 17.1 | 西北西 | 0.0 | 43.0 | 21.0 | 2018-01-03 | 1 | 10.0 | 0.95 |
3 | 4.8 | 9.5 | -1.0 | 0.0 | 9.4 | 5.3 | 北 | 2.9 | 北 | 9.7 | 北 | 0.0 | 41.0 | 25.0 | 2018-01-04 | 1 | 10.0 | 0.94 |
4 | 3.5 | 6.1 | 0.6 | 0.0 | 0.0 | 4.1 | 北北西 | 2.2 | 北 | 7.9 | 北北西 | 0.0 | 60.0 | 38.0 | 2018-01-05 | 1 | 10.0 | 0.00 |
日照時間rateと降水量から天気カテゴリーを作ります。
Estimate_Weather
を定義しますapply
関数とlambda
を組み合わせて、変数obs_dailyのデータフレーム1行ごとにEstimate_Weather
を適用します#---天気カテゴリーの推定をする関数
# 引数
# - sun_rate : 日照時間rate
# - prcp : 降水量
def Estimate_Weather(sun_rate, prcp):
if prcp > 0.0:
return "Rainy"
elif sun_rate < 0.2:
return "Cloudy"
else:
return "Sunny"
obs_daily["天気"] = obs_daily[["日照時間rate","降水量"]].apply(lambda x: Estimate_Weather(x[0],x[1]), axis=1)
obs_daily.drop(columns=["month","日照時間max"], inplace=True)
obs_daily.head()
平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | date | 日照時間rate | 天気 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 6.7 | 13.3 | 0.8 | 0.0 | 9.3 | 4.9 | 北北西 | 2.5 | 北北西 | 10.6 | 北北西 | 0.0 | 55.0 | 27.0 | 2018-01-01 | 0.93 | Sunny |
1 | 7.1 | 10.6 | 1.7 | 0.0 | 9.4 | 6.7 | 北西 | 3.2 | 北北西 | 12.0 | 北西 | 0.0 | 42.0 | 24.0 | 2018-01-02 | 0.94 | Sunny |
2 | 5.3 | 8.7 | 1.6 | 0.0 | 9.5 | 9.1 | 北西 | 5.1 | 北西 | 17.1 | 西北西 | 0.0 | 43.0 | 21.0 | 2018-01-03 | 0.95 | Sunny |
3 | 4.8 | 9.5 | -1.0 | 0.0 | 9.4 | 5.3 | 北 | 2.9 | 北 | 9.7 | 北 | 0.0 | 41.0 | 25.0 | 2018-01-04 | 0.94 | Sunny |
4 | 3.5 | 6.1 | 0.6 | 0.0 | 0.0 | 4.1 | 北北西 | 2.2 | 北 | 7.9 | 北北西 | 0.0 | 60.0 | 38.0 | 2018-01-05 | 0.00 | Cloudy |
これで日単位のコインランドリー稼働回数データと気象データが準備できましたので、dateカラムを軸に両者を結合します。
combo_data = pd.merge(laundry_daily, obs_daily, on="date", how="left")
combo_data.head()
date | コース1 | コース2 | コース3 | コース4 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | 日照時間rate | 天気 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018-09-01 | 7 | 10 | 3 | 24 | 26.5 | 30.7 | 22.8 | 1.0 | 1.7 | 7.6 | 南 | 4.1 | 南 | 10.6 | 北東 | 1.0 | 86.0 | 66.0 | 0.144068 | Rainy |
1 | 2018-09-02 | 14 | 14 | 3 | 77 | 23.1 | 24.5 | 21.3 | 5.5 | 0.0 | 4.3 | 北 | 2.4 | 北 | 7.5 | 北 | 4.0 | 95.0 | 87.0 | 0.000000 | Rainy |
2 | 2018-09-03 | 6 | 2 | 1 | 79 | 24.6 | 27.3 | 22.5 | 2.5 | 1.6 | 3.8 | 北 | 2.0 | 北 | 5.7 | 北 | 1.5 | 95.0 | 83.0 | 0.135593 | Rainy |
3 | 2018-09-04 | 12 | 6 | 2 | 42 | 27.3 | 30.2 | 24.8 | 8.0 | 2.2 | 20.0 | 南 | 10.1 | 南 | 28.2 | 南 | 7.5 | 90.0 | 75.0 | 0.186441 | Rainy |
4 | 2018-09-05 | 11 | 11 | 1 | 30 | 27.3 | 29.9 | 22.9 | 8.5 | 8.7 | 17.6 | 南南西 | 11.7 | 南南西 | 24.6 | 南南西 | 5.0 | 75.0 | 56.0 | 0.737288 | Rainy |
これで天気とコースごとの分析を行う準備ができました。
ではさっそく天気とコースごとにクロス集計し、稼働回数の平均値・標準偏差・中央値を求めてみましょう。
mean
関数(std
関数・median
関数)を適用loc[wx_category]
の部分)np.around
関数)# 平均値
course = ["コース1","コース2","コース3","コース4"]
wx_category = ["Sunny","Cloudy","Rainy"]
np.around( combo_data.groupby("天気")[course].mean().loc[wx_category] )
コース1 | コース2 | コース3 | コース4 | |
---|---|---|---|---|
天気 | ||||
Sunny | 12.0 | 6.0 | 4.0 | 27.0 |
Cloudy | 12.0 | 6.0 | 3.0 | 43.0 |
Rainy | 10.0 | 6.0 | 3.0 | 47.0 |
# 標準偏差
np.around( combo_data.groupby("天気")[course].std().loc[wx_category] )
コース1 | コース2 | コース3 | コース4 | |
---|---|---|---|---|
天気 | ||||
Sunny | 6.0 | 3.0 | 2.0 | 11.0 |
Cloudy | 6.0 | 3.0 | 2.0 | 18.0 |
Rainy | 5.0 | 3.0 | 2.0 | 25.0 |
# 中央値
np.around( combo_data.groupby("天気")[course].median().loc[wx_category] )
コース1 | コース2 | コース3 | コース4 | |
---|---|---|---|---|
天気 | ||||
Sunny | 11 | 6 | 3 | 24 |
Cloudy | 11 | 5 | 3 | 38 |
Rainy | 9 | 5 | 2 | 41 |
天気が悪くなるほど、乾燥コースの稼働回数が増えています。
その他のコースはあまり天気で差がありませんが、あえていうと洗濯乾燥(標準)コースは雨の日に少しだけ少ないようにも見えます。
また稼働回数が多いほどばらつきも大きくなっているようです。この様子を箱ひげ図で可視化してみます。
#---箱ひげ図を描画する関数
# 引数
# - data : 稼働回数と気象データ結合後のデータフレーム
# - param : x軸に使う気象要素
# - order : x軸の気象要素の順番
def box_plot(data, param, order=wx_category):
fig = plt.figure(figsize=(20,12))
for i, c in enumerate(course):
ax = fig.add_subplot(2,2,i+1)
ax = sns.boxplot(x=param, y=c, data=data, order=order)
ax.set_title(c, fontsize=15)
ax.set_xlabel(param, fontsize=15)
ax.set_ylabel("稼働回数", fontsize=15)
plt.show()
# 箱ひげ図
box_plot(combo_data, "天気")
天気カテゴリーによる分析では、雨天ほど利用者数(稼働回数)が多い傾向がありました。
そうすると日照時間と利用者数には逆相関があるかもしれません(日照時間が増えるほど稼働回数は減る)。
そこで、今度は連続値の変数である日照時間と利用者数との相関係数を計算してみましょう。Pandasのcorr
関数を使います。
なお、これより先は稼働回数が多い乾燥コースと、洗濯乾燥(標準)コースに絞って見ていきます。
course = ["コース1","コース4"]
combo_data[["日照時間",*course]].corr()
日照時間 | コース1 | コース4 | |
---|---|---|---|
日照時間 | 1.000000 | 0.093591 | -0.569132 |
コース1 | 0.093591 | 1.000000 | 0.309644 |
コース4 | -0.569132 | 0.309644 | 1.000000 |
日照時間とコース1はほとんど相関がありませんが、コース4は負の相関係数が0.5以上あり、それなりに相関がありそうです。
次に、日照時間と稼働回数はどちらも連続値のため、散布図を描いてばらつき具合を可視化してみましょう。
#---散布図を描画する関数
# 引数
# - data : 稼働回数と気象データ結合後のデータフレーム
# - param : x軸に使う気象要素
# - order : x軸の気象要素の順番
def scatter_plot(data, param):
fig = plt.figure(figsize=(20,12))
for i, c in enumerate(course):
ax = fig.add_subplot(2,2,i+1)
ax.scatter(data[param], combo_data[c], alpha=0.5)
ax.set_title(c, fontsize=15)
ax.set_xlabel(param, fontsize=15)
ax.set_ylabel("稼働回数", fontsize=15)
plt.show()
scatter_plot(combo_data, "日照時間")
→ 天気カテゴリーで得られた結果と同様の傾向が見られました。
今日の天気だけでなく、昨日・一昨日の天気も利用者数に影響を与えそうです。
例えば雨の日が数日連続して続くと、自宅で干せないため、コインランドリーで洗濯や乾燥をする人が増える、という仮説が立てられます。
そこで、当日を含めて雨が何日続いているか、すなわち連続降水日数をもとに分析を行います。
連続降水日数の求め方は、天気カラムの値をshift関数で1日ずつずらしたカラムを作り、何日連続でRainyの天気が続いているかを逐次計算していく方法をとっています。
#---連続降水日数を求める関数
# 引数
# - data : 天気カラムを含むデータフレーム
# - max_days : 連続降水日数として想定する最大の日数
def Continuous_Rainy_Days(data, max_days=30):
df = data[["天気"]].copy()
df["Rainy_Days"] = 0
wx_minus_days = []
for d in range(max_days):
df["wx_minus_" + str(d)] = df["天気"].shift(d).apply(lambda x: 1 if x == "Rainy" else 0)
wx_minus_days.append("wx_minus_" + str(d))
df["sum"] = df[wx_minus_days].sum(axis=1)
df["Rainy_Days"] = df["Rainy_Days"].mask(df["sum"] == d + 1, d + 1)
return df["Rainy_Days"].values
obs_daily["連続降水日数"] = Continuous_Rainy_Days(obs_daily)
obs_daily["連続降水日数"].value_counts()
0 660 1 168 2 88 3 38 4 19 5 8 6 5 8 2 7 2 9 2 10 2 11 2 18 1 12 1 13 1 14 1 15 1 16 1 17 1 19 1 Name: 連続降水日数, dtype: int64
どうも梅雨の時期に10日以上連続して雨天が続いたケースがあるようです。
今回はサンプル数のバランスもふまえて、1日連続と2日連続を1つにまとめます。値としては2日連続の場合に1
を代入します。
また3日以上連続も1つにまとめます。値としては4日以上連続の場合に3
を代入します。
これで連続降水日数を、0日連続|1〜2日連続|3日以上連続、の3つのカテゴリーに分けたことになります。
obs_daily["連続降水日数"] = obs_daily["連続降水日数"].mask(obs_daily["連続降水日数"]==2, 1)
obs_daily["連続降水日数"] = obs_daily["連続降水日数"].mask(obs_daily["連続降水日数"]>3, 3)
obs_daily["連続降水日数"].value_counts()
0 660 1 256 3 88 Name: 連続降水日数, dtype: int64
連続降水日数が求まりましたので、改めてコインランドリー稼働回数データと気象データを結合します。
combo_data = pd.merge(laundry_daily, obs_daily, on="date", how="left")
combo_data.head()
date | コース1 | コース2 | コース3 | コース4 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | 日照時間rate | 天気 | 連続降水日数 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018-09-01 | 7 | 10 | 3 | 24 | 26.5 | 30.7 | 22.8 | 1.0 | 1.7 | 7.6 | 南 | 4.1 | 南 | 10.6 | 北東 | 1.0 | 86.0 | 66.0 | 0.144068 | Rainy | 1 |
1 | 2018-09-02 | 14 | 14 | 3 | 77 | 23.1 | 24.5 | 21.3 | 5.5 | 0.0 | 4.3 | 北 | 2.4 | 北 | 7.5 | 北 | 4.0 | 95.0 | 87.0 | 0.000000 | Rainy | 1 |
2 | 2018-09-03 | 6 | 2 | 1 | 79 | 24.6 | 27.3 | 22.5 | 2.5 | 1.6 | 3.8 | 北 | 2.0 | 北 | 5.7 | 北 | 1.5 | 95.0 | 83.0 | 0.135593 | Rainy | 3 |
3 | 2018-09-04 | 12 | 6 | 2 | 42 | 27.3 | 30.2 | 24.8 | 8.0 | 2.2 | 20.0 | 南 | 10.1 | 南 | 28.2 | 南 | 7.5 | 90.0 | 75.0 | 0.186441 | Rainy | 3 |
4 | 2018-09-05 | 11 | 11 | 1 | 30 | 27.3 | 29.9 | 22.9 | 8.5 | 8.7 | 17.6 | 南南西 | 11.7 | 南南西 | 24.6 | 南南西 | 5.0 | 75.0 | 56.0 | 0.737288 | Rainy | 3 |
では連続降水日数とコースごとにクロス集計しましょう。
統計量としては中央値を使い、またコースも洗濯乾燥コース(標準)と乾燥コースに絞って見ていきます。
# 中央値
course = ["コース1","コース4"]
rainy_days = [0,1,3]
np.around( combo_data.groupby("連続降水日数")[course].median() )
コース1 | コース4 | |
---|---|---|
連続降水日数 | ||
0 | 11 | 26 |
1 | 9 | 39 |
3 | 10 | 53 |
# 箱ひげ図
box_plot(combo_data, "連続降水日数", rainy_days)
降水量が増えると利用者数にどのような影響があるでしょうか?
乾燥コースの利用者数が増えそうですが、逆にしっかり雨が降っていると、そもそも外出を避ける傾向もあるかもしれません。
そこで、まず降水量とコースごとの利用者数との相関係数を計算してみます。
combo_data[["降水量",*course]].corr()
降水量 | コース1 | コース4 | |
---|---|---|---|
降水量 | 1.000000 | -0.121698 | 0.232771 |
コース1 | -0.121698 | 1.000000 | 0.309644 |
コース4 | 0.232771 | 0.309644 | 1.000000 |
降水量と乾燥コースの相関係数はプラスの値になりましたが、相関係数自体は大きくありません。
次に、降水量と稼働回数はどちらも連続値のため、散布図を描いてばらつき具合を可視化してみましょう。
scatter_plot(combo_data, "降水量")
湿度が上がると利用者数にどのような影響があるでしょうか?
湿度が高いと洗濯物が乾きにくいため乾燥コースの利用者が増えそうです。
(ただし、単純に雨天の時は湿度が高いので、雨の影響を間接的に見ていることになるかもしれません。)
まず相対湿度とコースごとの利用者数との相関係数を計算してみます。
combo_data[["平均湿度","最小相対湿度",*course]].corr()
平均湿度 | 最小相対湿度 | コース1 | コース4 | |
---|---|---|---|---|
平均湿度 | 1.000000 | 0.919387 | -0.004535 | 0.510514 |
最小相対湿度 | 0.919387 | 1.000000 | -0.009615 | 0.567948 |
コース1 | -0.004535 | -0.009615 | 1.000000 | 0.309644 |
コース4 | 0.510514 | 0.567948 | 0.309644 | 1.000000 |
湿度と洗濯乾燥コースには相関がなさそうですが、乾燥コースは相関係数0.5以上あり、それなりに相関がありそうです。
また平均湿度より最小相対湿度の方が少しだけ相関係数が大きくなっています。
そこで最小相対湿度について、散布図を描いてばらつき具合を可視化してみます。
scatter_plot(combo_data, "最小相対湿度")
→ ただし、雨天の時に湿度が上がる様子が現れているだけかもしれないので、解釈には注意が必要です。
先ほどは空気の水分量を相対湿度で扱いましたが、空調の分野では絶対湿度、すなわち水蒸気量も扱うようです。
確かに同じ相対湿度でも、気温によって空気中の水蒸気量は異なります。
洗濯物を干す場合でも、空気があとどれくらいの水蒸気量を含むことができるか、という指標が洗濯物の乾き具合に影響しそうです。
つまり、空気が追加で含みうる水蒸気量が小さいと、洗濯物が乾きにくくて乾燥コースの利用者が増える、との仮説を立てられます。
そこで、平均気温と平均湿度から飽和水蒸気量と水蒸気量を計算し、両者の差分を取ることで「追加で含みうる水蒸気量」を求めてみます。
計算にはwxparamsを利用します。
また「追加で含みうる水蒸気量」では長いので「可蒸発量」という造語でカラム名をつけることとします。
import wxparams as wx
# 可蒸発量を計算
obs_daily["露点温度"] = wx.RH_to_Td(obs_daily["平均気温"], obs_daily["平均湿度"])
obs_daily["水蒸気量"] = wx.Absolute_Humidity(obs_daily["平均気温"], obs_daily["露点温度"])
obs_daily["飽和水蒸気量"] = wx.Absolute_Humidity(obs_daily["平均気温"], obs_daily["平均気温"])
obs_daily["可蒸発量"] = obs_daily["飽和水蒸気量"] - obs_daily["水蒸気量"]
# 不要なカラムを落としてマージ
obs_daily.drop(columns=["露点温度","飽和水蒸気量"], inplace=True)
combo_data = pd.merge(laundry_daily, obs_daily, on="date", how="left")
combo_data.head()
date | コース1 | コース2 | コース3 | コース4 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | 日照時間rate | 天気 | 連続降水日数 | 水蒸気量 | 可蒸発量 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018-09-01 | 7 | 10 | 3 | 24 | 26.5 | 30.7 | 22.8 | 1.0 | 1.7 | 7.6 | 南 | 4.1 | 南 | 10.6 | 北東 | 1.0 | 86.0 | 66.0 | 0.144068 | Rainy | 1 | 21.530870 | 3.505025 |
1 | 2018-09-02 | 14 | 14 | 3 | 77 | 23.1 | 24.5 | 21.3 | 5.5 | 0.0 | 4.3 | 北 | 2.4 | 北 | 7.5 | 北 | 4.0 | 95.0 | 87.0 | 0.000000 | Rainy | 1 | 19.632719 | 1.033301 |
2 | 2018-09-03 | 6 | 2 | 1 | 79 | 24.6 | 27.3 | 22.5 | 2.5 | 1.6 | 3.8 | 北 | 2.0 | 北 | 5.7 | 北 | 1.5 | 95.0 | 83.0 | 0.135593 | Rainy | 3 | 21.379729 | 1.125249 |
3 | 2018-09-04 | 12 | 6 | 2 | 42 | 27.3 | 30.2 | 24.8 | 8.0 | 2.2 | 20.0 | 南 | 10.1 | 南 | 28.2 | 南 | 7.5 | 90.0 | 75.0 | 0.186441 | Rainy | 3 | 23.555547 | 2.617283 |
4 | 2018-09-05 | 11 | 11 | 1 | 30 | 27.3 | 29.9 | 22.9 | 8.5 | 8.7 | 17.6 | 南南西 | 11.7 | 南南西 | 24.6 | 南南西 | 5.0 | 75.0 | 56.0 | 0.737288 | Rainy | 3 | 19.629623 | 6.543208 |
combo_data[["可蒸発量","水蒸気量",*course]].corr()
可蒸発量 | 水蒸気量 | コース1 | コース4 | |
---|---|---|---|---|
可蒸発量 | 1.000000 | 0.033651 | 0.078660 | -0.591520 |
水蒸気量 | 0.033651 | 1.000000 | 0.083705 | 0.211929 |
コース1 | 0.078660 | 0.083705 | 1.000000 | 0.309644 |
コース4 | -0.591520 | 0.211929 | 0.309644 | 1.000000 |
想定した通り、乾燥コースの稼働回数と可蒸発量には、それなりの相関がありそうです(負の相関)。
また洗濯乾燥(標準)コースと稼働回数には相関があまりなさそうです。
では散布図も描いて、ばらつき具合を可視化してみましょう。
scatter_plot(combo_data, "可蒸発量")
先ほどは単純に日積算の降水量と日照時間で天気を推定しました。
しかしコインランドリー利用者、すなわち洗濯する人の行動や気持ちを想像すると、午前中の気象条件の方が影響が大きいかもしれません。
そこで今度は時別の気象データから午前中の天気を推測してみます。
まず時別のアメダス観測値を読み込みます。
時別のアメダス観測値は複数のファイルに分かれているため、まずglob関数でファイル名のリストを取得します。
そして日別アメダス観測値を読み込む際に作ったRead_AMeDAS_CSV関数で読み込みを行います。
#---アメダス江戸川臨海の観測値
all_columns3 = ["年月日時","気温","気温_q","気温_i","降水量","降水量_q","降水量_i",
"日照時間","日照時間_q","日照時間_i","風速","風速_q","風向","風向_q","風速_i"]
select_columns3 = ["年月日時","気温","降水量","日照時間","風速","風向"]
# 全ファイルをまとめて読み込み
dirname3 = "weather"
filenames3 = "edogawarinkai_hourly_*"
files3 = glob.glob( os.path.join(dirname3,filenames3) )
print("読み込むファイル")
obs3 = pd.DataFrame()
for file in files3:
print(file)
tmp = Read_AMeDAS_CSV(file, all_columns3, select_columns3)
obs3 = pd.concat([obs3,tmp])
# 観測時刻順にソート
obs3["年月日時"] = pd.to_datetime(obs3["年月日時"])
obs3 = obs3.sort_values("年月日時").reset_index(drop=True)
読み込むファイル weather/edogawarinkai_hourly_20200701_20200930.csv weather/edogawarinkai_hourly_20180701_20180930.csv weather/edogawarinkai_hourly_20191001_20191231.csv weather/edogawarinkai_hourly_20181001_20181231.csv weather/edogawarinkai_hourly_20190701_20190930.csv weather/edogawarinkai_hourly_20190401_20190630.csv weather/edogawarinkai_hourly_20190101_20190331.csv weather/edogawarinkai_hourly_20180101_20180331.csv weather/edogawarinkai_hourly_20180401_20180630.csv weather/edogawarinkai_hourly_20200101_20200331.csv weather/edogawarinkai_hourly_20200401_20200630.csv
#---東京管区気象台の観測値
all_columns4 = ["年月日時","相対湿度","相対湿度_q","相対湿度_i"]
select_columns4 = ["年月日時","相対湿度"]
# 全ファイルをまとめて読み込み
dirname4 = "weather"
filenames4 = "tokyo_hourly_*"
files4 = glob.glob( os.path.join(dirname4,filenames4) )
print("読み込むファイル")
obs4 = pd.DataFrame()
for file in files4:
print(file)
tmp = Read_AMeDAS_CSV(file, all_columns4, select_columns4, start_index=5)
obs4 = pd.concat([obs4,tmp])
# 観測時刻順にソート
obs4["年月日時"] = pd.to_datetime(obs4["年月日時"])
obs4 = obs4.sort_values("年月日時").reset_index(drop=True)
読み込むファイル weather/tokyo_hourly_20190101_20191231.csv weather/tokyo_hourly_20180101_20181231.csv weather/tokyo_hourly_20200101_20200930.csv
#---合成
obs_hourly = pd.merge(obs3, obs4, on="年月日時", how="left")
obs_hourly.rename(columns={"年月日時":"datetime"}, inplace=True)
obs_hourly.head()
datetime | 気温 | 降水量 | 日照時間 | 風速 | 風向 | 相対湿度 | |
---|---|---|---|---|---|---|---|
0 | 2018-01-01 01:00:00 | 3.2 | 0.0 | NaN | 1.2 | 北北東 | 82.0 |
1 | 2018-01-01 02:00:00 | 2.1 | 0.0 | NaN | 2.0 | 北北東 | 83.0 |
2 | 2018-01-01 03:00:00 | 2.2 | 0.0 | NaN | 1.8 | 北北東 | 80.0 |
3 | 2018-01-01 04:00:00 | 1.8 | 0.0 | NaN | 1.5 | 北北東 | 85.0 |
4 | 2018-01-01 05:00:00 | 1.5 | 0.0 | NaN | 2.0 | 北 | 80.0 |
これで時別のアメダス観測値を読み込むことができました。
次は午前中の天気の推定です。
日天気の推定と同様に、午前中の日照時間rateと降水量を使って天気を推定します。
# 日付情報の処理
obs_hourly["date"] = obs_hourly["datetime"].dt.strftime("%Y-%m-%d")
obs_hourly["month"] = obs_hourly["datetime"].dt.month
obs_hourly["hour"] = obs_hourly["datetime"].dt.hour
#---午前中の天気を判定する
# データを6時〜12時に絞り、日付で集計して降水量と日照時間の合計を計算
obs_am = obs_hourly.query('6 < hour <= 12').groupby("date").agg({"日照時間":np.sum,
"降水量":np.sum,
"month":np.mean})
# 月ごとに日照時間の最大値を求め、日照時間rateを計算
obs_am = obs_am.rename(columns={"日照時間":"日照時間am","降水量":"降水量am"}).reset_index()
max_sun_am = obs_am.groupby("month")[["日照時間am"]].max().reset_index()
max_sun_am.columns = ["month","日照時間max_am"]
obs_am = pd.merge(obs_am, max_sun_am, on="month", how="left")
obs_am["日照時間rate_am"] = obs_am["日照時間am"] / obs_am["日照時間max_am"]
# 日照時間rate_amと降水量amから天気を推定
obs_am["天気am"] = obs_am[["日照時間rate_am","降水量am"]].apply(lambda x: Estimate_Weather(x[0],x[1]), axis=1)
obs_am = obs_am[["date","天気am","降水量am"]]
obs_am.head()
date | 天気am | 降水量am | |
---|---|---|---|
0 | 2018-01-01 | Sunny | 0.0 |
1 | 2018-01-02 | Sunny | 0.0 |
2 | 2018-01-03 | Sunny | 0.0 |
3 | 2018-01-04 | Sunny | 0.0 |
4 | 2018-01-05 | Cloudy | 0.0 |
午前中の天気が求まりましたので、変数obs_dailyに結合します。
また稼働回数データにも結合します。
# obs_dailyに午前中の天気を結合
obs_am["date"] = pd.to_datetime(obs_am["date"])
obs_daily = pd.merge(obs_daily, obs_am, on="date", how="left")
# 稼働回数と気象データを結合
combo_data = pd.merge(laundry_daily, obs_daily, on="date", how="left")
combo_data.head()
date | コース1 | コース2 | コース3 | コース4 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | 日照時間rate | 天気 | 連続降水日数 | 水蒸気量 | 可蒸発量 | 天気am | 降水量am | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018-09-01 | 7 | 10 | 3 | 24 | 26.5 | 30.7 | 22.8 | 1.0 | 1.7 | 7.6 | 南 | 4.1 | 南 | 10.6 | 北東 | 1.0 | 86.0 | 66.0 | 0.144068 | Rainy | 1 | 21.530870 | 3.505025 | Cloudy | 0.0 |
1 | 2018-09-02 | 14 | 14 | 3 | 77 | 23.1 | 24.5 | 21.3 | 5.5 | 0.0 | 4.3 | 北 | 2.4 | 北 | 7.5 | 北 | 4.0 | 95.0 | 87.0 | 0.000000 | Rainy | 1 | 19.632719 | 1.033301 | Rainy | 5.5 |
2 | 2018-09-03 | 6 | 2 | 1 | 79 | 24.6 | 27.3 | 22.5 | 2.5 | 1.6 | 3.8 | 北 | 2.0 | 北 | 5.7 | 北 | 1.5 | 95.0 | 83.0 | 0.135593 | Rainy | 3 | 21.379729 | 1.125249 | Rainy | 2.0 |
3 | 2018-09-04 | 12 | 6 | 2 | 42 | 27.3 | 30.2 | 24.8 | 8.0 | 2.2 | 20.0 | 南 | 10.1 | 南 | 28.2 | 南 | 7.5 | 90.0 | 75.0 | 0.186441 | Rainy | 3 | 23.555547 | 2.617283 | Rainy | 7.5 |
4 | 2018-09-05 | 11 | 11 | 1 | 30 | 27.3 | 29.9 | 22.9 | 8.5 | 8.7 | 17.6 | 南南西 | 11.7 | 南南西 | 24.6 | 南南西 | 5.0 | 75.0 | 56.0 | 0.737288 | Rainy | 3 | 19.629623 | 6.543208 | Rainy | 0.5 |
では、これまでと同様に午前中の天気とコースごとにクロス集計しましょう。
統計量としては中央値を使い、またコースも洗濯乾燥コース(標準)と乾燥コースに絞って見ていきます。
# 中央値
np.around( combo_data.groupby("天気am")[course].median().loc[wx_category] )
コース1 | コース4 | |
---|---|---|
天気am | ||
Sunny | 10.0 | 25.0 |
Cloudy | 10.0 | 40.0 |
Rainy | 9.0 | 49.0 |
# 箱ひげ図
box_plot(combo_data, "天気am", wx_category)
朝は晴れていたので洗濯物を干して出かけたが、急なゲリラ豪雨(夕立)で洗濯物が全滅…
なんてことがあると、乾燥のためコインランドリーを利用する、という人が増えそうです。
そこで今度は時別の気象データから夕立の有無を判別し、利用者数の違いを分析します。
夕立を次のように定義します。
#---夕立の有無を判別
# 7〜9月の午後の積算降水量
obs_pm = obs_hourly.query('7 <= month <= 9').query('12 < hour <= 19').groupby("date")[["降水量"]].sum()
obs_pm = obs_pm.rename(columns={"降水量":"降水量pm"}).reset_index()
# 変数obs_dailyに降水量pmを結合
obs_pm["date"] = pd.to_datetime(obs_pm["date"])
obs_daily = pd.merge(obs_daily, obs_pm, on="date", how="left")
# 夕立カラムを作成
obs_daily["夕立"] = 0
obs_daily["夕立"] = obs_daily["夕立"].mask((obs_daily["天気am"] == "Sunny") & (obs_daily["降水量pm"] > 0.0), 1)
obs_daily.drop(columns=["降水量pm"], inplace=True)
len( obs_daily.query('夕立 == 1') )
12
夕立の定義を満たす日は、12日しかないようです。
ちょっと少ない気もしますが、今までと同様に夕立の有無で集計したり、箱ひげ図を描いてみます。
以下の条件の日を抽出して比較を行います。
# 稼働回数と気象データを結合
combo_data = pd.merge(laundry_daily, obs_daily, on="date", how="left")
combo_data["month"] = combo_data["date"].dt.month
combo_data = combo_data.query('7 <= month <= 9').query('天気am == "Sunny"')
combo_data.head()
date | コース1 | コース2 | コース3 | コース4 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 最大風速の風向 | 平均風速 | 最多風向 | 最大瞬間風速 | 最大瞬間風速の風向 | 最大1時間降水量 | 平均湿度 | 最小相対湿度 | 日照時間rate | 天気 | 連続降水日数 | 水蒸気量 | 可蒸発量 | 天気am | 降水量am | 夕立 | month | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
5 | 2018-09-06 | 6 | 4 | 2 | 22 | 27.9 | 31.9 | 25.9 | 0.0 | 7.0 | 10.7 | 南西 | 7.7 | 南南西 | 13.5 | 南南西 | 0.0 | 67.0 | 49.0 | 0.593220 | Sunny | 0 | 18.126505 | 8.927980 | Sunny | 0.0 | 0 | 9 |
7 | 2018-09-08 | 9 | 12 | 7 | 35 | 27.6 | 30.1 | 26.4 | 0.0 | 10.7 | 13.0 | 南南西 | 9.5 | 南南西 | 16.5 | 南南西 | 0.0 | 74.0 | 57.0 | 0.906780 | Sunny | 0 | 19.691778 | 6.918733 | Sunny | 0.0 | 0 | 9 |
8 | 2018-09-09 | 10 | 10 | 3 | 21 | 27.5 | 29.8 | 26.1 | 0.0 | 10.9 | 11.0 | 南 | 8.4 | 南 | 14.2 | 南 | 0.0 | 76.0 | 52.0 | 0.923729 | Sunny | 0 | 20.112580 | 6.351341 | Sunny | 0.0 | 0 | 9 |
9 | 2018-09-10 | 6 | 2 | 1 | 27 | 25.6 | 29.5 | 21.6 | 37.0 | 3.7 | 9.0 | 南 | 5.1 | 南 | 11.4 | 南 | 20.0 | 89.0 | 66.0 | 0.313559 | Rainy | 1 | 21.189125 | 2.618881 | Sunny | 0.0 | 1 | 9 |
11 | 2018-09-12 | 8 | 9 | 1 | 30 | 21.4 | 24.9 | 18.7 | 0.0 | 3.9 | 6.0 | 北北東 | 3.6 | 北東 | 10.1 | 北北東 | 0.0 | 70.0 | 54.0 | 0.330508 | Sunny | 0 | 13.118176 | 5.622076 | Sunny | 0.0 | 0 | 9 |
# 中央値
shower = [0,1]
np.around( combo_data.groupby("夕立")[course].median() )
コース1 | コース4 | |
---|---|---|
夕立 | ||
0 | 12.0 | 24.0 |
1 | 14.0 | 29.0 |
# 箱ひげ図
box_plot(combo_data, "夕立", shower)
ここまでの分析結果を簡単にまとめます。
#---改めて稼働回数と気象データの結合
combo_data = pd.merge(laundry_daily, obs_daily, on="date", how="left")
#---暦情報のカラムを追加
# 月の情報
# → 7月の稼働回数が顕著に多いが、雨の影響と考えられ天気で説明できるため、月は説明変数に含めない
# combo_data["月"] = combo_data["date"].dt.month
# 曜日の情報
# 休日(=1)と平日(=0)に分けた変数とする
# → 今回は祝日は考慮しないが、祝日情報を入手して休日に含めても良いでしょう
combo_data["曜日"] = combo_data["date"].dt.weekday
combo_data["休日"] = combo_data["曜日"].apply(lambda x: 0 if x < 5 else 1)
# weekday = ["月","火","水","木","金","土","日"]
# combo_data["曜日"] = combo_data["曜日"].apply(lambda x: weekday[x])
#---不要なカラムを削除
# モデリングは代表してコース1とコース4で行う
# 直接使用しない気象要素も削除
drop_columns = ["コース2","コース3","曜日",
"最大風速の風向","最多風向","最大瞬間風速","最大瞬間風速の風向","最大1時間降水量"]
combo_data.drop(columns=drop_columns, inplace=True)
#---過去の実績を追加
combo_data["コース1前日"] = combo_data["コース1"].shift(1)
combo_data["コース1前々日"] = combo_data["コース1"].shift(2)
combo_data["コース4前日"] = combo_data["コース4"].shift(1)
combo_data["コース4前々日"] = combo_data["コース4"].shift(2)
#---欠損値の処理
# 欠損値を含む行を削除
combo_data = combo_data.dropna().reset_index(drop=True)
combo_data.head()
date | コース1 | コース4 | 平均気温 | 最高気温 | 最低気温 | 降水量 | 日照時間 | 最大風速 | 平均風速 | 平均湿度 | 最小相対湿度 | 日照時間rate | 天気 | 連続降水日数 | 水蒸気量 | 可蒸発量 | 天気am | 降水量am | 夕立 | 休日 | コース1前日 | コース1前々日 | コース4前日 | コース4前々日 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018-09-03 | 6 | 79 | 24.6 | 27.3 | 22.5 | 2.5 | 1.6 | 3.8 | 2.0 | 95.0 | 83.0 | 0.135593 | Rainy | 3 | 21.379729 | 1.125249 | Rainy | 2.0 | 0 | 0 | 14.0 | 7.0 | 77.0 | 24.0 |
1 | 2018-09-04 | 12 | 42 | 27.3 | 30.2 | 24.8 | 8.0 | 2.2 | 20.0 | 10.1 | 90.0 | 75.0 | 0.186441 | Rainy | 3 | 23.555547 | 2.617283 | Rainy | 7.5 | 0 | 0 | 6.0 | 14.0 | 79.0 | 77.0 |
2 | 2018-09-05 | 11 | 30 | 27.3 | 29.9 | 22.9 | 8.5 | 8.7 | 17.6 | 11.7 | 75.0 | 56.0 | 0.737288 | Rainy | 3 | 19.629623 | 6.543208 | Rainy | 0.5 | 0 | 0 | 12.0 | 6.0 | 42.0 | 79.0 |
3 | 2018-09-06 | 6 | 22 | 27.9 | 31.9 | 25.9 | 0.0 | 7.0 | 10.7 | 7.7 | 67.0 | 49.0 | 0.593220 | Sunny | 0 | 18.126505 | 8.927980 | Sunny | 0.0 | 0 | 0 | 11.0 | 12.0 | 30.0 | 42.0 |
4 | 2018-09-07 | 5 | 19 | 27.1 | 29.7 | 25.8 | 0.0 | 2.1 | 13.4 | 10.3 | 75.0 | 65.0 | 0.177966 | Cloudy | 0 | 19.413378 | 6.471126 | Cloudy | 0.0 | 0 | 0 | 6.0 | 11.0 | 22.0 | 30.0 |
#---CSV出力
# 出力したCSVをExcelで開いても文字化けしないよう、Shift-JISで出力します
combo_data.to_csv("data_forR.csv", index=False, encoding="shift-jis")
# 事前に配布したCSVファイルから前処理済みデータを読み込む場合は、以下のコメントアウトを外して実行して下さい。
# combo_data = pd.read_csv("prepared_data_forR.csv", encoding="shift-jis")
# combo_data["date"] = pd.to_datetime(combo_data["date"])
# Rの変数選択で最終的に残った特徴量に絞る
use_cols = ["date","コース4","最高気温","日照時間","最大風速","最小相対湿度",
"連続降水日数","天気am","降水量am","休日","コース4前日"]
combo_data = combo_data[use_cols]
Pythonの場合、カテゴリー変数はダミー変数化するなどの前処理が必要になります。
ここでは連続降水日数と月もカテゴリー変数として扱うため、型を文字列に変換します。
そしてカテゴリー変数をダミー変数化する処理をします。
# 文字列型へ変換
to_category = ["連続降水日数"]
for name in to_category:
combo_data[name] = combo_data[name].astype(str)
# ダミー変数化
combo_data = pd.get_dummies(combo_data)
combo_data.head()
date | コース4 | 最高気温 | 日照時間 | 最大風速 | 最小相対湿度 | 降水量am | 休日 | コース4前日 | 連続降水日数_0 | 連続降水日数_1 | 連続降水日数_3 | 天気am_Cloudy | 天気am_Rainy | 天気am_Sunny | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018-09-03 | 79 | 27.3 | 1.6 | 3.8 | 83.0 | 2.0 | 0 | 77.0 | 0 | 0 | 1 | 0 | 1 | 0 |
1 | 2018-09-04 | 42 | 30.2 | 2.2 | 20.0 | 75.0 | 7.5 | 0 | 79.0 | 0 | 0 | 1 | 0 | 1 | 0 |
2 | 2018-09-05 | 30 | 29.9 | 8.7 | 17.6 | 56.0 | 0.5 | 0 | 42.0 | 0 | 0 | 1 | 0 | 1 | 0 |
3 | 2018-09-06 | 22 | 31.9 | 7.0 | 10.7 | 49.0 | 0.0 | 0 | 30.0 | 1 | 0 | 0 | 0 | 0 | 1 |
4 | 2018-09-07 | 19 | 29.7 | 2.1 | 13.4 | 65.0 | 0.0 | 0 | 22.0 | 1 | 0 | 0 | 1 | 0 | 0 |
ダミー変数化したら、多重共線性を回避するため、各カテゴリー変数とも1種類のカテゴリーを削除する必要があります。
ここでは、Rの結果に倣ってカテゴリーの1種類を削除します。
exclude_cols = ["連続降水日数_0","天気am_Cloudy"]
combo_data.drop(columns=exclude_cols, inplace=True)
combo_data.head()
date | コース4 | 最高気温 | 日照時間 | 最大風速 | 最小相対湿度 | 降水量am | 休日 | コース4前日 | 連続降水日数_1 | 連続降水日数_3 | 天気am_Rainy | 天気am_Sunny | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2018-09-03 | 79 | 27.3 | 1.6 | 3.8 | 83.0 | 2.0 | 0 | 77.0 | 0 | 1 | 1 | 0 |
1 | 2018-09-04 | 42 | 30.2 | 2.2 | 20.0 | 75.0 | 7.5 | 0 | 79.0 | 0 | 1 | 1 | 0 |
2 | 2018-09-05 | 30 | 29.9 | 8.7 | 17.6 | 56.0 | 0.5 | 0 | 42.0 | 0 | 1 | 1 | 0 |
3 | 2018-09-06 | 22 | 31.9 | 7.0 | 10.7 | 49.0 | 0.0 | 0 | 30.0 | 0 | 0 | 0 | 1 |
4 | 2018-09-07 | 19 | 29.7 | 2.1 | 13.4 | 65.0 | 0.0 | 0 | 22.0 | 0 | 0 | 0 | 0 |
機械学習では、データを学習用・検証用・テスト用に分けて予測モデル開発を行います。
ここでは最初の1年分を学習用、残りの1年分を検証用とします。
(またテスト用データについては後述します。)
# 日付情報からtrainとvalidationに分割
train = combo_data.query('date < "2019-09-03"')
val = combo_data.query('"2019-09-03" <= date')
train.shape, val.shape
((362, 13), (362, 13))
ちょうど同じデータ数に分割されました。
次に予測ターゲット(y)と、特徴量(X)に分割します。
X_train = train.drop(columns=["date","コース4"])
y_train = train["コース4"]
X_val = val.drop(columns=["date","コース4"])
y_val = val["コース4"]
データの準備が整いましたので、線形回帰から試してみましょう。
# scikit-learnから線形回帰と評価関数もインポート
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, mean_absolute_error
# モデルの学習
lr = LinearRegression()
lr.fit(X_train, y_train)
# 予測
y_pred_lr = lr.predict(X_val)
# 評価
rmse = np.sqrt(mean_squared_error(y_val, y_pred_lr))
mae = mean_absolute_error(y_val, y_pred_lr)
print("RMSE : {:.3f}, MAE : {:.3f}".format(rmse,mae))
RMSE : 11.765, MAE : 8.840
from sklearn.ensemble import RandomForestRegressor
# モデルの学習
rf = RandomForestRegressor(n_estimators=500, n_jobs=-1, random_state=0)
rf.fit(X_train, y_train)
# 予測
y_pred_rf = rf.predict(X_val)
# 評価
rmse = np.sqrt(mean_squared_error(y_val, y_pred_rf))
mae = mean_absolute_error(y_val, y_pred_rf)
print("RMSE : {:.3f}, MAE : {:.3f}".format(rmse,mae))
RMSE : 12.169, MAE : 8.948
ランダムフォレストはハイパーパラメータの調整で予測精度が変わりますが、今回は木の数(n_estimators)を500本にしてみました。
また学習の再現性を確保するため、randam_stateを固定しています。
結果を見ると、MAE・RMSEともに線形回帰の方が精度が良い結果となっています。
どのような予測が出ていたのか、散布図で可視化してみます。
#-x軸・y軸の値の範囲
axis_range = (0,160)
#-散布図
fig = plt.figure(figsize=(13,6))
# 線形回帰の結果を可視化
ax = fig.add_subplot(1,2,1)
ax.scatter(y_val, y_pred_lr, alpha=0.5)
ax.plot(axis_range, axis_range, color="r")
ax.set_title("線形回帰", fontsize=15)
ax.set_xlabel("実績", fontsize=15)
ax.set_ylabel("予測", fontsize=15)
ax.set_xlim(axis_range)
ax.set_ylim(axis_range)
# ランダムフォレストの結果を可視化
ax = fig.add_subplot(1,2,2)
ax.scatter(y_val, y_pred_rf, alpha=0.5)
ax.plot(axis_range, axis_range, color="r")
ax.set_title("ランダムフォレスト", fontsize=15)
ax.set_xlabel("実績", fontsize=15)
ax.set_ylabel("予測", fontsize=15)
ax.set_xlim(axis_range)
ax.set_ylim(axis_range)
plt.show()
線形回帰に比べると、ランダムフォレストの方が大外しをするケースが多いため、RMSEが悪くなったと考えられます。
またランダムフォレストはアルゴリズムの特性上、実績の下限値付近を予測しにくい傾向が現れています。
今回、ランダムフォレストに入力した特徴量は線形回帰と同じものにしました。
しかし、ランダムフォレストの場合は多重共線性は気にしなくて良いので、入力できる全ての特徴量を利用しても良いのかもしれません。
ぜひ試してみて下さい。
では、各モデルがどのような学習をしたか見てみましょう。
まず線形回帰の回帰係数をデータフレームで書き出します。
# 回帰係数を表示
pd.DataFrame(lr.coef_, index=X_train.columns, columns=["回帰係数"])
回帰係数 | |
---|---|
最高気温 | -0.149203 |
日照時間 | -0.628070 |
最大風速 | -0.316856 |
最小相対湿度 | 0.261814 |
降水量am | -0.656658 |
休日 | 15.292475 |
コース4前日 | 0.366725 |
連続降水日数_1 | 0.547678 |
連続降水日数_3 | 6.362162 |
天気am_Rainy | 9.778133 |
天気am_Sunny | -4.518166 |
Rの線形回帰と全く同じにはなりませんが、同じような傾向(正負の符号)で回帰係数が求まっていると思います。
次にランダムフォレストで重要度を見てみましょう。
# feature importanceを表示
pd.DataFrame(rf.feature_importances_,
index=X_train.columns,
columns=["重要度"]).sort_values("重要度", ascending=False)
重要度 | |
---|---|
最小相対湿度 | 0.430535 |
コース4前日 | 0.252929 |
日照時間 | 0.098470 |
休日 | 0.085882 |
最大風速 | 0.047302 |
最高気温 | 0.046807 |
降水量am | 0.014925 |
天気am_Sunny | 0.009807 |
天気am_Rainy | 0.005055 |
連続降水日数_1 | 0.004269 |
連続降水日数_3 | 0.004018 |
import wxbcgrib as wg
import wxparams as wx
今回は2020年9月6日のコース4の稼働回数を、MSM 2020年9月5日12zベースを使って予測します。
MSM(GRIB2形式)は気象フォルダに入っています。
wxbcgribを使い、ウォッシュプラス富士見1丁目店の場所の気象予測データを取得します。
# GRIB2ファイルへのパスの設定
gpvdir = "weather"
gpvname = "Z__C_RJTD_20200905120000_MSM_GPV_Rjp_Lsurf_FH*_grib2.bin"
gpvpath = os.path.join(gpvdir, gpvname)
# 予測対象日の月と曜日を設定
# 2020年9月6日(日曜日)
# month = 9
weekday = 6
# 気象要素名のリスト
parameters = [
"PRMSL_meansealevel", # 海面気圧[Pa]
"PRES_surface", # 地上気圧[Pa]
"TMP_1D5maboveground", # 気温[K]
"RH_1D5maboveground", # 相対湿度[%]
"UGRD_10maboveground", # 東西風(U成分)[m/s]
"VGRD_10maboveground", # 南北風(V成分)[m/s]
"LCDC_surface", # 下層雲量[%]
"MCDC_surface", # 中層雲量[%]
"HCDC_surface", # 上層雲量[%]
"TCDC_surface", # 全雲量[%]
"APCP_surface", # 前1時間の積算降水量[mm]
"DSWRF_surface" # 前1時間の平均日射量[W/m^2]
]
# 店舗の緯度経度
latitude = 35.648286
longitude = 139.892668
# 地点データの切り出し
gpvdata = None
for prm in parameters:
var = wg.DS_from_grnc(gpvpath, prm)
if gpvdata is None:
gpvdata = pd.DataFrame(var.time, columns=["validtime"])
gpvdata[prm] = var.ts(latitude, longitude, fig=False, csv=False)
gpvdata.head()
validtime | PRMSL_meansealevel | PRES_surface | TMP_1D5maboveground | RH_1D5maboveground | UGRD_10maboveground | VGRD_10maboveground | LCDC_surface | MCDC_surface | HCDC_surface | TCDC_surface | APCP_surface | DSWRF_surface | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 2020-09-05 21:00:00 | 101281.023438 | 101262.046875 | 299.891693 | 88.879700 | -2.725462 | -2.341960 | 22.448538 | 28.615793 | 1.766652 | 28.615793 | 9.999000e+20 | 9.999000e+20 |
1 | 2020-09-05 22:00:00 | 101268.578125 | 101246.242188 | 299.698242 | 89.153275 | -3.222933 | -2.709334 | 24.325115 | 23.112719 | 0.000000 | 25.771294 | 2.793501e-02 | 0.000000e+00 |
2 | 2020-09-05 23:00:00 | 101284.054688 | 101264.921875 | 299.548584 | 90.293449 | -3.417754 | -2.419678 | 23.765482 | 24.950085 | 4.272910 | 26.934505 | 1.152631e-03 | 0.000000e+00 |
3 | 2020-09-06 00:00:00 | 101279.273438 | 101257.328125 | 299.394226 | 92.024345 | -3.062622 | -2.101271 | 32.290981 | 34.154152 | 2.364340 | 34.898289 | 8.624951e-02 | 0.000000e+00 |
4 | 2020-09-06 01:00:00 | 101246.820312 | 101224.078125 | 299.268341 | 93.517319 | -3.245892 | -2.059918 | 33.272507 | 30.556480 | 17.444317 | 43.338585 | 2.241600e-02 | 0.000000e+00 |
予測対象である9月6日の予測データに絞ります。
gpvdata = gpvdata.query('"2020-09-06 00:00:00" < validtime <= "2020-09-07 00:00:00"')
gpvdata.head()
validtime | PRMSL_meansealevel | PRES_surface | TMP_1D5maboveground | RH_1D5maboveground | UGRD_10maboveground | VGRD_10maboveground | LCDC_surface | MCDC_surface | HCDC_surface | TCDC_surface | APCP_surface | DSWRF_surface | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
4 | 2020-09-06 01:00:00 | 101246.820312 | 101224.078125 | 299.268341 | 93.517319 | -3.245892 | -2.059918 | 33.272507 | 30.556480 | 17.444317 | 43.338585 | 0.022416 | 0.0 |
5 | 2020-09-06 02:00:00 | 101223.796875 | 101205.039062 | 299.187805 | 93.283592 | -3.571578 | -2.123898 | 37.698940 | 26.272205 | 89.147011 | 89.147011 | 0.340161 | 0.0 |
6 | 2020-09-06 03:00:00 | 101212.867188 | 101190.218750 | 299.133636 | 89.409912 | -3.757829 | -2.134167 | 24.039164 | 23.890709 | 24.896355 | 28.097450 | 1.220781 | 0.0 |
7 | 2020-09-06 04:00:00 | 101219.039062 | 101202.484375 | 298.994659 | 86.910484 | -3.525291 | -2.252622 | 24.256243 | 20.730280 | 0.064835 | 24.256243 | 0.716338 | 0.0 |
8 | 2020-09-06 05:00:00 | 101230.304688 | 101211.531250 | 298.892426 | 85.702026 | -3.523757 | -2.427860 | 23.571444 | 16.607367 | 0.000000 | 23.571444 | 0.228985 | 0.0 |
これで9月6日対象の1時間ごとの予測データを取得できましたが、まだ予測はできません。
学習済みモデルを作った時と同じ特徴量(1日ごとの気象要素など)に変換する必要があります。
今回は1日分のテストデータを作るだけなので、少々強引に1つずつ特徴量を計算していきます。
まずテストデータを入れる「箱」を用意します。
# テストデータ
# とりあえず全部ゼロを入れておく
X_test = pd.DataFrame(np.zeros((1,X_train.shape[1])), columns=X_train.columns)
X_test
最高気温 | 日照時間 | 最大風速 | 最小相対湿度 | 降水量am | 休日 | コース4前日 | 連続降水日数_1 | 連続降水日数_3 | 天気am_Rainy | 天気am_Sunny | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
最高気温・最大風速・最小相対湿度は、それぞれ1時間ごとの値の最大値・最小値を計算して求めることができます。
# 最高気温(絶対温度[K]なので摂氏温度[C]に変換)
X_test.loc[0,"最高気温"] = gpvdata["TMP_1D5maboveground"].max() - 273.15
# 最小相対湿度
X_test.loc[0,"最小相対湿度"] = gpvdata["RH_1D5maboveground"].min()
# 最大風速は、wxparamsを使ってUV成分から風速を求め、その最大値を計算
wspd, wdir = wx.UV_to_SpdDir(gpvdata["UGRD_10maboveground"], gpvdata["VGRD_10maboveground"])
X_test.loc[0,"最大風速"] = wspd.max()
日照時間をどのように求めるかは工夫が必要です。
日射量を使う手もあるかもしれませんが、前1時間の平均値から求めるのは困難です。
ここでは全雲量を使い (100 - 全雲量) / 100
を日照時間の推定値として使用することとします。
# 日照時間の観測では、0.1時間ごとに「直達日射量が0.12kW/m^2以上」を基準として観測されている。
# そこで日射量が120W/m^2以上の時刻について、全雲量から日照時間の推定を行う。
X_test.loc[0,"日照時間"] = (100 - gpvdata.query("DSWRF_surface >= 120")["TCDC_surface"]).sum()/100
午前中の天気と降水量です。
まず数値予報は非常に小さい値の降水量を計算する場合があるので、0.1mm未満の微弱な降水量は切り捨ててゼロとします。
その後、午前7時〜午前12時の予測データを用いて、午前中の天気と降水量を計算します。
# 0.1mm未満の降水量はゼロにする
gpvdata["APCP_surface"] = gpvdata["APCP_surface"].apply(lambda x: 0.0 if 0.0 < x < 0.1 else x)
# 午前のデータを取り出す
gpvdata_am = gpvdata.query('"2020-09-06 06:00:00" < validtime <= "2020-09-06 12:00:00"').copy()
# 晴れ・曇りの推定
# 全雲量80%を超えたら曇りとして1、それ以外は晴れとして0とする
gpvdata_am["晴れor曇り"] = gpvdata_am["TCDC_surface"].apply(lambda x: 1 if x > 80 else 0)
# 午前中の降水量
X_test.loc[0,"降水量am"] = gpvdata_am["APCP_surface"].sum()
# 午前中の天気
# 降水量があれば雨とする
# 午前7時〜12時(6時間)の間の、1時間ごとの天気のうち、曇りが8割以下であれば晴れとする
if X_test.loc[0,"降水量am"] > 0.0:
X_test.loc[0,"天気am_Rainy"] = 1
elif gpvdata_am["晴れor曇り"].sum() / len(gpvdata_am["晴れor曇り"]) <= 0.8:
X_test.loc[0,"天気am_Sunny"] = 1
連続降水日数は、前日9月5日までの連続降水日数を使い、9月6日の予想天気が雨ならば+1することで求まります。
上で示した表を見れば、9月5日までの連続降水日数がわかります。
# 9月5日までの連続降水日数
rainy_days = 1
# 9月6日が雨の予想なら+1、雨でなければ0にリセット
if gpvdata["APCP_surface"].sum() > 0.0:
rainy_days += 1
else:
rainy_days = 0
# 連続降水日数:2日→1日、3日以上→3日に変換
if rainy_days == 2:
rainy_days = 1
elif rainy_days > 3:
rainy_days = 3
# 連続降水日数が0日でなければ、その日数のカラムを1にする
if rainy_days > 0:
X_test.loc[0, "連続降水日数_{}".format(rainy_days)] = 1.0
残りは前日9月5日の稼働回数の実績を代入し、9月と日曜日のカラムを1にします。
# 前日の実績
X_test.loc[0,"コース4前日"] = 21
# 休日 or 平日
if weekday < 5:
X_test.loc[0,"休日"] = 0
else:
X_test.loc[0,"休日"] = 1
# 完成したX_testを出力
X_test
最高気温 | 日照時間 | 最大風速 | 最小相対湿度 | 降水量am | 休日 | コース4前日 | 連続降水日数_1 | 連続降水日数_3 | 天気am_Rainy | 天気am_Sunny | |
---|---|---|---|---|---|---|---|---|---|---|---|
0 | 29.07757 | 4.799683 | 4.321568 | 72.233566 | 1.57674 | 1.0 | 21.0 | 1.0 | 0.0 | 1.0 | 0.0 |
やっとテストデータが完成しました!
では2つの学習済みモデルに投入して、予測を出力しましょう。
#--テストデータを学習済みモデルに投入して予測を行う
# 線形回帰
y_pred_lr = lr.predict(X_test)
# ランダムフォレスト
y_pred_rf = rf.predict(X_test)
print("9月6日のコース4実績 : 64")
print("線形回帰の予測 : ", np.round(y_pred_lr[0]))
print("ランダムフォレストの予測 : ", np.round(y_pred_rf[0]))
9月6日のコース4実績 : 64 線形回帰の予測 : 57.0 ランダムフォレストの予測 : 57.0
線形回帰・ランダムフォレストともに、良い予測ができたと言えるのではないでしょうか。
本PBLでは気象観測データを使って学習モデルを作り、気象予測データを入力してターゲットの予測を得ました。
このようなアプローチをPPM (Perfect Prognosis Method)と呼びます。
これに対し学習モデルを作る段階から過去の気象予測データを使い、気象予測データを入力してターゲットの予測を得る方法もあります。
このようなアプローチをMOS (Model Output Statistics)と呼びます。
気象データは同じ気象要素でも、観測データと予測データという性質の異なる2種類のデータが存在します。
観測データと予測データ、PPMとMOSの違いを理解して使い分ける(使いこなす)ことが重要です。
これまで、コインランドリーの利用者数と気象状況との関係性を分析し、様々な知見を得るとともに、予測まで行うことができました。
これらの結果をもとにして、三密回避&機会ロス削減のためのピークシフト施策を考えてみましょう。