import glob,os,csv
import numpy as np
import rasterio
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
SRC="tiles"
F=12
W=360*120//F;H=70*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
    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((35-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))
    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("plateaux_ceinture_soleil.csv",encoding="utf-8")))
fig,ax=plt.subplots(figsize=(22,6.4),dpi=100)
fig.patch.set_facecolor("white")
land=np.where(np.isnan(mos)|(mos<=0),np.nan,mos)
ax.imshow(np.where(np.isnan(land),0.0,1.0),extent=(-180,180,-35,35),cmap="Greys",vmin=-0.2,vmax=4,origin="upper",aspect="auto")
ax.imshow(land,extent=(-180,180,-35,35),cmap="gray",vmin=-500,vmax=3500,origin="upper",aspect="auto",alpha=0.85)
col={"A":"#e63900","B":"#f2b705","C":"#3f7fc4"}
lab={"A":"classe A : 2 100 h de tracker et plus","B":"classe B : 1 800 à 2 100 h","C":"classe C : moins de 1 800 h (hors V14 Soleil)"}
for cl in "CBA":
    xs=[float(r["lon"]) for r in rows if r["classe"]==cl];ys=[float(r["lat"]) for r in rows if r["classe"]==cl]
    ax.scatter(xs,ys,s=2.2 if cl!="C" else 1.2,c=col[cl],lw=0,alpha=0.85 if cl!="C" else 0.5,label=f"{lab[cl]} : {len(xs)} plateaux")
ax.set_xlim(-180,180);ax.set_ylim(-35,35)
ax.set_xticks(range(-180,181,30));ax.set_yticks(range(-30,31,10))
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 range(-30,31,10)],fontsize=9,color="#444")
ax.grid(color="#999",lw=0.4,alpha=0.5)
ax.plot([-120,-120+1000/111.32],[-32.5,-32.5],color="#111",lw=3);ax.text(-120+500/111.32,-30.8,"1 000 km à l'équateur",ha="center",fontsize=9,color="#111")
ax.legend(loc="lower left",bbox_to_anchor=(0.005,0.1),fontsize=9,frameon=True,markerscale=5)
ax.set_title("Les 21 925 plateaux de la ceinture du soleil (35° N à 35° S) de 7 km² et plus, à pente de 5 % au plus, dominant un versant de 300 m à moins de 3 km : relief GMTED2010 à 1 km, productible du Global Solar Atlas. www.mulo.world, 23/10/2026, CC BY 4.0",fontsize=11,color="#222",pad=10)
plt.tight_layout()
plt.savefig("carte_plateaux_ceinture_soleil.png",dpi=100,facecolor="white")
print("ok")