# numpyをnpという別名でインポートします．
import numpy as np
# matplotlibをpltという別名でインポートします．
import matplotlib.pyplot as plt
# Seabornをインポートします．
import seaborn as sns
# datetimeは日時データを処理する際に便利なメソッドです．インポートします．
from datetime import date, datetime, timedelta
# pandasをpdという別名でインポートします．
import pandas as pd
# matplotlibで時系列図を作成するときには以下をインポートします
from pandas.plotting import register_matplotlib_converters
# これを登録しておきます．
register_matplotlib_converters()
# sklearn(scikit-learn)は機械学習関連のライブラリーです．インポートします．
from sklearn import linear_model
%precision 3
%matplotlib inline

# 気象庁アメダスの気温の時系列データを読み込んで，
# DataFrameに割り当てる関数
def readamedas(filename,skipline): 
    # ここがポイント！
    # pandasのread_csvというメソッドでcsvファイルを読み込みます．
    # 引数として，
    # [0]入力ファイル名
    # [1]エンコーディング
    # [2]読み飛ばす行数，
    # [3]column名
    # [4]datetime型で読み込むcolumn名
    # [5]indexとするcolumn名
    # を与える．
    df=pd.read_csv(
       filename,
       encoding='Shift_JIS',
       skiprows=skipline,
       names=['date_time','temp','dummy1','dummy2'],
       parse_dates={'datetime':['date_time']}, 
       index_col='datetime'
       )
    return df

# 気象庁AMeDAS（東京）の気温（℃）のcsvファイルの名前を指定します．
filename1='amedas2017.csv'
# csvファイルの上から5行目まではデータではないため呼び飛ばします．
skipline1=5
# 気象庁AMeDAS（東京）の気温（℃）のcsvファイルを読み込んで，
# pandasのDataFrame（tepco）に割り当てる関数を呼び出します．
amedas=readamedas(filename1,skipline1)
# DataFrame(amedas)の中のdummy1とdummy2の列を削除する．
amedas=amedas.drop(['dummy1','dummy2'],axis=1)
amedas

# 東京電力の電力消費量の時系列データを読み込んで，
# DataFrameに割り当てる関数
def readtepco(filename,skipline): 
    # ここがポイント！
    # pandasのread_csvというメソッドでcsvファイルを読み込みます．
    # 引数として，
    # [0]入力ファイル名
    # [1]エンコーディング
    # [2]読み飛ばす行数，
    # [3]column名
    # [4]datetime型で読み込むcolumn名
    # [5]indexとするcolumn名
    # を与える．
    df=pd.read_csv(
       filename,
       encoding='Shift_JIS',
       skiprows=skipline,
       names=['date','time','power'],
       parse_dates={'datetime':['date','time']},
       index_col='datetime'
       )
    return df

# 東京電力の電力消費量（10^6 kW）のcsvファイルの名前を指定します．
filename2='tepco2017.csv'
# csvファイルの上から3行目まではデータではないため呼び飛ばします．
skipline2=3
# 東京電力の電力消費量（10^6 kW）のcsvファイルを読み込んで，
# pandasのDataFrame（tepco）に割り当てる関数を呼び出します．
tepco=readtepco(filename2,skipline2)
tepco

# ここがポイント！
# 2つのDataFrame（amedasとtepco）を結合します．
data=pd.merge(amedas, tepco, how="outer", left_index=True, right_index=True)
data

# ２つの時系列データからNaNを取り除いて，
# DataFrameに割り当てる関数
def dataprocess0(df,name1,name2):
    # ここがポイント！
    # 2つのデータのどちらかにNaNがあれば，その行を削除します．
    df1=df.dropna(how='any')
    # ここがポイント！
    # 処理した2つのname1とname2のデータを取り出し，
    # 2つのデータ間の相関係数を計算します．
    corr=df1[name1].corr(df1[name2])
    # 処理したデータと相関係数を返します．
    return df1,corr

# データ前処理0：　dataprocess0関数を呼び出して，
#             どちらかにNaNがある行を取り除きます．
data0,corr=dataprocess0(data,'temp','power')
# 処理された2つのデータの間の相関係数を表示します．
print("Dataprocess0... Corr=",corr)
data0

Dataprocess0... Corr= 0.005663728756090262

# ２つの時系列データから時系列図を作成する関数
def timeseries(df,name1,name2): 
    # dfのインデックス（時間）をXとする
    X=df.index
    # dfのname1列を指定してデータを取り出し，numpy配列で値をY1に与える．
    Y1=df.loc[:,[name1]].values
    # dfのname1列を指定してデータを取り出し，numpy配列で値をY2に与える．
    Y2=df.loc[:,[name2]].values
    # 時系列図の大きさを指定
    plt.figure(figsize=(20, 10))
    # 1つ目(name1)のグラフを1行1列の1つ目に
    ax1=plt.subplot(1,1,1)
    # 2つ目(name2)のグラフのx軸を共有する
    ax2=ax1.twinx()
    # 1つ目(name1)の時系列 
    ax1.plot(X,Y1,color='blue',label=name1)
    # 2つ目(name2)の時系列 
    ax2.plot(X,Y2,color='red',label=name2)
    # グラフのタイトル
    ax1.set_title("Timeseries:"+name1+" and "+name2)
    # x軸のラベル
    ax1.set_xlabel('Time')
    # y軸（左側の第1軸）のラベル
    ax1.set_ylabel('Temperature [degree C]')
    # y軸（右側の第2軸）のラベル
    ax2.set_ylabel('Electric Power [$10^6$kW]')
    # 1つ目(name1)の凡例（左上に置く） 
    ax1.legend(loc='upper left')
    # 2つ目(name1)の凡例（右上に置く）
    ax2.legend(loc='upper right')
    return

# 処理されたデータを用いて時系列図を作成します．
timeseries(data0,'temp','power')

# ２つの時系列データから散布図を作成する関数
def scatter(df,name1,name2):
    # ここがポイント！
    # scikit-learnの線形回帰モデルのクラスを呼び出す
    clf = linear_model.LinearRegression()
    # 説明変数Xにはname1を割り当てる（numpy配列）
    X=df.loc[:,[name1]].values
    # 説明変数Yにはname2を割り当てる（numpy配列）
    Y=df.loc[:,[name2]].values
    # ここがポイント！
    # Y=aX+bの予測モデルを作成する
    clf.fit(X,Y)
    # ここがポイント！
    # 回帰係数a
    slope=clf.coef_
    # ここがポイント！
    # 切片b
    intercept=clf.intercept_
    # ここがポイント！
    # 決定係数R2（回帰直線の当てはまりの良さ）
    r2=clf.score(X,Y)
    # 文字列"Y=aX+b (R2=r2)"
    equation="  y = "+str('{:.1f}'.format(slope[0][0]))+" x +"+str('{:.0f}'.format(intercept[0]))+" (R2="+str('{:.2f}'.format(r2))+")"
    # 散布図の大きさを指定
    plt.figure(figsize=(8, 8))
    # 散布図のプロット
    plt.plot(X, Y, 'o')
    # ここがポイント！
    # 散布図上に回帰直線を引く
    plt.plot(X, clf.predict(X))
    # 文字列"Y=aX+b (R2=r2)"を図の左上に置く
    plt.text(np.nanmin(X), np.nanmax(Y), equation)
    # グラフのタイトル
    plt.title("Scatter diagram:"+name1+" and "+name2)
    # x軸のラベル
    plt.xlabel('Temperature [degree C]')
    # y軸のラベル
    plt.ylabel('Electric Power [$10^6$kW]')    
    return

# 処理されたデータを用いて散布図を作成します．
scatter(data0,'temp','power')

# ２つの時系列データから夏の期間（7/1から9/30まで）を取り出して，
# DataFrameに割り当てる関数
def dataprocess1(df,name1,name2):
    # ここがポイント！
    # 夏の期間（7/1から9/30まで）のデータを取り出します．
    df1=df[(df.index.date >= date(2017,7,1) ) & (df.index.date <= date(2017,9,30))]
    # 処理した2つのname1とname2のデータを取り出し，
    # 2つのデータ間の相関係数を計算します．
    corr=df1[name1].corr(df1[name2])
    # 処理したデータと相関係数を返します．．    
    return df1,corr

# データ前処理１：　dataprocess1関数を呼び出して，
#             夏の期間（7/1から9/30まで）だけを取り出します．
data1,corr=dataprocess1(data0,'temp','power')
# 処理された2つのデータの間の相関係数を表示します．
print("Dataprocess1... Corr=",corr)
data1

Dataprocess1... Corr= 0.7510366253230192

# 処理されたデータを用いて時系列図を作成します．
timeseries(data1,'temp','power')
# 処理されたデータを用いて散布図を作成します．
scatter(data1,'temp','power')

# ２つの時系列データから12～15時だけを取り出して，
# DataFrameに割り当てる関数
def dataprocess2(df,name1,name2):
    # ここがポイント！
    # 14時のデータだけを取り出します．
    df1=df[(df.index.hour >= 12)&((df.index.hour <= 15))]
    # 処理した2つのname1とname2のデータを取り出し，
    # 2つのデータ間の相関係数を計算します．
    corr=df1[name1].corr(df1[name2])
    # 処理したデータと相関係数を返します．
    return df1,corr

# データ前処理2：　dataprocess2関数を呼び出して，
#             12～15時だけを取り出します．
data2,corr=dataprocess2(data1,'temp','power')
# 処理された2つのデータの間の相関係数を表示します．
print("Dataprocess2... Corr=",corr)
data2

Dataprocess2... Corr= 0.7821019235290689

# 処理されたデータを用いて時系列図を作成します．
timeseries(data2,'temp','power')
# 処理されたデータを用いて散布図を作成します．
scatter(data2,'temp','power')

# ２つの時系列データから土曜日・日曜日のデータを取り除いて，
# DataFrameに割り当てる関数
def dataprocess3(df,name1,name2):
    # ここがポイント！
    # 土曜日（weekday==5）と日曜日（weekday==6）を取り除く．
    df1=df[(df.index.weekday != 5) & (df.index.weekday != 6) ]
    # 処理した2つのname1とname2のデータを取り出し，
    # 2つのデータ間の相関係数を計算します．
    corr=df1[name1].corr(df1[name2])
    # 処理したデータと相関係数を返します．
    return df1,corr

# データ前処理3：　dataprocess3関数を呼び出して，
#             土曜日・日曜日を取り除きます．
data3,corr=dataprocess3(data2,'temp','power')
# 処理された2つのデータの間の相関係数を表示します．
print("Dataprocess3... Corr=",corr)
data3

Dataprocess3... Corr= 0.8804545257470542

# 処理されたデータを用いて時系列図を作成します．
timeseries(data3,'temp','power')
# 処理されたデータを用いて散布図を作成します．
scatter(data3,'temp','power')

# ２つの時系列データからお盆や祝日のデータを取り除いて，
# DataFrameに割り当てる関数
def dataprocess4(df,name1,name2):
    # ここがポイント！
    # お盆休暇（8/7から8/20まで）以外残す
    df1=df[(df.index.date < date(2017,8,7) ) | (df.index.date > date(2017,8,20))] 
    # ここがポイント！
    # 海の日以外を残す
    df1=df1[(df1.index.date != date(2017,7,17))] 
    # ここがポイント！
    # 山の日以外を残す
    df1=df1[(df1.index.date != date(2017,8,11))] 
    # ここがポイント！
    # 敬老の日以外を残す
    df1=df1[(df1.index.date != date(2017,9,18))] 
    # ここがポイント！
    # 秋分の日以外を残す
    df1=df1[(df1.index.date != date(2017,9,23))] 
    # 処理した2つのname1とname2のデータを取り出し，
    # 2つのデータ間の相関係数を計算します．
    corr=df1[name1].corr(df1[name2])
    # 処理したデータと相関係数を返します．
    return df1,corr

# データ前処理4：　dataprocess4関数を呼び出して，
#             お盆・祝日を取り除きます．
data4,corr=dataprocess4(data3,'temp','power')
# 処理された2つのデータの間の相関係数を表示します．
print("Dataprocess4... Corr=",corr)
data4

Dataprocess4... Corr= 0.9133959244892919

# 処理されたデータを用いて時系列図を作成します．
timeseries(data4,'temp','power')
# 処理されたデータを用いて散布図を作成します．
scatter(data4,'temp','power')

	temp
datetime
2017-01-01 01:00:00	5.1
2017-01-01 02:00:00	4.1
2017-01-01 03:00:00	4.0
2017-01-01 04:00:00	3.0
2017-01-01 05:00:00	3.6
...	...
2017-12-31 20:00:00	2.8
2017-12-31 21:00:00	3.3
2017-12-31 22:00:00	2.7
2017-12-31 23:00:00	2.7
2018-01-01 00:00:00	2.3

	power
datetime
2017-01-01 00:00:00	2783
2017-01-01 01:00:00	2634
2017-01-01 02:00:00	2520
2017-01-01 03:00:00	2438
2017-01-01 04:00:00	2389
...	...
2017-12-31 19:00:00	3527
2017-12-31 20:00:00	3443
2017-12-31 21:00:00	3335
2017-12-31 22:00:00	3200
2017-12-31 23:00:00	3141

	temp	power
datetime
2017-01-01 00:00:00	NaN	2783.0
2017-01-01 01:00:00	5.1	2634.0
2017-01-01 02:00:00	4.1	2520.0
2017-01-01 03:00:00	4.0	2438.0
2017-01-01 04:00:00	3.0	2389.0
...	...	...
2017-12-31 20:00:00	2.8	3443.0
2017-12-31 21:00:00	3.3	3335.0
2017-12-31 22:00:00	2.7	3200.0
2017-12-31 23:00:00	2.7	3141.0
2018-01-01 00:00:00	2.3	NaN

	temp	power
datetime
2017-01-01 01:00:00	5.1	2634.0
2017-01-01 02:00:00	4.1	2520.0
2017-01-01 03:00:00	4.0	2438.0
2017-01-01 04:00:00	3.0	2389.0
2017-01-01 05:00:00	3.6	2394.0
...	...	...
2017-12-31 19:00:00	4.0	3527.0
2017-12-31 20:00:00	2.8	3443.0
2017-12-31 21:00:00	3.3	3335.0
2017-12-31 22:00:00	2.7	3200.0
2017-12-31 23:00:00	2.7	3141.0

	temp	power
datetime
2017-07-01 00:00:00	23.9	2730.0
2017-07-01 01:00:00	23.0	2570.0
2017-07-01 02:00:00	23.3	2491.0
2017-07-01 03:00:00	22.7	2488.0
2017-07-01 04:00:00	23.1	2478.0
...	...	...
2017-09-30 19:00:00	19.2	3057.0
2017-09-30 20:00:00	18.7	2929.0
2017-09-30 21:00:00	18.7	2765.0
2017-09-30 22:00:00	18.2	2626.0
2017-09-30 23:00:00	18.6	2464.0

1. アメダス気象データ分析チャレンジ！（Python版）　診断的分析¶

1.1 気象データ分析とは？¶

1.2 使用するデータ¶

1.2.1 気象庁AMeDASの気温データ¶

1.2.2 東京電力の電力消費量データ¶

1.3 Pythonプログラムの構成と実行¶

1.3.1 ライブラリのインポート¶

1.3.2 気象庁AMeDASの気温データ（csv形式）を読む¶

1.3.3 東京電力の電力消費量データ（csv形式）を読む¶

1.3.4 2つのデータを結合する¶

1.3.5 データの整理（NaNを取り除く）¶

1.3.6 時系列図を作成する¶

1.3.7 散布図と回帰直線を作成する¶

1.3.8 データの整理と分析（仮説１：夏の期間のデータの抽出）¶

1.3.9 データの整理と分析（仮説２：12～15時のデータの抽出）¶

1.3.10 データの整理と分析（仮説３：休日のデータの削除）¶

1.3.11 データの整理と分析（仮説４：お盆・祝日のデータの削除）¶

1.4 まとめ¶

1.5 著作権について¶

	temp	power
datetime
2017-07-01 12:00:00	24.1	3382.0
2017-07-01 13:00:00	23.0	3340.0
2017-07-01 14:00:00	22.9	3297.0
2017-07-01 15:00:00	22.7	3220.0
2017-07-02 12:00:00	28.7	3363.0
...	...	...
2017-09-29 15:00:00	23.5	3458.0
2017-09-30 12:00:00	22.2	2981.0
2017-09-30 13:00:00	23.6	2995.0
2017-09-30 14:00:00	23.8	2979.0
2017-09-30 15:00:00	23.1	2948.0

	temp	power
datetime
2017-07-03 12:00:00	31.9	4375.0
2017-07-03 13:00:00	31.3	4578.0
2017-07-03 14:00:00	31.0	4664.0
2017-07-03 15:00:00	30.4	4611.0
2017-07-04 12:00:00	28.3	4150.0
...	...	...
2017-09-28 15:00:00	23.7	3579.0
2017-09-29 12:00:00	22.7	3342.0
2017-09-29 13:00:00	23.2	3463.0
2017-09-29 14:00:00	23.5	3485.0
2017-09-29 15:00:00	23.5	3458.0

1. アメダス気象データ分析チャレンジ！（Python版） 診断的分析¶

1.1 気象データ分析とは？¶

1.2 使用するデータ¶

1.2.1 気象庁AMeDASの気温データ¶

1.2.2 東京電力の電力消費量データ¶

1.3 Pythonプログラムの構成と実行¶

1.3.1 ライブラリのインポート¶

1.3.2 気象庁AMeDASの気温データ（csv形式）を読む¶

1.3.3 東京電力の電力消費量データ（csv形式）を読む¶

1.3.4 2つのデータを結合する¶

1.3.5 データの整理（NaNを取り除く）¶

1.3.6 時系列図を作成する¶

1.3.7 散布図と回帰直線を作成する¶

1.3.8 データの整理と分析（仮説１：夏の期間のデータの抽出）¶

1.3.9 データの整理と分析（仮説２：12～15時のデータの抽出）¶

1.3.10 データの整理と分析（仮説３：休日のデータの削除）¶

1.3.11 データの整理と分析（仮説４：お盆・祝日のデータの削除）¶

1.4 まとめ¶

1.5 著作権について¶

1. アメダス気象データ分析チャレンジ！（Python版）　診断的分析¶