import csv,sys,collections
import rasterio
GPPD=sys.argv[1] if len(sys.argv)>1 else "gppd.csv"
PV=sys.argv[2] if len(sys.argv)>2 else "PVOUT.tif"
SEUIL_MW=float(sys.argv[3]) if len(sys.argv)>3 else 500.0
SOLEIL=[(1800,871.3),(1900,917.4),(2000,959.8),(2100,998.5),(2200,1033.7),(2300,1065.7),(2400,1094.8),(2500,1121.5),(2600,1145.6),(2700,1167.1)]
def interp(t,h):
    h=min(max(h,t[0][0]),t[-1][0])
    for (a,x),(b,y) in zip(t,t[1:]):
        if a<=h<=b: return x+(y-x)*(h-a)/(b-a)
    return t[-1][1]
HOOVER_GWH_PAR_MWC=4268.0/2000.0;HOOVER_H=2375.0;HOOVER_MEUR_PAR_MWC=1500.0/2000.0;PV_PAR_MW_POSTE=1.0
V14S_MEUR_PAR_GWH=0.43;V14S_GWH=1094.0;V14_GWH=885.0;V14_MEUR=364.0
rows=[x for x in csv.DictReader(open(GPPD,encoding="utf-8")) if x["primary_fuel"]=="Hydro" and float(x["capacity_mw"] or 0)>=SEUIL_MW and x["latitude"] and x["longitude"]]
def production(x):
    v=[float(x[k]) for k in ("generation_gwh_2017","generation_gwh_2018","generation_gwh_2019") if x.get(k)]
    if v: return sum(v)/len(v),"mesuree"
    if x.get("estimated_generation_gwh_2017"): return float(x["estimated_generation_gwh_2017"]),"estimee"
    return None,None
coords=[(float(x["longitude"]),float(x["latitude"])) for x in rows]
with rasterio.open(PV) as src:
    pv=[v[0] for v in src.sample(coords)]
out=[]
for x,p in zip(rows,pv):
    cap=float(x["capacity_mw"]);g,src_g=production(x)
    if g is None: continue
    cf=g/(cap*8.76)
    if cf>0.45 or cf<=0: continue
    heures=float(p)*365*1.25 if p==p and p>0 else 0.0
    facteur=1.0 if cf<=0.30 else 0.5
    pv_mwc=round(cap*PV_PAR_MW_POSTE*facteur)
    if heures>=1800:
        typ="V21 Soleil";gwh=pv_mwc*HOOVER_GWH_PAR_MWC*interp(SOLEIL,heures)/interp(SOLEIL,HOOVER_H);meur=pv_mwc*HOOVER_MEUR_PAR_MWC
        econ=gwh*V14S_MEUR_PAR_GWH-meur;sites=gwh/V14S_GWH
    else:
        n=cap*facteur/300.0;typ="V21 vent (V14 sans lac neuf)";gwh=n*V14_GWH;meur=n*V14_MEUR;econ=0.0;sites=0.0;pv_mwc=0
    out.append(dict(pays=x["country"],nom=x["name"],lat=float(x["latitude"]),lon=float(x["longitude"]),poste_mw=cap,gwh_barrage=round(g),facteur_de_charge=round(cf,2),source=src_g,heures_tracker=round(heures),type=typ,pv_mwc=pv_mwc,gwh_v21=round(gwh),meur_v21=round(meur),economie_meur=round(econ),plateaux_evites=round(sites,1)))
out.sort(key=lambda r:-r["gwh_v21"])
with open(f"v21_retenues_{int(SEUIL_MW)}mw.csv","w",newline="",encoding="utf-8") as f:
    w=csv.DictWriter(f,fieldnames=list(out[0].keys()));w.writeheader();w.writerows(out)
S=lambda l,k:sum(r[k] for r in l)
sol=[r for r in out if r["type"]=="V21 Soleil"];vent=[r for r in out if r["type"]!="V21 Soleil"]
print(f"Barrages hydro >= {SEUIL_MW:.0f} MW dans la base : {len(rows)} ; de pointe (facteur de charge 45 % au plus, production connue) : {len(out)}, dont {len(sol)} au soleil (1 800 h et plus) et {len(vent)} temperes")
print(f"V21 Soleil : {S(sol,'pv_mwc'):.0f} MWc, {S(sol,'gwh_v21')/1000:.0f} TWh, {S(sol,'meur_v21')/1000:.0f} Md€, economie sur les V14 Soleil du plan {S(sol,'economie_meur')/1000:.0f} Md€, plateaux evites {S(sol,'plateaux_evites'):.0f} ; mesurees {sum(1 for r in sol if r['source']=='mesuree')}")
print(f"V21 vent : {S(vent,'gwh_v21')/1000:.0f} TWh pour {S(vent,'meur_v21')/1000:.0f} Md€ ({len(vent)} barrages, {sum(round(r['meur_v21']/V14_MEUR) for r in vent)} V14 sans lac neuf)")
print(f"Compte prudent (deux barrages sur trois ont une retenue d'aval) : soleil {S(sol,'gwh_v21')/1000*2/3:.0f} TWh, {S(sol,'meur_v21')/1000*2/3:.0f} Md€, economie {S(sol,'economie_meur')/1000*2/3:.0f} Md€, plateaux evites {S(sol,'plateaux_evites')*2/3:.0f} ; vent {S(vent,'gwh_v21')/1000*2/3:.0f} TWh")
pays=collections.defaultdict(lambda:[0,0.0,0.0,0,0.0,0.0,0.0])
for r in out:
    t=pays[r["pays"]];t[0]+=1
    if r["type"]=="V21 Soleil": t[1]+=r["gwh_v21"]/1000;t[2]+=r["meur_v21"]/1000;t[3]+=1;t[5]+=r["economie_meur"]/1000;t[6]+=r["plateaux_evites"]
    else: t[4]+=r["gwh_v21"]/1000
for k,t in sorted(pays.items(),key=lambda x:-(x[1][1]+x[1][4]))[:30]: print(k,t[0],"soleil",t[3],round(t[1]),"TWh",round(t[2]),"Md€","econ",round(t[5]),"plateaux",round(t[6]),"; vent",round(t[4]),"TWh")
print("--- 30 premieres V21 Soleil")
for r in sol[:30]: print(r["pays"],r["nom"],r["poste_mw"],r["gwh_barrage"],r["facteur_de_charge"],r["source"],r["heures_tracker"],r["pv_mwc"],r["gwh_v21"],r["meur_v21"],r["economie_meur"])
print("--- 15 premieres temperees")
for r in vent[:15]: print(r["pays"],r["nom"],r["poste_mw"],r["gwh_barrage"],r["facteur_de_charge"],r["source"],r["heures_tracker"],r["gwh_v21"],r["meur_v21"])