# -*- coding: utf-8 -*-
"""砂塵の全世界地図（GEOS-FP DUEXTTAU）。全球に値がある場なので v4 §3 の枠は ±90 の 2:1。"""
import sys, os
import numpy as np
sys.path.insert(0, '/home/claude/dust'); sys.path.insert(0, '/home/claude/forest'); sys.path.insert(0, '/home/claude/quakes')
import dust as Dm, forest1 as F, q1 as Q
from PIL import Image
from scipy import ndimage
OUT = '/home/claude/dust/out'
PAPER = np.array([255,255,255], np.float32); LG = np.array([238,238,235], np.float32)
def scene(bbox, size, ss):
    sc = F.Scene(bbox, (size[0]*ss, size[1]*ss), proj='pc', pad=0.0); lon, lat, ok = sc.lonlat_grid()
    land = F.sample(Q.land_grid_all().astype(np.float32), lon, lat, ok) > 0.5
    return sc, lon, lat, ok, land
def base_of(sc, land, ss, ink):
    img = np.zeros((sc.H, sc.W, 3), np.float32); img[:] = PAPER; img[land] = LG
    coast = (land ^ ndimage.binary_erosion(land, iterations=ss)).astype(np.float32)
    F.NAVY = ink; F.WHITE = PAPER
    return np.array(F.comp_rel(Image.fromarray(img.astype(np.uint8)), coast, 0.22, 0.16)).astype(np.float32)
def render(name, t, bbox=(-180,-90,180,90), size=(1350,675), ss=3, lim=1.2, gamma=0.8, ink=Dm.INK, key='DUEXTTAU', smooth=1.0):
    sc, lon, lat, ok, land = scene(bbox, size, ss)
    base = base_of(sc, land, ss, ink)
    a = ndimage.gaussian_filter(np.nan_to_num(Dm.grid(key, t)), smooth)
    fr = F.sample_bilinear(a, lon, lat, ok)
    al = (np.clip(fr/lim, 0, 1)**gamma)[..., None]
    im = Image.fromarray(np.clip(base*(1-al) + ink*al, 0, 255).astype(np.uint8)).resize(size, Image.LANCZOS)
    im.save(f'{OUT}/{name}.png'); print(name, flush=True); return im
if __name__ == '__main__':
    t = '20250601_0130'
    ims = [render('W1', t),                                                     # ±90 2:1
           render('W2', t, bbox=(-180,-60,180,75), size=(1350, int(round(1350*135/360)))),   # 砂塵が及ぶ範囲
           render('W3', t, lim=0.8, gamma=1.0),                                 # 濃く
           render('W4', t, key='TOTEXTTAU', lim=1.2)]                           # 全エアロゾル（砂塵＋海塩＋炭素＋硫酸塩）
    y = 0; s = Image.new('RGB', (1350, sum(im.size[1]+10 for im in ims)), (200,200,200))
    for im in ims: s.paste(im, (0, y)); y += im.size[1]+10
    s.save(f'{OUT}/Wsheet.png')