How to avoid Triads

from datetime import timezone as tz
from datetime import timedelta as td
from zoneinfo import ZoneInfo

import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import pandas as pd
from mpl_toolkits.axes_grid1 import make_axes_locatable
from scipy.ndimage import gaussian_filter1d


# There's no offset this time of year, but good practise to
# handle timezones properly
london = ZoneInfo("Europe/London")
df = pd.read_csv("triads.csv", parse_dates=["Start Time"])
df["Start Time"] = df["Start Time"].dt.tz_convert(london)
df.head(9)
# Seriously, just data processing here, skip this cell.
dfs = pd.read_csv("sunsets.csv", parse_dates=["date"])
dfs["date"] = dfs["date"].dt.tz_localize(tz.utc).dt.tz_convert(london)
dfs["sunset_float"] = (
    dfs["Sunset"].str.split(":", expand=True)[0].astype(float)
    + dfs["Sunset"].str.split(":", expand=True)[1].astype(float) / 60
)
dfs["dusk_float"] = (
    dfs["Dusk"].str.split(":", expand=True)[0].astype(float)
    + dfs["Dusk"].str.split(":", expand=True)[1].astype(float) / 60
)

# The data we have is a bit noisy, so a super small smooth will make it look nicer.
dfs["sunset_float"] = gaussian_filter1d(dfs["sunset_float"], 2, mode="nearest")
dfs["dusk_float"] = gaussian_filter1d(dfs["dusk_float"], 2, mode="nearest")
dfs = dfs.iloc[:-4]


# Create a start time all in the same year, and turn the time into an hours float for plotting
df["Date2022"] = df["Start Time"].apply(
    lambda x: x.replace(year=2021) if x.month > 10 else x.replace(year=2022)
)
df["Time"] = df["Start Time"].dt.hour + df["Start Time"].dt.minute / 60

fig, ax = plt.subplots()

# Plot sunset, dusk and shading
ax.plot(dfs["date"], dfs["sunset_float"], label="Sunset", lw=1)
ax.plot(dfs["date"], dfs["dusk_float"], label="Dusk", lw=1)
ax.fill_between(
    dfs["date"],
    dfs["sunset_float"],
    dfs["dusk_float"],
    color="#ea580c",
    alpha=0.1,
    label="Twilight",
)

# Plot actual triads
h = ax.scatter(
    df["Date2022"],
    df["Time"],
    c=df["Demand (MW)"],
    s=df["Start Time"].dt.year - 1988,
    label="Demand",
    zorder=2,
)

# Add colorbar
cax = make_axes_locatable(ax).append_axes("right", size="3%", pad=0.1)
fig.colorbar(h, cax=cax, label="Demand (MW)")

# For fun, lets count up the number of occurences instead of adding a histogram trace
for t, v in df["Time"].value_counts().items():
    ax.annotate(
        f"({v})", (df["Date2022"].max() + td(days=7), t), ha="right", va="center", fontsize=10
    )


df_weather = (
    pd.read_csv("london_weather.csv.gz")
    .rename(columns={"dt_iso": "date"})
    .drop(
        columns=[
            "feels_like",
            "pressure",
            "clouds_all",
            "weather_description",
            "snow_1h",
            "snow_3h",
        ]
    )
)
df_weather["date"] = pd.to_datetime(df_weather["date"], format="%Y-%m-%d %H:%M:%S %z UTC")
df_weather = df_weather.set_index("date").sort_index(ascending=False).fillna(0)

# Create a map of triads
df["hour"] = df["Start Time"].dt.floor("60T")
df_triad = df[["hour"]].assign(triad=1).set_index("hour")
df_all = df_weather.join(df_triad)

results = []
for days in [1, 3, 7]:
    df_weather_smoothed = df_weather.rolling(24 * days).median()
    df_all = df_all.join((df_weather - df_weather_smoothed).add_suffix(f"_change_over_{days}_days"))
df_all = df_all.fillna(0)
df_all.columns = [c.replace("_", " ").title().replace("1 Days", "1 Day") for c in df_all.columns]

df_demand = pd.read_csv("demand.csv.gz", parse_dates=["date_start"], index_col="date_start")
# The report does 5min intervals, so we'll aggregate to a SP
df_demand = (
    df_demand.resample("30T").mean().rename(columns={"rolling_demand": "Rolling Demand"}) / 1e3
)

# Combine our data
dfc = df_demand.join(df.assign(Triad=1).set_index("Start Time")[["Triad"]], how="left").fillna(0)
dfc = dfc.join(df_all.drop(columns="Triad"), how="left").ffill(limit=1)
dfc = dfc["2021-11-15":"2022-02-01"]


fig, (ax, ax2) = plt.subplots(
    figsize=(12, 7), nrows=2, sharex=True, gridspec_kw={"height_ratios": [3, 1], "hspace": 0.1}
)
ax.plot(dfc.index, dfc["Rolling Demand"], label="Demand", c="#047857", lw=0.5)
triads = dfc.loc[dfc["Triad"] == 1, "Rolling Demand"]
ax.scatter(triads.index, triads, c="#f59e0b", s=10, label="Triads", zorder=2)
for t in triads.index:
    ax.axvspan(t - td(days=10), t + td(days=10), color="#fcd34d", alpha=0.2, zorder=1)
ax2.plot(dfc.index, dfc["Temp"], label="Temperature", c="#0ea5e9", lw=0.5)

df_final = (
    dfc.assign(triad_risk=lambda x: x["Rolling Demand"] > x["Rolling Demand"].quantile(0.99))
    .loc[lambda x: x.index.hour.isin([16, 17, 18, 19])]
    .assign(Weekday=lambda x: x.index.dayofweek < 5)
    .assign(Weekend=lambda x: x.index.dayofweek >= 5)
    .assign(
        SP=lambda x: "SP" + (1 + 2 * (x.index.hour + x.index.minute / 60)).astype(int).astype(str)
    )
)
# Adding the one hot encoding for day of week and SP
df_final = pd.concat(
    [
        df_final.reset_index(),
        df_final.index.day_name().str.get_dummies().to_frame().reset_index(drop=True),
    ],
    axis=1,
)
df_final = pd.concat([df_final, df_final["SP"].str.get_dummies().reset_index(drop=True)], axis=1)

# Removing some of the weather features that just clutter the plot
df_final = df_final[
    [c for c in df_final.columns if "Pressure" not in c and "Wind" not in c and "Humidity" not in c]
]


corrs2 = (
    df_final.corrwith(df_final["triad_risk"], method="spearman", numeric_only=True)
    .dropna()
    .sort_values()
    .drop(labels=["Rolling Demand", "Triad", "triad_risk"])
)
fig, ax = plt.subplots(figsize=(9, 9))
corrs2.plot.barh(ax=ax, facecolor="#7c3aed")
ax.axvline(0, c="k", lw=1)
ax.set_xlabel("Spearman Correlation Coefficient");