import glob,os,sys,math,csv
import numpy as np
import rasterio
from scipy import ndimage
SRC=sys.argv[1] if len(sys.argv)>1 else "."
OUT=sys.argv[2] if len(sys.argv)>2 else "plateaux_bruts.csv"
LATMAX=35.0
SLOPE_MAX=0.05
AREA_MIN_KM2=7.0
DROP_MIN=300.0
DROP_RADIUS_PX=3
rows=[]
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.abs(lats)<=LATMAX
    if not keep.any(): continue
    pxy=111.32*abs(dlat)
    pxx=(111.32*abs(dlon)*np.cos(np.radians(lats)))[:,None]
    pxx2=np.broadcast_to(pxx,(H,W))
    zf=np.where(np.isnan(z),0.0,z)
    land=~np.isnan(z)&(z>-100)
    gy=(np.roll(zf,-1,0)-np.roll(zf,1,0))/(2*pxy*1000.0)
    gx=(np.roll(zf,-1,1)-np.roll(zf,1,1))/(2*pxx2*1000.0)
    slope=np.sqrt(gx*gx+gy*gy)
    flat=(slope<=SLOPE_MAX)&land&(z>0)
    flat[~keep,:]=False
    zmin=ndimage.minimum_filter(np.where(land,zf,1e6),size=2*DROP_RADIUS_PX+1)
    drop=zf-zmin
    rim=flat&(drop>=DROP_MIN)
    near=ndimage.binary_dilation(rim,structure=np.ones((3,3)),iterations=DROP_RADIUS_PX)
    flat=flat&near
    lab,n=ndimage.label(flat,structure=np.ones((3,3)))
    if n==0: continue
    area_px=pxy*pxx2
    areas=ndimage.sum(area_px,lab,index=np.arange(1,n+1))
    maxdrop=ndimage.maximum(np.where(flat,drop,0),lab,index=np.arange(1,n+1))
    rimkm=ndimage.sum(np.where(rim,np.sqrt(area_px),0),lab,index=np.arange(1,n+1))
    meanz=ndimage.mean(zf,lab,index=np.arange(1,n+1))
    cy,cx=zip(*ndimage.center_of_mass(flat,lab,index=np.arange(1,n+1))) if n>0 else ([],[])
    for i in range(n):
        if areas[i]>=AREA_MIN_KM2 and maxdrop[i]>=DROP_MIN:
            la=lat0+dlat*(cy[i]+0.5);lo=lon0+dlon*(cx[i]+0.5)
            rows.append((round(la,3),round(lo,3),round(float(areas[i]),1),round(float(meanz[i])),round(float(maxdrop[i])),round(float(rimkm[i]),1),os.path.basename(f)[:7]))
    print(os.path.basename(f)[:7],"composantes",n,"retenues",sum(1 for i in range(n) if areas[i]>=AREA_MIN_KM2 and maxdrop[i]>=DROP_MIN),flush=True)
with open(OUT,"w",newline="") as fo:
    w=csv.writer(fo);w.writerow(["lat","lon","surface_km2","altitude_m","chute_m","rebord_km","tuile"]);w.writerows(rows)
print("total plateaux",len(rows))