PVGIS EU Science Hub :
simuler 8 760 heures de production
solaire via API REST

import requests
import pandas as pd
from io import StringIO

PVGIS_BASE = "https://re.jrc.ec.europa.eu/api/v5_2"


def appeler_pvgis_8760h(
    lat: float,
    lon: float,
    kwc: float,
    angle: float,
    azimut: float,         # 0=Sud, -90=Est, 90=Ouest
    pertes: float = 14.0,
    montage: str = "building",
) -> pd.DataFrame:
    """
    Appelle PVGIS seriescalc (v5.2) et retourne un DataFrame de 8 760 lignes.
    
    Colonnes retournées :
        datetime  — Index (UTC, résolution 1h, année de référence 2022)
        P         — Puissance instantanée (W)
        Gb        — Irradiance faisceau plan incliné (W/m²)
        Gd        — Irradiance diffuse plan incliné (W/m²)
        Gr        — Irradiance réfléchie sol (W/m²)
        T2m       — Température ambiante 2m (°C)
    """
    params = {
        "lat":           lat,
        "lon":           lon,
        "peakpower":     kwc,
        "loss":          pertes,
        "angle":         angle,
        "aspect":        azimut,
        "mountingplace": montage,
        "pvtechchoice":  "crystSi",
        "raddatabase":   "PVGIS-SARAH3",
        "outputformat":  "json",
        "browser":       0,
    }

    resp = requests.get(
        PVGIS_BASE + "/seriescalc",
        params=params,
        timeout=60,
    )
    resp.raise_for_status()
    data = resp.json()

    # Les données horaires sont dans data['outputs']['hourly']
    hourly = data["outputs"]["hourly"]
    df = pd.DataFrame(hourly)

    # Convertir la colonne 'time' (format "YYYYDDHHMM") en DatetimeIndex UTC
    df["datetime"] = pd.to_datetime(df["time"], format="%Y%d%m:%H%M", utc=True)
    df = df.set_index("datetime").drop(columns=["time"])
    df.index.name = "datetime"

    return df   # 8 760 lignes × 6 colonnes

def calculer_production_annuelle(df: pd.DataFrame) -> dict:
    """
    À partir du DataFrame 8760h retourné par appeler_pvgis_8760h :
    - Calcule la production totale annuelle (kWh)
    - Agrège par mois (kWh/mois)
    - Identifie les mois de pointe et de creux
    """
    # P est en W → énergie par heure = P * 1h / 1000 = kWh
    df["kwh"] = df["P"] / 1000

    production_annuelle = round(df["kwh"].sum(), 0)

    # Agrégation mensuelle
    mensuel = (
        df["kwh"]
        .resample("MS")   # MS = Month Start (début de mois)
        .sum()
        .round(1)
    )

    MOIS_FR = ["Jan","Fév","Mar","Avr","Mai","Jun",
               "Jul","Aoû","Sep","Oct","Nov","Déc"]

    return {
        "production_kwh_an":   production_annuelle,
        "production_par_mois": dict(zip(MOIS_FR, mensuel.values)),
        "mois_pointe":        mensuel.idxmax().strftime("%B"),
        "mois_creux":         mensuel.idxmin().strftime("%B"),
        "ratio_hiver_ete":    round(mensuel.iloc[:3].mean() / mensuel.iloc[6:9].mean(), 2),
    }


def p50_p90(production_kwh_an: float, incertitude_pct: float = 6.0) -> dict:
    """
    Calcule P50 et P90 à partir de la production annuelle PVGIS.
    
    Convention bancaire : P90 = P50 × (1 - 1.28 × sigma)
    Incertitude par défaut : ±6% (recommandation ACEA/Société Générale CIB)
    """
    import math
    sigma = incertitude_pct / 100
    p50   = production_kwh_an
    p90   = p50 * (1 - 1.28 * sigma)  # Loi normale, 90e percentile

    return {
        "P50_kwh_an": round(p50, 0),
        "P90_kwh_an": round(p90, 0),
        "P90_P50_ratio_pct": round((p90 / p50 - 1) * 100, 1),
        "incertitude_1sigma_pct": incertitude_pct,
    }

from itertools import product
import time


def optimiser_inclinaison(
    lat: float, lon: float, kwc: float,
    angles: list[int] = [0,10,15,20,25,30,35],
    azimuts: list[int] = [-45,0,45],
) -> dict:
    """
    Teste n angles × m azimuts et retourne la combinaison optimale.
    
    Nota : ~21 appels API (throttling 1 requête/s pour respecter les CGU PVGIS).
    """
    resultats = []

    for angle, azimut in product(angles, azimuts):
        df = appeler_pvgis_8760h(lat, lon, kwc, angle, azimut)
        prod = df["kwh"].sum() if "kwh" in df. else df["P"].sum()/1000
        resultats.append({
            "angle": angle, "azimut": azimut, "production_kwh_an": prod
        })
        time.sleep(1.1)  # ← Respecter la limite PVGIS (1 req/s)

    return max(resultats, key=lambda x: x["production_kwh_an"])

def bilan_pv_complet(
    lat: float, lon: float, kwc: float,
    angle: float, azimut: float,
    access_token_enedis: str,
    usage_point_id: str,
    prix_kwh_achat: float = 0.2516,   # TRV HP 2024
    tarif_injection: float = 0.1326,   # OA CRE 2024 < 36 kWc
) -> dict:

    # 1 – Production horaire PVGIS
    df_pvgis = appeler_pvgis_8760h(lat, lon, kwc, angle, azimut)
    df_pvgis["kwh"] = df_pvgis["P"] / 1000

    # 2 – Courbe de charge Enedis (rééchantillonnée à 1h)
    from datetime import date
    df_enedis = recuperer_courbe_charge(
        access_token_enedis, usage_point_id,
        date_debut=date(2023,1,1), date_fin=date(2023,12,31)
    )
    df_enedis_1h = df_enedis["wh"].resample("1H").sum() / 1000  # → kWh/h
    df_enedis_1h = df_enedis_1h.rename("kwh_conso")

    # 3 – Croisement
    bilan = pd.DataFrame({
        "prod":  df_pvgis["kwh"],
        "conso": df_enedis_1h,
    }).dropna()

    bilan["auto"]     = bilan[["prod","conso"]].min(axis=1)
    bilan["surplus"]  = (bilan["prod"] - bilan["conso"]).clip(0)
    bilan["soutire"]  = (bilan["conso"] - bilan["prod"]).clip(0)

    # 4 – Calcul financier
    eco_auto    = bilan["auto"].sum() * prix_kwh_achat
    revenu_inj  = bilan["surplus"].sum() * tarif_injection
    facture_abo = bilan["soutire"].sum() * prix_kwh_achat

    return {
        "production_kwh_an":          round(bilan["prod"].sum(),0),
        "consommation_kwh_an":        round(bilan["conso"].sum(),0),
        "taux_autoconsommation_pct":  round(bilan["auto"].sum()/bilan["prod"].sum()*100,1),
        "taux_autosuffisance_pct":    round(bilan["auto"].sum()/bilan["conso"].sum()*100,1),
        "economies_autoconso_eur_an": round(eco_auto,0),
        "revenus_injection_eur_an":   round(revenu_inj,0),
        "gain_total_eur_an":          round(eco_auto + revenu_inj,0),
    }