import glob,os,sys,math,csv,gzip,json
import numpy as np
import rasterio
from scipy import ndimage
from skimage.morphology import skeletonize
SRC=sys.argv[1] if len(sys.argv)>1 else "tiles"
VENT=sys.argv[2] if len(sys.argv)>2 else "gwa150_001deg.bin.gz"
OUT=sys.argv[3] if len(sys.argv)>3 else "cretes_brutes.csv"
BANDES=[(35.0,70.0),(-55.0,-35.0)]
DROP_MIN=300.0
DROP_RADIUS_PX=3
DOMINE=20.0
DOMINE2=150.0
ALT_MAX=2500.0
DROP_MAX=900.0
LONG_MIN=12.0
BLOC=1.0
def charge_vent(p):
    b=gzip.open(p).read()
    n=int.from_bytes(b[:4],"little");h=json.loads(b[4:4+n]);o=4+n;t={}
    for k in h["keys"]:
        t[k]=np.frombuffer(b,dtype=np.uint8,count=1000000,offset=o).reshape(1000,1000);o+=1000000
    return t
V=charge_vent(VENT)
def vent(la,lo):
    R=np.floor((90.0-la)*100).astype(np.int64);C=np.floor((lo+180.0)*100).astype(np.int64)
    k=(R//1000)*100+C//1000;out=np.full(len(R),np.nan)
    for kk in np.unique(k):
        a=V.get(int(kk))
        if a is None: continue
        m=k==kk;v=a[R[m]%1000,C[m]%1000].astype(float);v[v==0]=np.nan;out[m]=v/10.0
    return out
def exclu(la,lo):
    return (la>59.5)&(lo>-60)&(lo<-11)&~((la>62.5)&(la<67.5)&(lo>-25.5)&(lo<-12.5))
rows=[];bilan=[]
for f in sorted(glob.glob(os.path.join(SRC,"*gmted_mea300.tif"))):
    with rasterio.open(f) as r:
        z=r.read(1).astype(np.float32);tr=r.transform;nod=r.nodata
    z[z==nod]=np.nan
    H,W=z.shape
    lat0=tr[5];dlat=tr[4];lon0=tr[2];dlon=tr[0]
    lats=lat0+dlat*(np.arange(H)+0.5)
    keep=np.zeros(H,bool)
    for a,b in BANDES: keep|=(lats>=a)&(lats<=b)
    if not keep.any(): continue
    pxy=111.32*abs(dlat)
    pxx=111.32*abs(dlon)*np.cos(np.radians(lats))
    land=~np.isnan(z)&(z>-100)
    zf=np.where(land,z,np.nan)
    zmin=ndimage.minimum_filter(np.where(land,z,1e6),size=2*DROP_RADIUS_PX+1)
    drop=np.where(land,z-zmin,0)
    zp=np.pad(np.where(land,z,-1e6),2,mode="edge")
    c=zp[2:-2,2:-2]
    def sh(dy,dx): return zp[2+dy:zp.shape[0]-2+dy,2+dx:zp.shape[1]-2+dx]
    faite=np.zeros_like(land)
    for dy,dx in ((1,0),(0,1),(1,1),(1,-1)):
        faite|=(c-sh(dy,dx)>=DOMINE)&(c-sh(-dy,-dx)>=DOMINE)&(c-sh(2*dy,2*dx)>=DOMINE2)&(c-sh(-2*dy,-2*dx)>=DOMINE2)
    crete=land&(drop>=DROP_MIN)&(drop<=DROP_MAX)&faite&(z<=ALT_MAX)
    crete[~keep,:]=False
    crete=ndimage.binary_closing(crete,structure=np.ones((3,3)))&land
    crete[~keep,:]=False
    sq=skeletonize(crete)
    lab,n=ndimage.label(crete,structure=np.ones((3,3)))
    if n==0: continue
    yy,xx=np.nonzero(sq)
    la=lat0+dlat*(yy+0.5);lo=lon0+dlon*(xx+0.5)
    ok=~exclu(la,lo)
    yy,xx,la,lo=yy[ok],xx[ok],la[ok],lo[ok]
    l=lab[yy,xx]
    side=np.sqrt(pxy*pxx[yy])*1.08
    v=vent(la,lo)
    zz=z[yy,xx];dd=drop[yy,xx]
    by=np.floor(la/BLOC).astype(np.int64);bx=np.floor(lo/BLOC).astype(np.int64)
    ob=ndimage.find_objects(lab)
    diam={}
    for i,s in enumerate(ob):
        if s is None: continue
        hy=(s[0].stop-s[0].start)*pxy;cy=lat0+dlat*((s[0].start+s[0].stop)/2)
        hx=(s[1].stop-s[1].start)*111.32*abs(dlon)*math.cos(math.radians(cy))
        diam[i+1]=math.hypot(hy,hx)
    big=np.array([diam.get(int(k),0)>=LONG_MIN for k in range(n+1)])
    m=big[l]
    yy,xx,la,lo,l,side,v,zz,dd,by,bx=[a[m] for a in (yy,xx,la,lo,l,side,v,zz,dd,by,bx)]
    ld=np.array([diam.get(int(k),0) for k in range(n+1)])
    key=l.astype(np.int64)*1000000+(by+100)*1000+(bx+200)
    order=np.argsort(key,kind="stable")
    key=key[order]
    yy,xx,la,lo,l,side,v,zz,dd=[a[order] for a in (yy,xx,la,lo,l,side,v,zz,dd)]
    cut=np.flatnonzero(np.diff(key))+1
    st=np.r_[0,cut];en=np.r_[cut,len(key)]
    t=os.path.basename(f)[:7];nt=0
    for s0,e0 in zip(st,en):
        L=float(side[s0:e0].sum())
        if L<LONG_MIN: continue
        vv=v[s0:e0];ss=side[s0:e0];fin=np.isfinite(vv)
        L65=float(ss[fin&(vv>=6.5)].sum());L7=float(ss[fin&(vv>=7.0)].sum());L6=float(ss[fin&(vv>=6.0)].sum())
        vm=float(np.nanmean(vv)) if fin.any() else float("nan")
        w=np.where(fin,np.nan_to_num(vv),0)+1e-6
        rows.append(dict(lat=round(float(np.average(la[s0:e0],weights=w)),3),lon=round(float(np.average(lo[s0:e0],weights=w)),3),longueur_km=round(L,1),etendue_km=round(float(ld[l[s0]]),1),altitude_m=round(float(zz[s0:e0].mean())),chute_max_m=round(float(dd[s0:e0].max())),vent_150m_moyen=round(vm,2) if vm==vm else "",vent_150m_max=round(float(np.nanmax(vv)),2) if fin.any() else "",part_sans_vent_pct=round(100.0*(~fin).sum()/len(vv)),km_6ms=round(L6,1),km_65ms=round(L65,1),km_7ms=round(L7,1),empl_A=int(L7//12),empl_AB=int(L65//12),empl_ABC=int(L6//12),emplacements=max(1,int(L//12)),tuile=t))
        nt+=1
    print(t,"composantes",n,"crêtes retenues",nt,"mailles de faîte",len(yy),"dont sans vent",int(np.isnan(v).sum()),flush=True)
    bilan.append((t,n,nt,len(yy),int(np.isnan(v).sum())))
with open(OUT,"w",newline="",encoding="utf-8") as fo:
    w=csv.DictWriter(fo,fieldnames=list(rows[0].keys()));w.writeheader();w.writerows(rows)
with open("bilan_tuiles_cretes.csv","w",newline="") as fo:
    w=csv.writer(fo);w.writerow(["tuile","composantes","cretes","mailles_faite","mailles_sans_vent"]);w.writerows(bilan)
print("total crêtes",len(rows),"emplacements A",sum(r["empl_A"] for r in rows),"A+B",sum(r["empl_AB"] for r in rows),"A+B+C",sum(r["empl_ABC"] for r in rows))