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"
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=[r for r in csv.DictReader(open("plateaux_ceinture_soleil.csv",encoding="utf-8")) if r["classe"]!="C"]
land=np.where(np.isnan(mos)|(mos<=0),np.nan,mos)
fig,(a1,a2)=plt.subplots(2,1,figsize=(22,12.6),dpi=100)
fig.patch.set_facecolor("white")
def fond(ax):
    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)
    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([95,95+1000/111.32],[-31,-31],color="#111",lw=3);ax.text(95+500/111.32,-29.3,"1 000 km à l'équateur",ha="center",fontsize=9,color="#111")
def S(l): return sum(int(r["sites"]) for r in l)
fond(a1)
col={"libre":"#2e8b3a","partage":"#f2b705","occupe":"#d62828"}
lab={"libre":"libre : moins de 20 % des mailles à plus de 50 habitants par km², moins de 1 % de bâti","partage":"partagé : 20 à 60 % des mailles peuplées, moins de 3 % de bâti","occupe":"occupé : plus de 60 % des mailles peuplées ou plus de 3 % de bâti"}
for k in ("libre","partage","occupe"):
    l=[r for r in rows if r["occupation"]==k]
    a1.scatter([float(r["lon"]) for r in l],[float(r["lat"]) for r in l],s=2.4,c=col[k],lw=0,alpha=0.85,label=f"{lab[k]} : {len(l)} plateaux, {S(l)} emplacements")
a1.legend(loc="lower left",bbox_to_anchor=(0.015,0.01),fontsize=9,frameon=True,markerscale=5)
a1.set_title("L'occupation des sols (GHS-SMOD et GHS-BUILT-S 2020, 1 km) sur les 16 474 plateaux de classe A et B",fontsize=11,color="#222",pad=8)
fond(a2)
colp={"sans_pluie":"#d62828","sec":"#f28c28","tempere":"#2e8b3a","arrose":"#1f77b4","pacha":"#5e3c99"}
labp={"sans_pluie":"moins de 75 mm par an : eau d'appoint obligatoire","sec":"75 à 250 mm : autonome, impluvium et réserve","tempere":"250 à 1 000 mm","arrose":"1 000 à 1 500 mm","pacha":"1 500 mm et plus : site de Pacha"}
for k in ("tempere","arrose","pacha","sec","sans_pluie"):
    l=[r for r in rows if r["pluie_classe"]==k]
    a2.scatter([float(r["lon"]) for r in l],[float(r["lat"]) for r in l],s=2.4,c=colp[k],lw=0,alpha=0.85,label=f"{labp[k]} : {len(l)} plateaux, {S(l)} emplacements")
a2.legend(loc="lower left",bbox_to_anchor=(0.015,0.01),fontsize=9,frameon=True,markerscale=5)
a2.set_title("La pluie annuelle (WorldClim 2.1, 1970-2000, 9 km) sur les mêmes plateaux",fontsize=11,color="#222",pad=8)
fig.suptitle("Les plateaux de la ceinture du soleil croisés avec les gens et la pluie : 39 068 emplacements de V14 Soleil, 34 856 libres ou partagés, 32 262 en compte prudent. www.mulo.world, 26/10/2026, CC BY 4.0",fontsize=12.5,color="#222",weight="bold",y=0.995)
plt.tight_layout(rect=(0,0,1,0.975))
plt.savefig("carte_croisement_plateaux_occupation_pluie.png",dpi=100,facecolor="white")
print("ok")