In [1]:
import os
import random

import pandas as pd
from pandas.plotting import register_matplotlib_converters
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter
import matplotlib.dates as mdates
%matplotlib inline
register_matplotlib_converters()
np.random.seed = 42

Nahrání metadat

Zde nejdříve nahraji metadata (sazbu zákazníka, území atd.) a vyfiltruji data pro dané území.

In [2]:
date_range = slice("2017-01-13 00:00:00","2017-01-14 23:00:00")
AREA = "kunratice"
df_meta = pd.read_csv("../../scripts/data/metadata.csv", usecols=["id", "id.1", "tdd_class", "area"])
df_meta_kunr = df_meta.loc[((df_meta["area"] == AREA))  & ((df_meta["tdd_class"] == 4))]
def f(x):
    x = x.split("_")
    return f"{x[0].strip()} {x[1].strip()}" 
ids = df_meta_kunr["id.1"].apply(f)

Příprava TDD

Zde si připravím datetime index a data pro TDD 4. Zdroj dat pro TDD 4 pro rok 2017 - zdroj

In [3]:
tdd_index = pd.date_range(start='1/1/2017 01:00:00', end="3/1/2017 01:00:00", freq='H')
TDD_04 = pd.read_csv("TDD_04_JAN_FEB.txt")
TDD_04 = TDD_04.set_index(tdd_index)
TDD_04.head()
Out[3]:
norm_consumption
2017-01-01 01:00:00 0.47816
2017-01-01 02:00:00 0.42856
2017-01-01 03:00:00 0.34875
2017-01-01 04:00:00 0.31430
2017-01-01 05:00:00 0.29346

Nahrání dat

Dále si nahraji data ze všech MONTHS

In [4]:
def load_data(folder_path: str, month_set: set) -> pd.DataFrame:
    dfs = []
    all_data = os.listdir(folder_path)
    for file in all_data:
        for month in month_set:
            if file.rfind(month) != -1:
                df_part = pd.read_csv(f"{folder_path}{file}", usecols=['DeviceID', 'Timestamp', 'ActiveEnergyImport'])
                dfs.append(df_part)
    return dfs
In [5]:
MONTHS = {"2017-02", "2017-01",}
dfs = load_data(f"../../scripts/data/{AREA}/", MONTHS)
df = pd.concat(dfs, ignore_index=True)
df.tail(2)
Out[5]:
DeviceID Timestamp ActiveEnergyImport
7245023 ODBERATEL 84 22.01.2017 5:59:00 164460.0
7245024 ODBERATEL 84 22.01.2017 6:00:00 164461.0

Úprava datetime indexu

Nejdříve je potřeba zbavit sloupec Timestamp bílých mezer a dalších znaků (\n atd). Dále se zbavím všech prázdných hodnot.

In [6]:
df.loc[:, "Timestamp"] = df["Timestamp"].apply(lambda x: str(x).strip())
df["Timestamp"] = pd.to_datetime(df["Timestamp"], format="%d.%m.%Y %H:%M:%S")
df = df.set_index("Timestamp")
df = df.dropna()
df.head(2)
Out[6]:
DeviceID ActiveEnergyImport
Timestamp
2017-02-09 08:01:00 TRAFOSTANICE 1 2474940.0
2017-02-09 08:02:00 TRAFOSTANICE 1 2474950.0

Analýza dat

Nejprve si vyfiltruji data podle id hodnot, které mám k dispozici z metadat.

In [7]:
df = df.loc[df["DeviceID"].isin(list(ids))]
df["DeviceID"] = df["DeviceID"].apply(lambda x: "_".join(x.split()+[AREA]))

Dále si udělám skupiny podle sloupce DeviceID. To znamená, že každá skupina bude mít všechny časové záznamy pro dané DeviceID (tudíž pro daného zákazníka).

In [8]:
groups = df.groupby("DeviceID")

Funkce modify_group

Tato funkce pro danou skupinu spočítá sloupec Consumption (tedy spotřeba) a norm_consumption (tedy normovaná spotřeba pro daný časový úsek). Nejprve se ale přesampluje časová řada na specifikovaný interval (H == hodina, 15T == 15 minut).

In [9]:
def modify_group(name, group):
    """
    Adds Consumption and norm_consumption to group
    """
    group.index = group.index.sort_values()
    # resample by hour
    # resample by 15 min -> "15T
    group = group.resample("H").max() 
    group.loc[:, "Consumption"] = group["ActiveEnergyImport"].diff()
    group.loc[group["Consumption"] < 0, "Consumption"] = 0
    group.loc[group["Consumption"] > 5000, "Consumption"] = 0
    max_val = group["Consumption"].max()
    group["norm_consumption"] = group["Consumption"] / max_val
          
    return group

Průběh spotřebitele do grafu

Funkce plot_consumption vynáší do grafu průběh n spotřebitelů pro daný časový úsek (date_range).

In [10]:
def plot_consumption(groups, date_range, n):
    plt.figure(figsize=(20,10))
    ax = plt.gca()
    ax.tick_params(axis = 'both', which = 'major', labelsize = 16)
    i = 0

    date_range = slice("2017-01-02 00:00:00","2017-01-09 23:00:00")
    TDD_test = TDD_04[date_range] 

    for i, (name, group) in enumerate(groups):

        group = group[date_range]
        group = modify_group(name, group)
        plt.plot(group.index, group["Consumption"] ,label=name, linewidth=2)
        print(group["Consumption"].sum())
        if i == n:
              break

    ax.xaxis.set_major_locator(mdates.HourLocator(interval=6))
    ax.xaxis.set_major_formatter(DateFormatter("%d.%m. %H h"))
    ax.yaxis.set_label_text("Spotřeba [Wh]", fontdict={"fontsize": 15})
    ax.xaxis.set_label_text("Datum [den]", fontdict={"fontsize": 15})
    plt.title(f"Průběh spotřeby dvou odběratelů v zóně {AREA.title()}", fontdict={"fontsize": 20})
    plt.xticks(rotation=90)
    plt.legend(fontsize=16)
    plt.show()
    # plt.savefig("graph_two_customers_Kunratice.png") # does not work!!
In [11]:
plot_consumption(groups, date_range, 1)
89879.0
92620.0

Vykreslení korelace do grafu

Tato funkce vykresluje do grafu všechny spotřebitele, jejichž korelace je nad corr_limit. Pokud je better_than == False, tak se vykreslí do grafu průběhy hodnot horší než daný corr_limit.

In [12]:
def plot_corr(groups: pd.core.groupby.generic.DataFrameGroupBy, date_range: tuple, corr_limit: float, better_than: bool = True) -> plt.plot:
    f, ax = plt.subplots(figsize=(20,10))
    ax.tick_params(axis = 'both', which = 'major', labelsize = 16)
    i = 0
    corrs = []
    TDD_test = TDD_04[date_range]
    ax.plot(TDD_test.index, TDD_test["norm_consumption"], label="TDD 4", linewidth=3)

    for name, group in groups:
        group = group[date_range]
        group = modify_group(name, group)
        corr = TDD_test["norm_consumption"].corr(group["norm_consumption"])
        if better_than:
            if corr > corr_limit:
                ax.plot(group.index, group["norm_consumption"] ,label=f"{name}; corr={corr:.2f}")
        else:
            if corr < corr_limit:
                ax.plot(group.index, group["norm_consumption"] ,label=f"{name}; corr={corr:.2f}")

    ax.xaxis.set_major_locator(mdates.HourLocator(interval=2))
    ax.xaxis.set_major_formatter(DateFormatter("%d.%m. %Hh"))
    ax.yaxis.set_label_text("Poměr [-]", fontdict={"fontsize": 15})
    ax.xaxis.set_label_text("Datum [den]", fontdict={"fontsize": 15})
    plt.title(f"Průběhy odběratelů jejichž korelace s TDD 4 {'≥' if better_than else '≤'} {corr_limit} \n v zóně {AREA.title()}", fontdict={"fontsize": 20})
    
    ax.legend(fontsize=16)
    f.autofmt_xdate()
    return ax
In [13]:
date_range = slice("2017-01-13 00:00:00","2017-01-14 23:00:00")
test_plot = plot_corr(groups, date_range, 0.5)
# plot_corr(groups, date_idx, -0.5, False)
In [14]:
plot_corr(groups, date_range, -0.3, False)
Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x134857358>
In [15]:
def add_cols_to_group(grouped, date_range):
    dffs = []
    for i, (name, group) in enumerate(groups):
        group = group[date_range]
        group = modify_group(name, group)
        dffs.append(group)
    
    dff = pd.concat(dffs)
    dff = dff.drop_duplicates().dropna()
    return dff
In [16]:
date_range = slice("2017-01-13 00:00:00","2017-01-14 23:00:00")
new_df = add_cols_to_group(groups, date_range)
a = new_df.groupby(new_df.index).sum().reset_index().set_index("Timestamp")
a.loc[:, "area"] = AREA

Porovnání průběhů za území s TDD

Zde vezmu normované TDD 4 pro daný časový úsek a porovnám se součtem spotřeby všech zákazníků. Ty poté normuji dle maxima v daném časovém úseku

In [17]:
f, ax = plt.subplots(figsize=(20,10))

ax.tick_params(axis = 'both', which = 'major', labelsize = 16)
i = 0


corrs = []
TDD_test = TDD_04[date_range]
ax.plot(TDD_test.index, TDD_test["norm_consumption"], label="TDD 04", linewidth=3)
ax.plot(a.index, a["Consumption"] / a["Consumption"].max(), label=f"Vsichni odberatele v zone {AREA.title()}")
ax.legend(fontsize=16)
ax.xaxis.set_major_locator(mdates.HourLocator(interval=8))
ax.xaxis.set_major_formatter(DateFormatter("%d.%m. %Hh"))
plt.title(f"Porovnani normovaneho prubehu vsech odberatelu v zone {AREA.title()} s TDD 4", fontdict={"fontsize": 20})
ax.yaxis.set_label_text("Pomer [-]", fontdict={"fontsize": 15})
ax.xaxis.set_label_text("Datum [den]", fontdict={"fontsize": 15})
corr = TDD_test["norm_consumption"].corr(a["norm_consumption"])
ax.annotate(f"corr = {round(corr,2)}", xy=(0.3,0.2), xycoords="figure fraction", size=24)

f.autofmt_xdate()
In [18]:
a.to_pickle(f"{AREA}_tdd4_celek.pkl")
In [19]:
b = pd.read_pickle("barrandov_tdd4_celek.pkl")
c = pd.read_pickle("smichov_tdd4_celek.pkl")
In [20]:
b.head()
Out[20]:
ActiveEnergyImport Consumption norm_consumption area
Timestamp
2017-01-13 01:00:00 11877920.0 1881.0 1.828473 barrandov
2017-01-13 02:00:00 11879892.0 1972.0 1.835366 barrandov
2017-01-13 03:00:00 11881705.0 1813.0 1.771781 barrandov
2017-01-13 04:00:00 11883564.0 1859.0 1.730932 barrandov
2017-01-13 05:00:00 11885636.0 2072.0 2.667290 barrandov
In [21]:
def plot_all_zones(dfs, date_range):
    TDD_test = TDD_04[date_range]
    f, ax = plt.subplots(figsize=(20,10))
    ax.plot(TDD_test.index, TDD_test["norm_consumption"], label="TDD 4", linewidth=3)
    for df in dfs:
        ax.tick_params(axis = 'both', which = 'major', labelsize = 16)
        i = 0
        corrs = []
        df["Consumption"] = df["Consumption"].astype(float)
        df["norm_consumption"] = df["norm_consumption"].astype(float)
        corr = TDD_test["norm_consumption"].corr(df["norm_consumption"])
        ax.plot(df.index, df["Consumption"] / df["Consumption"].max(), label=f"Odběratelé v zóně {df['area'].unique()[0].title()}; corr={round(corr, 2)}")
        ax.legend(fontsize=16)
        ax.xaxis.set_major_locator(mdates.HourLocator(interval=2))
        ax.xaxis.set_major_formatter(DateFormatter("%d.%m. %Hh"))
    
        ax.yaxis.set_label_text("Pomer [-]", fontdict={"fontsize": 15})
        ax.xaxis.set_label_text("Datum [den]", fontdict={"fontsize": 15})
        
    
        f.autofmt_xdate()
    plt.title(f"Porovnání normovaného průběhu všech odběratelů ve všech dostupných zónách s TDD 4", fontdict={"fontsize": 20})
In [22]:
date_range = slice("2017-01-13 00:00:00","2017-01-14 23:00:00")
plot_all_zones([a, b, c], date_range)