Spaces:

SolarWine
/

api

Sleeping

api / scripts /refresh_energy_data.py

Deploy: 2026 sensor migration + redesign + bucket B endpoints

13fc29d verified 6 months ago

9.24 kB

"""Refresh Data/energy_weather_merged.csv by joining TB Plant production
with IMS station-43 weather.
Pulls hourly `production` (Wh) from the TB Plant asset over the requested
date range, resamples the cached IMS CSV to hourly, joins on UTC timestamp,
and adds the time / solar-geometry features the EnergyPredictor consumes.
Usage:
    python -m scripts.refresh_energy_data                   # Dec 2025 → today
    python -m scripts.refresh_energy_data --from 2026-01-01 --to 2026-05-18
Replaces the legacy Air1-sourced workflow — that device no longer exists
in the 2026 fleet.
"""

from __future__ import annotations

import argparse
import sys
from datetime import datetime, timedelta, timezone
from pathlib import Path

import numpy as np
import pandas as pd

_PROJECT_ROOT = Path(__file__).resolve().parent.parent
sys.path.insert(0, str(_PROJECT_ROOT))

try:
    from dotenv import load_dotenv
    load_dotenv(_PROJECT_ROOT / ".env")
except ImportError:
    pass

from src.data.thingsboard_client import ThingsBoardClient

_LAT_RAD = np.radians(30.85)  # Sde Boker
_IMS_CSV = _PROJECT_ROOT / "Data" / "ims" / "ims_merged_15min.csv"
_OUT_CSV = _PROJECT_ROOT / "Data" / "energy_weather_merged.csv"


def _solar_sin_elevation(doy: int, hour_utc: int) -> float:
    dec = np.radians(23.45 * np.sin(np.radians(360 / 365 * (doy - 81))))
    ha = np.radians(15 * (hour_utc + 2 - 12))   # UTC+2 ≈ local solar
    return float(max(0.0, np.sin(_LAT_RAD) * np.sin(dec)
                      + np.cos(_LAT_RAD) * np.cos(dec) * np.cos(ha)))


_TRACKER_DEVICES = ["Tracker501", "Tracker502", "Tracker503", "Tracker509"]


def fetch_tb_tracker_angles(start: datetime, end: datetime) -> pd.DataFrame:
    """Pull hourly mean tracker angle across the 4 trackers.
    Adds two columns:
      - tracker_angle_mean : mean angle (°) across reporting trackers
      - tracker_angle_std  : std dev across trackers (fleet disagreement)
                             — proxy for "one tracker stowed for intervention"
    """
    client = ThingsBoardClient()
    step = timedelta(days=7)
    per_tracker: dict[str, pd.DataFrame] = {}

    for name in _TRACKER_DEVICES:
        frames = []
        cursor = start
        while cursor < end:
            cursor_end = min(cursor + step, end)
            try:
                df = client.get_timeseries(
                    name, ["angle"],
                    start=cursor, end=cursor_end,
                    limit=10_000, interval_ms=3_600_000, agg="AVG",
                )
            except Exception:
                df = pd.DataFrame()
            if not df.empty:
                frames.append(df)
            cursor = cursor_end
        if frames:
            tdf = pd.concat(frames).sort_index()
            tdf.index = pd.to_datetime(tdf.index, utc=True).floor("h")
            tdf = tdf[~tdf.index.duplicated(keep="last")]
            per_tracker[name] = tdf.rename(columns={"angle": name})

    if not per_tracker:
        # No tracker data — return empty so the join becomes a no-op.
        return pd.DataFrame(columns=["tracker_angle_mean", "tracker_angle_std"])

    wide = pd.concat(per_tracker.values(), axis=1)
    return pd.DataFrame({
        "tracker_angle_mean": wide.mean(axis=1, skipna=True),
        "tracker_angle_std":  wide.std(axis=1, skipna=True).fillna(0.0),
    }, index=wide.index)


def fetch_tb_production(start: datetime, end: datetime) -> pd.DataFrame:
    """Pull hourly Plant production (Wh, summed) from ThingsBoard.
    The Plant asset emits `production` every ~5 min as the energy generated
    in that interval. We sum to hourly buckets via TB's SUM aggregation
    (NOT AVG — AVG of 5-min samples is ~12× smaller than the hourly total).
    Returns a DataFrame indexed by hourly UTC timestamp with column
    `production` in Wh.
    """
    client = ThingsBoardClient()
    frames = []
    cursor = start
    step = timedelta(days=7)
    while cursor < end:
        cursor_end = min(cursor + step, end)
        df = client.get_asset_timeseries(
            "Plant",
            keys=["production"],
            start=cursor,
            end=cursor_end,
            limit=10_000,
            interval_ms=3_600_000,
            agg="SUM",
        )
        if not df.empty:
            frames.append(df)
        cursor = cursor_end

    if not frames:
        raise RuntimeError(f"No Plant production data returned for {start} → {end}")

    combined = pd.concat(frames).sort_index()
    combined.index = pd.to_datetime(combined.index, utc=True)
    # TB returns bucket-center timestamps (e.g. 12:30 for the 12:00 hour);
    # floor to hour so we join cleanly against IMS. Dedupe by index after
    # the floor (not by column values — many zero-production hours would
    # otherwise be collapsed to a single row).
    combined.index = combined.index.floor("h")
    combined = combined[~combined.index.duplicated(keep="last")]
    return combined


def load_ims_hourly() -> pd.DataFrame:
    """Load the cached IMS 15-min CSV and resample to hourly."""
    if not _IMS_CSV.exists():
        raise FileNotFoundError(
            f"{_IMS_CSV} not found. Run scripts.download_ims_data first."
        )
    df = pd.read_csv(_IMS_CSV)
    df["timestamp_utc"] = pd.to_datetime(df["timestamp_utc"], utc=True)
    df = df.set_index("timestamp_utc").sort_index()
    hourly = df.resample("1h").mean()
    return hourly


def build_merged(production: pd.DataFrame, weather: pd.DataFrame,
                 trackers: pd.DataFrame | None = None) -> pd.DataFrame:
    """Inner-join production and weather on hourly timestamp and add features.
    TB SUM aggregation omits zero-production hours; we reindex to a contiguous
    hourly span and fill missing production with 0 so the model learns the
    "GHI=0 → production=0" mapping for nighttime.
    """
    # Reindex production to contiguous hourly grid (nighttime → 0 kWh)
    full_idx = pd.date_range(production.index.min(), production.index.max(),
                             freq="1h", tz="UTC")
    production = production.reindex(full_idx).fillna(0.0)
    production.index.name = "timestamp_utc"

    merged = production.join(weather, how="inner")
    if trackers is not None and not trackers.empty:
        merged = merged.join(trackers, how="left")
        # Backfill short gaps then forward-fill so brief telemetry holes don't
        # introduce NaN; long gaps remain NaN and get dropped at training.
        merged[["tracker_angle_mean", "tracker_angle_std"]] = (
            merged[["tracker_angle_mean", "tracker_angle_std"]].bfill(limit=2).ffill(limit=2)
        )

    # Rename IMS columns to match the legacy EnergyPredictor feature names.
    merged = merged.rename(columns={
        "air_temperature_c": "airTemperature",
        "ghi_w_m2": "GSR",
        "wind_speed_ms": "windSpeed",
        "rh_percent": "airHumidity",
    })

    # Time features
    merged["hour"] = merged.index.hour
    merged["month"] = merged.index.month
    merged["day_of_year"] = merged.index.dayofyear
    merged["hour_sin"] = np.sin(2 * np.pi * merged["hour"] / 24)
    merged["hour_cos"] = np.cos(2 * np.pi * merged["hour"] / 24)
    merged["doy_sin"] = np.sin(2 * np.pi * merged["day_of_year"] / 365)
    merged["doy_cos"] = np.cos(2 * np.pi * merged["day_of_year"] / 365)

    # Solar geometry
    merged["sin_elevation"] = [
        _solar_sin_elevation(int(doy), int(h))
        for doy, h in zip(merged["day_of_year"], merged["hour"])
    ]
    merged["clearness"] = np.where(
        merged["sin_elevation"] > 0.05,
        np.minimum(merged["GSR"] / (merged["sin_elevation"] * 1000), 1.5),
        0.0,
    )

    # Target: kWh per hour
    merged["production_kwh"] = merged["production"] / 1000.0

    return merged


def main() -> None:
    p = argparse.ArgumentParser(description="Refresh Data/energy_weather_merged.csv from TB + IMS.")
    p.add_argument("--from", dest="from_date", default="2025-12-01")
    p.add_argument("--to", dest="to_date", default=None,
                   help="End date (UTC, inclusive). Default: today.")
    args = p.parse_args()

    start = datetime.fromisoformat(args.from_date).replace(tzinfo=timezone.utc)
    end = (datetime.fromisoformat(args.to_date).replace(tzinfo=timezone.utc)
           if args.to_date else datetime.now(tz=timezone.utc))

    print(f"Pulling TB Plant production {start.date()} → {end.date()} ...")
    production = fetch_tb_production(start, end)
    print(f"  {len(production):,} hourly rows from TB Plant asset")

    print("Pulling TB tracker angles ...")
    trackers = fetch_tb_tracker_angles(start, end)
    print(f"  {len(trackers):,} hourly rows from {len(_TRACKER_DEVICES)} trackers")

    print("Loading IMS weather ...")
    weather = load_ims_hourly()
    print(f"  {len(weather):,} hourly rows from IMS")

    merged = build_merged(production, weather, trackers)
    merged.index.name = "timestamp_utc"
    merged.to_csv(_OUT_CSV)
    print(f"\nWrote {len(merged):,} rows × {len(merged.columns)} cols → {_OUT_CSV}")
    print(f"  range: {merged.index.min()} → {merged.index.max()}")
    print(f"  non-null production_kwh: {merged['production_kwh'].notna().sum():,}")


if __name__ == "__main__":
    main()

	"""Refresh Data/energy_weather_merged.csv by joining TB Plant production
	with IMS station-43 weather.

	Pulls hourly `production` (Wh) from the TB Plant asset over the requested
	date range, resamples the cached IMS CSV to hourly, joins on UTC timestamp,
	and adds the time / solar-geometry features the EnergyPredictor consumes.

	Usage:
	python -m scripts.refresh_energy_data # Dec 2025 → today
	python -m scripts.refresh_energy_data --from 2026-01-01 --to 2026-05-18

	Replaces the legacy Air1-sourced workflow — that device no longer exists
	in the 2026 fleet.
	"""

	from __future__ import annotations

	import argparse
	import sys
	from datetime import datetime, timedelta, timezone
	from pathlib import Path

	import numpy as np
	import pandas as pd

	_PROJECT_ROOT = Path(__file__).resolve().parent.parent
	sys.path.insert(0, str(_PROJECT_ROOT))

	try:
	from dotenv import load_dotenv
	load_dotenv(_PROJECT_ROOT / ".env")
	except ImportError:
	pass

	from src.data.thingsboard_client import ThingsBoardClient

	_LAT_RAD = np.radians(30.85) # Sde Boker
	_IMS_CSV = _PROJECT_ROOT / "Data" / "ims" / "ims_merged_15min.csv"
	_OUT_CSV = _PROJECT_ROOT / "Data" / "energy_weather_merged.csv"


	def _solar_sin_elevation(doy: int, hour_utc: int) -> float:
	dec = np.radians(23.45 * np.sin(np.radians(360 / 365 * (doy - 81))))
	ha = np.radians(15 * (hour_utc + 2 - 12)) # UTC+2 ≈ local solar
	return float(max(0.0, np.sin(_LAT_RAD) * np.sin(dec)
	+ np.cos(_LAT_RAD) * np.cos(dec) * np.cos(ha)))


	_TRACKER_DEVICES = ["Tracker501", "Tracker502", "Tracker503", "Tracker509"]


	def fetch_tb_tracker_angles(start: datetime, end: datetime) -> pd.DataFrame:
	"""Pull hourly mean tracker angle across the 4 trackers.

	Adds two columns:
	- tracker_angle_mean : mean angle (°) across reporting trackers
	- tracker_angle_std : std dev across trackers (fleet disagreement)
	— proxy for "one tracker stowed for intervention"
	"""
	client = ThingsBoardClient()
	step = timedelta(days=7)
	per_tracker: dict[str, pd.DataFrame] = {}

	for name in _TRACKER_DEVICES:
	frames = []
	cursor = start
	while cursor < end:
	cursor_end = min(cursor + step, end)
	try:
	df = client.get_timeseries(
	name, ["angle"],
	start=cursor, end=cursor_end,
	limit=10_000, interval_ms=3_600_000, agg="AVG",
	)
	except Exception:
	df = pd.DataFrame()
	if not df.empty:
	frames.append(df)
	cursor = cursor_end
	if frames:
	tdf = pd.concat(frames).sort_index()
	tdf.index = pd.to_datetime(tdf.index, utc=True).floor("h")
	tdf = tdf[~tdf.index.duplicated(keep="last")]
	per_tracker[name] = tdf.rename(columns={"angle": name})

	if not per_tracker:
	# No tracker data — return empty so the join becomes a no-op.
	return pd.DataFrame(columns=["tracker_angle_mean", "tracker_angle_std"])

	wide = pd.concat(per_tracker.values(), axis=1)
	return pd.DataFrame({
	"tracker_angle_mean": wide.mean(axis=1, skipna=True),
	"tracker_angle_std": wide.std(axis=1, skipna=True).fillna(0.0),
	}, index=wide.index)


	def fetch_tb_production(start: datetime, end: datetime) -> pd.DataFrame:
	"""Pull hourly Plant production (Wh, summed) from ThingsBoard.

	The Plant asset emits `production` every ~5 min as the energy generated
	in that interval. We sum to hourly buckets via TB's SUM aggregation
	(NOT AVG — AVG of 5-min samples is ~12× smaller than the hourly total).

	Returns a DataFrame indexed by hourly UTC timestamp with column
	`production` in Wh.
	"""
	client = ThingsBoardClient()
	frames = []
	cursor = start
	step = timedelta(days=7)
	while cursor < end:
	cursor_end = min(cursor + step, end)
	df = client.get_asset_timeseries(
	"Plant",
	keys=["production"],
	start=cursor,
	end=cursor_end,
	limit=10_000,
	interval_ms=3_600_000,
	agg="SUM",
	)
	if not df.empty:
	frames.append(df)
	cursor = cursor_end

	if not frames:
	raise RuntimeError(f"No Plant production data returned for {start} → {end}")

	combined = pd.concat(frames).sort_index()
	combined.index = pd.to_datetime(combined.index, utc=True)
	# TB returns bucket-center timestamps (e.g. 12:30 for the 12:00 hour);
	# floor to hour so we join cleanly against IMS. Dedupe by index after
	# the floor (not by column values — many zero-production hours would
	# otherwise be collapsed to a single row).
	combined.index = combined.index.floor("h")
	combined = combined[~combined.index.duplicated(keep="last")]
	return combined


	def load_ims_hourly() -> pd.DataFrame:
	"""Load the cached IMS 15-min CSV and resample to hourly."""
	if not _IMS_CSV.exists():
	raise FileNotFoundError(
	f"{_IMS_CSV} not found. Run scripts.download_ims_data first."
	)
	df = pd.read_csv(_IMS_CSV)
	df["timestamp_utc"] = pd.to_datetime(df["timestamp_utc"], utc=True)
	df = df.set_index("timestamp_utc").sort_index()
	hourly = df.resample("1h").mean()
	return hourly


	def build_merged(production: pd.DataFrame, weather: pd.DataFrame,
	trackers: pd.DataFrame \| None = None) -> pd.DataFrame:
	"""Inner-join production and weather on hourly timestamp and add features.

	TB SUM aggregation omits zero-production hours; we reindex to a contiguous
	hourly span and fill missing production with 0 so the model learns the
	"GHI=0 → production=0" mapping for nighttime.
	"""
	# Reindex production to contiguous hourly grid (nighttime → 0 kWh)
	full_idx = pd.date_range(production.index.min(), production.index.max(),
	freq="1h", tz="UTC")
	production = production.reindex(full_idx).fillna(0.0)
	production.index.name = "timestamp_utc"

	merged = production.join(weather, how="inner")
	if trackers is not None and not trackers.empty:
	merged = merged.join(trackers, how="left")
	# Backfill short gaps then forward-fill so brief telemetry holes don't
	# introduce NaN; long gaps remain NaN and get dropped at training.
	merged[["tracker_angle_mean", "tracker_angle_std"]] = (
	merged[["tracker_angle_mean", "tracker_angle_std"]].bfill(limit=2).ffill(limit=2)
	)

	# Rename IMS columns to match the legacy EnergyPredictor feature names.
	merged = merged.rename(columns={
	"air_temperature_c": "airTemperature",
	"ghi_w_m2": "GSR",
	"wind_speed_ms": "windSpeed",
	"rh_percent": "airHumidity",
	})

	# Time features
	merged["hour"] = merged.index.hour
	merged["month"] = merged.index.month
	merged["day_of_year"] = merged.index.dayofyear
	merged["hour_sin"] = np.sin(2 * np.pi * merged["hour"] / 24)
	merged["hour_cos"] = np.cos(2 * np.pi * merged["hour"] / 24)
	merged["doy_sin"] = np.sin(2 * np.pi * merged["day_of_year"] / 365)
	merged["doy_cos"] = np.cos(2 * np.pi * merged["day_of_year"] / 365)

	# Solar geometry
	merged["sin_elevation"] = [
	_solar_sin_elevation(int(doy), int(h))
	for doy, h in zip(merged["day_of_year"], merged["hour"])
	]
	merged["clearness"] = np.where(
	merged["sin_elevation"] > 0.05,
	np.minimum(merged["GSR"] / (merged["sin_elevation"] * 1000), 1.5),
	0.0,
	)

	# Target: kWh per hour
	merged["production_kwh"] = merged["production"] / 1000.0

	return merged


	def main() -> None:
	p = argparse.ArgumentParser(description="Refresh Data/energy_weather_merged.csv from TB + IMS.")
	p.add_argument("--from", dest="from_date", default="2025-12-01")
	p.add_argument("--to", dest="to_date", default=None,
	help="End date (UTC, inclusive). Default: today.")
	args = p.parse_args()

	start = datetime.fromisoformat(args.from_date).replace(tzinfo=timezone.utc)
	end = (datetime.fromisoformat(args.to_date).replace(tzinfo=timezone.utc)
	if args.to_date else datetime.now(tz=timezone.utc))

	print(f"Pulling TB Plant production {start.date()} → {end.date()} ...")
	production = fetch_tb_production(start, end)
	print(f" {len(production):,} hourly rows from TB Plant asset")

	print("Pulling TB tracker angles ...")
	trackers = fetch_tb_tracker_angles(start, end)
	print(f" {len(trackers):,} hourly rows from {len(_TRACKER_DEVICES)} trackers")

	print("Loading IMS weather ...")
	weather = load_ims_hourly()
	print(f" {len(weather):,} hourly rows from IMS")

	merged = build_merged(production, weather, trackers)
	merged.index.name = "timestamp_utc"
	merged.to_csv(_OUT_CSV)
	print(f"\nWrote {len(merged):,} rows × {len(merged.columns)} cols → {_OUT_CSV}")
	print(f" range: {merged.index.min()} → {merged.index.max()}")
	print(f" non-null production_kwh: {merged['production_kwh'].notna().sum():,}")


	if __name__ == "__main__":
	main()