import csv,math,sys,collections
import numpy as np
import rasterio
from rasterio.warp import transform
WC=sys.argv[1] if len(sys.argv)>1 else "wc"
SMOD=sys.argv[2] if len(sys.argv)>2 else "ghs/GHS_SMOD_E2020_GLOBE_R2023A_54009_1000_V1_0.tif"
BUILT=sys.argv[3] if len(sys.argv)>3 else "ghs/GHS_BUILT_S_E2020_GLOBE_R2023A_54009_1000_V1_0.tif"
rows=list(csv.DictReader(open("plateaux_ceinture_soleil.csv",encoding="utf-8")))
lons=[float(r["lon"]) for r in rows];lats=[float(r["lat"]) for r in rows]
pluie=np.zeros(len(rows))
for m in range(1,13):
    with rasterio.open(f"{WC}/wc2.1_5m_prec_{m:02d}.tif") as src:
        a=src.read(1).astype(np.float32);nod=src.nodata;tr=src.transform
    a[a==nod]=np.nan
    for i,(lo,la) in enumerate(zip(lons,lats)):
        c=int((lo-tr[2])/tr[0]);r=int((la-tr[5])/tr[4])
        w=a[max(r-1,0):r+2,max(c-1,0):c+2]
        v=np.nanmean(w) if np.isfinite(w).any() else np.nan
        pluie[i]+=v if v==v else 0.0
xs,ys=transform("EPSG:4326","ESRI:54009",lons,lats)
def fenetre(src,x,y,demi):
    c=int((x-src.transform[2])/src.transform[0]);r=int((y-src.transform[5])/src.transform[4])
    r0=max(r-demi,0);c0=max(c-demi,0)
    return src.read(1,window=((r0,r+demi+1),(c0,c+demi+1)))
hab=np.zeros(len(rows));vil=np.zeros(len(rows));bati=np.zeros(len(rows))
with rasterio.open(SMOD) as s,rasterio.open(BUILT) as b:
    for i,r in enumerate(rows):
        demi=max(2,int(math.ceil(math.sqrt(float(r["surface_km2"]))/2)))
        w=fenetre(s,xs[i],ys[i],demi).astype(np.int16)
        terre=w>10
        n=max(int(terre.sum()),1)
        hab[i]=100.0*((w>=12)&terre).sum()/n
        vil[i]=100.0*((w>=13)&terre).sum()/n
        wb=fenetre(b,xs[i],ys[i],demi).astype(np.float32)
        wb[wb<0]=np.nan
        bati[i]=100.0*np.nanmean(wb)/1e6 if np.isfinite(wb).any() else 0.0
def occupation(h,v,b):
    if h<=20 and b<1: return "libre"
    if h<=60 and b<3: return "partage"
    return "occupe"
def pluie_classe(p):
    return "sans_pluie" if p<75 else "sec" if p<250 else "tempere" if p<1000 else "arrose" if p<1500 else "pacha"
out=[]
for i,r in enumerate(rows):
    out.append(dict(r,pluie_mm=round(float(pluie[i])),habite_pct=round(float(hab[i])),villages_pct=round(float(vil[i])),bati_pct=round(float(bati[i]),1),occupation=occupation(hab[i],vil[i],bati[i]),pluie_classe=pluie_classe(pluie[i])))
with open("plateaux_ceinture_soleil.csv","w",newline="",encoding="utf-8") as f:
    w=csv.DictWriter(f,fieldnames=list(out[0].keys()));w.writeheader();w.writerows(out)
AB=[r for r in out if r["classe"]!="C"]
def S(l): return sum(int(r["sites"]) for r in l)
def pres(r): return int(r["dist_grande_ville_km"])<=150
print("A+B sites",S(AB),"pres d'une ville",S([r for r in AB if pres(r)]))
for occ in ("libre","partage","occupe"):
    l=[r for r in AB if r["occupation"]==occ]
    print("occupation",occ,"plateaux",len(l),"sites",S(l),"pres d'une ville",S([r for r in l if pres(r)]),"dont A",S([r for r in l if pres(r) and r["classe"]=="A"]),"a moins de 50 km",S([r for r in l if int(r["dist_grande_ville_km"])<=50]))
for pc in ("sans_pluie","sec","tempere","arrose","pacha"):
    l=[r for r in AB if r["pluie_classe"]==pc]
    print("pluie",pc,"plateaux",len(l),"sites",S(l),"dont libres ou partages",S([r for r in l if r["occupation"]!="occupe"]),"dont libres ou partages pres d'une ville",S([r for r in l if r["occupation"]!="occupe" and pres(r)]))
U=[r for r in AB if r["occupation"]!="occupe"]
print("utilisables (A+B, libres ou partages) : sites",S(U),"pres d'une ville",S([r for r in U if pres(r)]),"dont avec pluie >= 75 mm",S([r for r in U if pres(r) and int(r["pluie_mm"])>=75]),"dont libres",S([r for r in U if pres(r) and r["occupation"]=="libre"]),"dont libres avec pluie",S([r for r in U if pres(r) and r["occupation"]=="libre" and int(r["pluie_mm"])>=75]))
for taille,lo,hi in (("7-20",7,20),("20-70",20,70),("70+",70,1e9)):
    l=[r for r in U if lo<=float(r["surface_km2"])<hi]
    print("taille",taille,"plateaux",len(l),"sites",S(l))
tot=collections.defaultdict(lambda:[0,0,0,0,0,0,0])
for o in AB:
    t=tot[o["pays"]];s=int(o["sites"]);t[0]+=1;t[1]+=s
    if o["occupation"]!="occupe":
        t[2]+=s
        if pres(o): t[3]+=s
        if int(o["pluie_mm"])>=75 and pres(o): t[4]+=s
    if o["occupation"]=="libre": t[5]+=s
    if o["pluie_classe"]=="pacha": t[6]+=s
with open("plateaux_ceinture_soleil_pays.csv","w",newline="",encoding="utf-8") as f:
    w=csv.writer(f);w.writerow(["pays","plateaux_A_B","sites_A_B","sites_libres_ou_partages","dont_a_moins_de_150_km_d_une_ville","dont_avec_pluie_75_mm","sites_libres","sites_pluie_1500_mm_pacha"])
    for k,t in sorted(tot.items(),key=lambda x:-x[1][3]): w.writerow([k]+t)
cand=[r for r in out if r["classe"]=="A" and pres(r) and int(r["sites"])>=3 and r["occupation"]!="occupe"]
cand.sort(key=lambda r:-(int(r["sites"])*int(r["heures_tracker"])*(1.0 if r["occupation"]=="libre" else 0.7)))
top=[];cnt=collections.Counter()
for r in cand:
    if cnt[r["pays"]]>=6: continue
    cnt[r["pays"]]+=1;top.append(r)
    if len(top)==100: break
with open("plateaux_ceinture_soleil_top100.csv","w",newline="",encoding="utf-8") as f:
    w=csv.DictWriter(f,fieldnames=list(out[0].keys()));w.writeheader();w.writerows(top)
print("candidats",len(cand),"top 100 : sites",S(top),"libres",sum(1 for r in top if r["occupation"]=="libre"),"sans pluie",sum(1 for r in top if r["pluie_classe"]=="sans_pluie"))
old=[(float(r["lat"]),float(r["lon"])) for r in csv.DictReader(open("top100_avant.csv",encoding="utf-8"))] if len(sys.argv)>4 else []
if old:
    new=set((float(r["lat"]),float(r["lon"])) for r in top)
    print("anciens du top 100 encore presents",sum(1 for o in old if o in new))
    anc={(float(r["lat"]),float(r["lon"])):r for r in out}
    print("anciens ecartes :",collections.Counter(anc[o]["occupation"] for o in old if o not in new))