import glob,os,csv,sys
import numpy as np
import rasterio
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
SRC=sys.argv[1] if len(sys.argv)>1 else "tiles"
TITRE=sys.argv[2] if len(sys.argv)>2 else ""
F=12
N,S=72,-58
W=360*120//F;H=(N-S)*120//F
mos=np.full((H,W),np.nan,dtype=np.float32)
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
    if os.path.basename(f).startswith("30N090W"):
        b=z[240:360,360:480];b[b==0]=183
    hh,ww=z.shape
    zs=z[:hh//F*F,:ww//F*F].reshape(hh//F,F,ww//F,F)
    zd=np.nanmean(zs,axis=(1,3))
    top=tr[5];left=tr[2]
    r0=int(round((N-top)*120/F));c0=int(round((left+180)*120/F))
    rr=slice(max(r0,0),min(r0+zd.shape[0],H));cc=slice(max(c0,0),min(c0+zd.shape[1],W))
    if rr.stop<=rr.start or cc.stop<=cc.start: continue
    src=zd[max(0,-r0):max(0,-r0)+(rr.stop-rr.start),max(0,-c0):max(0,-c0)+(cc.stop-cc.start)]
    cur=mos[rr,cc]
    mos[rr,cc]=np.where(np.isnan(cur),src,cur)
rows=list(csv.DictReader(open("cretes_vent_v14.csv",encoding="utf-8")))
land=np.where(np.isnan(mos)|(mos<=0),np.nan,mos)
fig,(a1,a2)=plt.subplots(2,1,figsize=(22,17),dpi=100)
fig.patch.set_facecolor("white")
def fond(ax):
    ax.imshow(np.where(np.isnan(land),0.0,1.0),extent=(-180,180,S,N),cmap="Greys",vmin=-0.2,vmax=4,origin="upper",aspect="auto")
    ax.imshow(land,extent=(-180,180,S,N),cmap="gray",vmin=-500,vmax=3500,origin="upper",aspect="auto",alpha=0.85)
    ax.axhspan(-35,35,color="white",alpha=0.55,lw=0)
    ax.set_xlim(-180,180);ax.set_ylim(S,N)
    ax.set_xticks(range(-180,181,30));ax.set_yticks([-55,-35,0,35,50,70])
    ax.set_xticklabels([f"{abs(x)}°{'O' if x<0 else 'E' if x>0 else ''}" for x in range(-180,181,30)],fontsize=9,color="#444")
    ax.set_yticklabels([f"{abs(y)}°{'S' if y<0 else 'N' if y>0 else ''}" for y in [-55,-35,0,35,50,70]],fontsize=9,color="#444")
    ax.grid(color="#999",lw=0.4,alpha=0.5)
    x0=-160;ax.plot([x0,x0+1000/(111.32*np.cos(np.radians(-45)))],[-47,-47],color="#111",lw=3);ax.text(x0+500/(111.32*np.cos(np.radians(-45))),-45,"1 000 km à 45°",ha="center",fontsize=9,color="#111")
def S_(l,k): return sum(int(r[k]) for r in l)
def taille(l,k): return [1.2+0.9*min(int(r[k]),12) for r in l]
fond(a1)
col={"A":"#5e3c99","B":"#1f77b4","C":"#f2b705","D":"#d9d9d9","sans_vent":"#ff00ff"}
lab={"A":"classe A, 7 m/s et plus à 150 m","B":"classe B, 6,5 à 7 m/s, la V14 nominale","C":"classe C, 6 à 6,5 m/s, la mixte seulement","D":"classe D, moins de 6 m/s, hors plan","sans_vent":"sans donnée de vent"}
for k in ("D","C","B","A","sans_vent"):
    l=[r for r in rows if r["classe"]==k]
    if not l: continue
    a1.scatter([float(r["lon"]) for r in l],[float(r["lat"]) for r in l],s=taille(l,"emplacements"),c=col[k],lw=0,alpha=0.85,label=f"{lab[k]} : {len(l)} crêtes")
a1.legend(loc="lower left",bbox_to_anchor=(0.43,0.02),fontsize=9.5,frameon=True,markerscale=2)
a1.set_title("Le vent moyen à 150 m (Global Wind Atlas 3) sur les crêtes de 12 km et plus qui dominent un versant de 300 m",fontsize=11,color="#222",pad=8)
fond(a2)
colo={"libre":"#2e8b3a","partage":"#f2b705","occupe":"#d62828"}
labo={"libre":"libre : moins de 20 % de mailles à plus de 50 habitants par km²","partage":"partagé : 20 à 60 % de mailles peuplées","occupe":"occupé : au-delà"}
AB=[r for r in rows if int(r["empl_AB"])>0]
for k in ("libre","partage","occupe"):
    l=[r for r in AB if r["occupation"]==k]
    a2.scatter([float(r["lon"]) for r in l],[float(r["lat"]) for r in l],s=taille(l,"empl_AB"),c=colo[k],lw=0,alpha=0.85,label=f"{labo[k]} : {len(l)} crêtes, {S_(l,'empl_AB')} emplacements A et B")
a2.legend(loc="lower left",bbox_to_anchor=(0.36,0.02),fontsize=9.5,frameon=True,markerscale=2)
a2.set_title("L'occupation des sols (GHS-SMOD 2020, 1 km) sur les crêtes qui ont au moins 12 km à 6,5 m/s",fontsize=11,color="#222",pad=8)
fig.suptitle(TITRE,fontsize=12.5,color="#222",weight="bold",y=0.995)
plt.tight_layout(rect=(0,0,1,0.975))
plt.savefig("carte_cretes_vent_v14.png",dpi=100,facecolor="white")
print("ok")