"""Session 13 - Agricultural data analysis with Pandas & NumPy.

Goal: describe the climate/soil of the KFS farm and detect drought periods.
Run:  python s13_explore_farm_data.py
Outputs (outputs/): s13_weekly_report.csv, s13_drought_periods.csv, s13_overview.png
"""
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

pd.set_option("display.width", 160)
pd.set_option("display.max_columns", 20)

from common import OUT, RAINFED_DRY as REFILL_POINT, load_hourly

df = load_hourly()
print("Rows:", len(df), "| from", df.index.min(), "to", df.index.max())
print("\n== Descriptive statistics ==")
print(df.describe().T[["mean", "std", "min", "max"]].round(2))

# ---- 1. Grouping by time: monthly climate profile
monthly = df.resample("MS").agg({"air_temp_c": "mean", "air_humidity_pct": "mean",
                                 "precipitation_mm": "sum", "et0_mm": "sum",
                                 "soil_moisture_0_7cm": "mean"}).round(2)
monthly["water_balance_mm"] = (monthly["precipitation_mm"] - monthly["et0_mm"]).round(1)
print("\n== Monthly profile ==")
print(monthly.head(12))

# ---- 2. Diurnal cycle: when is it hottest? (helps choose irrigation time)
diurnal = df.groupby(df.index.hour)[["air_temp_c", "air_humidity_pct", "et0_mm"]].mean().round(2)
print("\n== Hour of hottest air (mean):", int(diurnal["air_temp_c"].idxmax()), ":00")

# ---- 3. Weekly report for one season (summer 2024)
summer = df.loc["2024-06-01":"2024-08-31"]
weekly = summer.resample("W").agg(temp_mean=("air_temp_c", "mean"), temp_max=("air_temp_c", "max"),
                                  rain=("precipitation_mm", "sum"), et0=("et0_mm", "sum"),
                                  soil_mean=("soil_moisture_0_7cm", "mean"),
                                  soil_min=("soil_moisture_0_7cm", "min")).round(3)
weekly.to_csv(OUT / "s13_weekly_report.csv")

# ---- 4. Drought periods: consecutive hours with soil moisture below the refill point
dry = df["soil_moisture_0_7cm"] < REFILL_POINT
block_id = (dry != dry.shift()).cumsum()
periods = (df[dry].groupby(block_id[dry])
           .agg(start=("soil_moisture_0_7cm", lambda s: s.index.min()),
                end=("soil_moisture_0_7cm", lambda s: s.index.max()),
                hours=("soil_moisture_0_7cm", "size"),
                min_vwc=("soil_moisture_0_7cm", "min")))
periods = periods[periods["hours"] >= 48].sort_values("hours", ascending=False)
periods.to_csv(OUT / "s13_drought_periods.csv", index=False)
print(f"\n== Drought periods (>=48 h below {REFILL_POINT} m3/m3): {len(periods)}")
print(periods.head(5).to_string(index=False))

# ---- 5. Figure
fig, ax = plt.subplots(3, 1, figsize=(11, 8), sharex=True)
daily = df.resample("D").mean(numeric_only=True)
ax[0].plot(daily.index, daily["air_temp_c"], color="#e07a1f"); ax[0].set_ylabel("Air temp °C")
ax[1].bar(daily.index, df["precipitation_mm"].resample("D").sum(), color="#2a6fdb"); ax[1].set_ylabel("Rain mm")
ax[2].plot(daily.index, daily["soil_moisture_0_7cm"], color="#13787b")
ax[2].axhline(REFILL_POINT, ls="--", color="crimson", label="refill point"); ax[2].legend()
ax[2].set_ylabel("Soil VWC m³/m³")
fig.suptitle("KFS farm - climate & soil overview (2023-2024)")
fig.tight_layout(); fig.savefig(OUT / "s13_overview.png", dpi=120)
print("\nSaved outputs to", OUT)
print("Correlation soil moisture vs ET0 (daily):",
      np.round(daily["soil_moisture_0_7cm"].corr(df["et0_mm"].resample("D").sum()), 3))