{"cells":[{"metadata":{},"cell_type":"markdown","source":"## Missing hours in weather data\nExploring the data a bit you can see that there are several missing timestamps on the given weather data files (air temperatures, etc)\n\nI explore here few ways to fill those gaps using interpolation per site (weather station)\n\nThis kernels should provide an easy way to change your pipeline to use these files, so it outputs the same file structure (gzipped) and you only need to change the data source of your kernel to use the cleaned up data."},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# !pip install --upgrade numpy==1.17.3\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os, gc\nimport random\nimport datetime\n\nfrom tqdm import tqdm_notebook as tqdm\n\n# matplotlib and seaborn for plotting\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\nfrom pandas.plotting import register_matplotlib_converters\nregister_matplotlib_converters()\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import KFold, StratifiedKFold\nfrom sklearn import preprocessing\nfrom sklearn.metrics import mean_squared_error, mean_squared_log_error","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"path = '../input/ashrae-energy-prediction'\n# Input data files are available in the \"../input/\" directory.\nfor dirname, _, filenames in os.walk(path):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check the meter averages per weekday"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"building = pd.read_csv(f'{path}/building_metadata.csv')\n\ntrain = pd.read_csv(f'{path}/train.csv', parse_dates=['timestamp'])\ntrain = train[train.meter == 0]  # electricity only\ndel train['meter']\ntrain = train.merge(building, on='building_id', how='left')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"ematrix = train.groupby(['site_id', 'timestamp']).meter_reading.mean().to_frame('reading').reset_index()\nematrix = ematrix.pivot(index='timestamp', columns='site_id', values='reading')\nematrix.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"# norm\nematrix = ematrix - ematrix.mean()\nematrix = ematrix / ematrix.std()\nematrix.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"ematrix.reset_index(inplace=True)\nematrix['weekday'] = ematrix.timestamp.dt.weekday\nematrix['hour'] = ematrix.timestamp.dt.hour\nematrix.set_index('timestamp', inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"hourly_stats = ematrix.groupby(['weekday', 'hour']).mean()\nax = hourly_stats.plot()\nax.figure.set_size_inches(18, 4)\n_ = ax.set_title('Average readings per hour per weekday')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"building.groupby('site_id').primary_use.agg(lambda x:x.value_counts()[:3].to_dict()).to_dict()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check missing temperature data points"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#load training data 2016\nweather = pd.read_csv(f'{path}/weather_train.csv', parse_dates=['timestamp'])\n# pivot to plot\nwmatrix = weather.pivot(index='timestamp', columns='site_id', values='air_temperature')\n# site with largest amount of missing data points\nsite_id = wmatrix.count().idxmin()\n# plot perid\nstart_date, end_date = datetime.date(2016, 1, 1), datetime.date(2016, 1, 9)\nf,ax = plt.subplots(figsize=(18,6))\n_ = wmatrix.loc[start_date:end_date, site_id].plot(ax=ax, label=f'Jan 2016 site:{site_id}')\n\n# load test data 2017-2018\nweather = pd.read_csv(f'{path}/weather_test.csv', parse_dates=['timestamp'])\n\n# shift 2017 to 2016\nweather.timestamp = weather.timestamp - datetime.timedelta(365)\nwmatrix = weather.pivot(index='timestamp', columns='site_id', values='air_temperature')\n_ = wmatrix.loc[start_date:end_date, site_id].plot(ax=ax, label=f'Jan 2017 site:{site_id}', alpha=0.5)\n\n# shift 2018 to 2016\nweather.timestamp = weather.timestamp - datetime.timedelta(365)\nwmatrix = weather.pivot(index='timestamp', columns='site_id', values='air_temperature')\n_ = wmatrix.loc[start_date:end_date, site_id].plot(ax=ax, label=f'Jan 2018 site:{site_id}', alpha=0.5)\n\n_ = plt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Some different interpolation strategies\nCubic interpolation is generally the best option but the data is a bit too noisy for it, so 50/50% mix with linear interpolation seems better"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"#load training data 2016\nweather = pd.read_csv(f'{path}/weather_train.csv', parse_dates=['timestamp'])\n# pivot to plot\nwmatrix = weather.pivot(index='timestamp', columns='site_id', values='air_temperature')\n# site with largest amount of missing data points\nsite_id = wmatrix.count().idxmin()\n# plot perid\nstart_date, end_date = datetime.date(2016, 1, 1), datetime.date(2016, 1, 9)\nf,ax = plt.subplots(figsize=(18,6))\n\n# load test data 2017-2018\nweather_test = pd.read_csv(f'{path}/weather_test.csv', parse_dates=['timestamp'])\n\n# shift 2017 to 2016\nweather_test.timestamp = weather_test.timestamp - datetime.timedelta(365)\nwtmatrix = weather_test.pivot(index='timestamp', columns='site_id', values='air_temperature')\n_ = wtmatrix.loc[start_date:end_date, site_id].plot(ax=ax, label=f'Jan 2017 site:{site_id}', alpha=0.4)\n\n# shift 2018 to 2016\nweather_test.timestamp = weather_test.timestamp - datetime.timedelta(365)\nwtmatrix = weather_test.pivot(index='timestamp', columns='site_id', values='air_temperature')\n_ = wtmatrix.loc[start_date:end_date, site_id].plot(ax=ax, label=f'Jan 2018 site:{site_id}', alpha=0.4)\n\n\n# def fill_with_avg(wmatrix, w=12):\n#     return wmatrix.fillna(wmatrix.rolling(window=w, win_type='gaussian', center=True, min_periods=1).mean(std=2))\n\ndef fill_with_po3(wmatrix):\n    return wmatrix.fillna(wmatrix.interpolate(method='polynomial', order=3))\n\ndef fill_with_lin(wmatrix):\n    return wmatrix.fillna(wmatrix.interpolate(method='linear'))\n\ndef fill_with_mix(wmatrix):\n    wmatrix = (wmatrix.fillna(wmatrix.interpolate(method='linear', limit_direction='both')) +\n               wmatrix.fillna(wmatrix.interpolate(method='polynomial', order=3, limit_direction='both'))\n              ) * 0.5\n    # workaround: fill last NANs with neighbour\n    assert wmatrix.count().min() >= len(wmatrix)-1 # only the first item is missing\n    return wmatrix.fillna(wmatrix.iloc[1])         # fill with second item\n\n\n_ = fill_with_lin(wmatrix).loc[start_date:end_date, site_id].plot(ax=ax, label=f'linear Jan 2016 site:{site_id}', alpha=0.5)\n_ = fill_with_po3(wmatrix).loc[start_date:end_date, site_id].plot(ax=ax, label=f'cubic Jan 2016 site:{site_id}', alpha=0.5)\n_ = fill_with_mix(wmatrix).loc[start_date:end_date, site_id].plot(ax=ax, label=f'mix Jan 2016 site:{site_id}', alpha=0.5)\n_ = wmatrix.loc[start_date:end_date, site_id].plot(ax=ax, label=f'Jan 2016 site:{site_id}')\n\n_ = plt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### same for dew_temperature"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# pivot to plot\ncol = 'dew_temperature'\nwmatrix = weather.pivot(index='timestamp', columns='site_id', values=col)\n# site with largest amount of missing data points\nsite_id = wmatrix.count().idxmin()\n# plot perid\nstart_date, end_date = datetime.date(2016, 1, 1), datetime.date(2016, 1, 12)\nf,ax = plt.subplots(figsize=(18,6))\n\n_ = fill_with_mix(wmatrix).loc[start_date:end_date, site_id].plot(ax=ax, label=f'mix Jan 2016 site:{site_id}', alpha=0.5)\n_ = wmatrix.loc[start_date:end_date, site_id].plot(ax=ax, label=f'Jan 2016 site:{site_id}')\n\n_ = plt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Save the data for later use"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"def fill_temps(weather):\n    df = None\n    for col in ['air_temperature', 'dew_temperature']:\n        filled = fill_with_mix(weather.pivot(index='timestamp', columns='site_id', values=col))\n        filled = filled.sort_index().unstack().to_frame(col)\n        if df is None:\n            df = filled\n        else:\n            df[col] = filled[col]\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"for src in ['train', 'test']:\n    weather = pd.read_csv(f'{path}/weather_{src}.csv', parse_dates=['timestamp'])\n    wf = fill_temps(weather)\n    wf = wf.reset_index().merge(weather[['site_id', 'timestamp', 'cloud_coverage', 'precip_depth_1_hr', 'wind_direction', 'wind_speed']],\n                           how='left', on=['site_id', 'timestamp']).set_index(['site_id', 'timestamp'])\n    for col in ['cloud_coverage', 'precip_depth_1_hr', 'wind_direction', 'wind_speed']:\n        wf.loc[wf[col] < 0, col] = 0\n        wf.fillna(0, inplace=True)\n    wf.to_csv(f'weather_{src}.csv.gz', compression='gzip', float_format='%g')\n!ls -lah *.gz\nwf.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Double check the results"},{"metadata":{"trusted":true},"cell_type":"code","source":"test = pd.read_csv(f'weather_{src}.csv.gz', parse_dates=['timestamp'])\nax = test[test.site_id == 7].set_index('timestamp').groupby('site_id').air_temperature.plot()\n_ = plt.legend()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"ax = weather[weather.site_id == 7].set_index('timestamp').groupby('site_id').air_temperature.plot()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true,"collapsed":true},"cell_type":"code","source":"# buildings\nfrom sklearn.preprocessing import LabelEncoder\nbuilding = pd.read_csv(f'{path}/building_metadata.csv')\nbuilding.primary_use = LabelEncoder().fit_transform(building.primary_use)\ncols =  ['square_feet', 'year_built', 'floor_count']\nbuilding.fillna(building[cols].mean(), inplace=True)\nbuilding.to_csv(f'building_metadata.csv.gz', index=False, compression='gzip', float_format='%g')\n!ls -lh *.gz\nbuilding.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Visualise all meter readings"},{"metadata":{},"cell_type":"markdown","source":"### First a bad case"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"train = pd.read_csv(f'{path}/train.csv', parse_dates=['timestamp'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# pivot to plot\nematrix = train[(train.meter == 0) & (train.building_id <= 104)]. \\\n                pivot(index='timestamp', columns='building_id', values='meter_reading')\n# site with largest amount of missing data points\nbuilding_id = 0 #ematrix.count().idxmin()\n# plot perid\nstart_date, end_date = datetime.date(2016, 1, 1), datetime.date(2016, 12, 15)\nf,ax = plt.subplots(figsize=(18,6))\n\n_ = ematrix.loc[start_date:end_date, building_id].plot(ax=ax, label=f'Jan 2016 site:{building_id}')\n_ = plt.legend()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### All meters"},{"metadata":{"trusted":true},"cell_type":"code","source":"# Load data\ntrain = pd.read_csv(f'{path}/train.csv')\n# Plot missing values per building/meter\nfig, ax = plt.subplots(1, 4, figsize=(20,30))\nfor meter in range(4):\n    df = train[train.meter == meter]\n    missmap = pd.DataFrame(index=[i for i in range(train.building_id.nunique())])\n    missmap = missmap.merge(df.pivot(index='building_id', columns='timestamp', values='meter_reading'),\n                            how='left', left_index=True, right_index=True)\n    missmap = np.sign(missmap)  # -1, 0 or +1\n    ax[meter].set_title(f'Meter {meter}')\n    sns.heatmap(missmap, cmap='Paired', ax=ax[meter], cbar=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Remove zeros from meter 0 and save data"},{"metadata":{"trusted":true},"cell_type":"code","source":"for src in ['train', 'test', 'sample_submission']:\n    df = pd.read_csv(f'{path}/{src}.csv')\n    if src is 'train':\n        df.drop(index=df[(df.meter_reading <= 0) &\n#                          (df.timestamp <= '2016-06') &\n                         (df.meter == 0)].index,\n                inplace=True)\n    df.to_csv(f'{src}.csv.gz', index=False, compression='gzip', float_format='%g')\n!ls -hl *.gz\nwf.describe()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# PCA weather data"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"import matplotlib.patheffects as PathEffects\n\n# Utility function to visualize the outputs of PCA and t-SNE\ndef scatter(x, colors=np.arange(0,16)):\n    # choose a color palette with seaborn.\n    num_classes = len(np.unique(colors))\n    palette = np.array(sns.color_palette(\"hls\", num_classes))\n\n    # create a scatter plot.\n    f = plt.figure(figsize=(8, 8))\n    ax = plt.subplot(aspect='equal')\n    sc = ax.scatter(x[:,0], x[:,1], s=80, c=palette[colors.astype(np.int)])\n    ax.axis('off')\n    ax.axis('tight')\n\n    # add the labels for each digit corresponding to the label\n    txts = []\n    for i in range(num_classes):\n        # Position of each label at median of data points.\n        xtext, ytext = np.median(x[colors == i], axis=0)[:2]\n        txt = ax.text(xtext, ytext, str(i), fontsize=14, alpha=0.4)\n        txt.set_path_effects([\n            PathEffects.Stroke(linewidth=5, foreground=\"w\"),\n            PathEffects.Normal()])\n        txts.append(txt)\n\n    return f, ax, sc, txts","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":false},"cell_type":"code","source":"# load and norm\nwmatrix = fill_with_mix(weather.pivot(index='timestamp', columns='site_id', values='air_temperature'))\nwmatrix = wmatrix - wmatrix.values.mean()\n\n# # dew temps\n# dewmatrix = fill_with_mix(weather.pivot(index='timestamp', columns='site_id', values='dew_temperature'))\n# dewmatrix = dewmatrix - dewmatrix.values.mean()\n# # concat\n# wmatrix = pd.concat([wmatrix, dewmatrix], axis=0)\n\nX = wmatrix.values.T\nX = X / np.linalg.norm(X)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# from sklearn.manifold import TSNE\n\n# tsne_result = TSNE(random_state=42).fit_transform(X)\n# _ = scatter(tsne_result)\n# # useless :/","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.decomposition import PCA\n\npca = PCA(n_components=2)\npca_result = pca.fit_transform(X)\npd.DataFrame(pca_result, columns=['pca1', 'pca2'])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# rotate and plot (cold-up, hot-down)\n_ = scatter(-pca_result[:,::-1])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Sites sorted by mean temperature in 2016\nweather[['site_id', 'air_temperature', 'cloud_coverage', 'dew_temperature', 'precip_depth_1_hr', 'wind_direction', 'wind_speed']] \\\n.groupby('site_id').agg([np.mean, np.max, np.min]).sort_values(by=('air_temperature','mean'), ascending=False)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Available output files"},{"metadata":{"trusted":true},"cell_type":"code","source":"!ls -hl *.gz","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}