Source code for pyedautils.data_prep.solar_influence

"""Detect direct solar influence on a room temperature sensor.

Given a room temperature series, the local outdoor temperature and the
global radiation at the same location, three signals are computed:

1. **Pearson correlation** between global radiation and the excess
   temperature ``dT = t_room - t_outdoor`` (hourly means). A high positive
   correlation indicates direct sun exposure.
2. **Time-of-day peak**: the local hour with the highest mean ``dT``
   during sunny hours, hinting at the sensor orientation (morning ≈ east,
   noon ≈ south, evening ≈ west).
3. **Steep-rise events**: hours with a strong room-temperature rise,
   strong radiation and a large excess temperature — the typical pattern
   of sun hitting a sensor directly.
"""

from __future__ import annotations

from typing import Dict, Optional, Union

import numpy as np
import pandas as pd


def _ns(x):
    """Return a nanosecond-resolution DatetimeIndex from a Series or Index.

    Different datetime units (e.g. ``[us]`` from Parquet) fail to align in
    :func:`pandas.concat`; nanoseconds is the common ground. No-op on
    pandas < 2.0. Works for tz-naive and tz-aware input.
    """
    idx = pd.DatetimeIndex(pd.to_datetime(x))
    as_unit = getattr(idx, "as_unit", None)
    return as_unit("ns") if as_unit is not None else idx


def _hourly(data: Union[pd.Series, pd.DataFrame]) -> pd.Series:
    """Hourly mean of a ``[timestamp, value]`` DataFrame or a Series."""
    if isinstance(data, pd.Series):
        s = data.copy()
        s.index = _ns(s.index)
    else:
        df = data.copy()
        df.columns = ["timestamp", "value"]
        df = df.set_index(_ns(df["timestamp"]))
        s = df["value"]
    return s.resample("1h").mean()


[docs] def analyze_solar_influence( room: Union[pd.Series, pd.DataFrame], outdoor: Union[pd.Series, pd.DataFrame], radiation: Union[pd.Series, pd.DataFrame], *, local_tz: Optional[str] = None, daytime_radiation: float = 50.0, sunny_radiation: float = 100.0, event_radiation: float = 200.0, event_gradient: float = 1.5, event_dt: float = 3.0, min_hours: int = 48, ) -> Dict[str, Optional[float]]: """Quantify direct solar influence on a room temperature sensor. Args: room: Room temperature, ``[timestamp, value]`` or Series [°C]. outdoor: Outdoor temperature, ``[timestamp, value]`` or Series [°C]. radiation: Global radiation, ``[timestamp, value]`` or Series [W/m²]. local_tz: Optional timezone (e.g. ``"Europe/Zurich"``) used only for the ``peak_hour_local``. Input is assumed UTC; if *None* the peak hour is reported in the input's own clock. daytime_radiation: Radiation above which an hour counts as daytime (denominator of the event rate). Default 50 W/m². sunny_radiation: Radiation above which an hour counts as sunny (time-of-day profile). Default 100 W/m². event_radiation: Radiation threshold for a steep-rise event. Default 200 W/m². event_gradient: Room-temperature rise threshold [K/h] for an event. Default 1.5. event_dt: Excess-temperature threshold [°C] for an event. Default 3.0. min_hours: Minimum overlapping hours required; below this all metrics are *None*. Default 48. Returns: dict: ``n_hours``, ``pearson_r``, ``peak_hour_local``, ``n_solar_events``, ``n_daytime_hours``, ``event_rate``. Values are *None* when fewer than *min_hours* overlapping hours are available. """ room_h = _hourly(room).rename("t_room") out_h = _hourly(outdoor).rename("t_out") sol_h = _hourly(radiation).rename("solar") joined = pd.concat([room_h, out_h, sol_h], axis=1, sort=True).dropna() if len(joined) < min_hours: return { "n_hours": int(len(joined)), "pearson_r": None, "peak_hour_local": None, "n_solar_events": None, "n_daytime_hours": None, "event_rate": None, } joined["dT"] = joined["t_room"] - joined["t_out"] joined["dT_dt"] = joined["t_room"].diff() # 1. Pearson correlation (scipy-free). pearson_r = float(np.corrcoef(joined["solar"], joined["dT"])[0, 1]) # 2. Time-of-day profile during sunny hours. sunny = joined[joined["solar"] > sunny_radiation] if not sunny.empty: idx = sunny.index if local_tz is not None: tz_idx = idx.tz_localize("UTC").tz_convert(local_tz) \ if idx.tz is None else idx.tz_convert(local_tz) hours = tz_idx.hour else: hours = idx.hour hourly_dt = sunny.groupby(hours)["dT"].mean() peak_hour = int(hourly_dt.idxmax()) else: peak_hour = None # 3. Steep-rise events. events = joined[ (joined["dT_dt"] > event_gradient) & (joined["solar"] > event_radiation) & (joined["dT"] > event_dt) ] n_daytime = int((joined["solar"] > daytime_radiation).sum()) n_events = int(len(events)) event_rate = n_events / n_daytime if n_daytime > 0 else 0.0 return { "n_hours": int(len(joined)), "pearson_r": round(pearson_r, 3), "peak_hour_local": peak_hour, "n_solar_events": n_events, "n_daytime_hours": n_daytime, "event_rate": round(event_rate, 4), }