# -*- coding: utf-8 -*-
"""サハラの砂塵（GEOS-FP 3 時間ごと、0.25°、DUEXTTAU＝砂塵の光学的厚さ）。面×観測の時間（事象の尺度＝日）。"""
import sys, os, glob
import numpy as np
sys.path.insert(0, '/home/claude/forest'); sys.path.insert(0, '/home/claude/quakes')
import forest1 as F, sylvania as S, q1 as Q
from PIL import Image
from scipy import ndimage
OUT = '/home/claude/dust/out'
INK = np.array([168, 96, 34], np.float32)
BBOX = (-110, -5, 25, 45); SIZE = (1350, int(round(1350*50/135))); SS = 3
def grid(k, t):
    d = np.load(f'/home/claude/dust/data/{t}.npz'); a = d[k]           # lat -90..90 (南が先), lon -180..180
    return a[::-1]                                                        # 北を上に
def scene(ss=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 render(name, t, paper, landc, ink=INK, lim=1.2, gamma=0.8, dark=False, ss=SS, smooth=1.0):
    sc, lon, lat, ok, land = scene(ss)
    a = ndimage.gaussian_filter(np.nan_to_num(grid('DUEXTTAU', t)), smooth)
    fr = F.sample_bilinear(a, lon, lat, ok)
    img = np.zeros((sc.H, sc.W, 3), np.float32); img[:] = paper; img[land] = landc
    coast = (land ^ ndimage.binary_erosion(land, iterations=ss)).astype(np.float32)
    F.NAVY = ink; F.WHITE = paper
    img = np.array(F.comp_rel(Image.fromarray(img.astype(np.uint8)), coast, 0.22, 0.16)).astype(np.float32)
    al = np.clip(fr/lim, 0, 1)**gamma
    al = al[..., None]
    img = img*(1-al) + ink*al
    im = Image.fromarray(np.clip(img,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'
    WHITE = np.array([255,255,255],np.float32); LG = np.array([238,238,235],np.float32)
    NAVY = np.array(F.r.col(12,270,16),np.float32); NL = np.array(F.r.col(18,268,14),np.float32)
    DARK = np.array([38,36,34],np.float32); DL = np.array([52,50,47],np.float32)
    render('D1', t, WHITE, LG)                                   # 白紙・霞
    render('D2', t, NAVY, NL, ink=np.array([214,150,86],np.float32))         # ネイビー紙・明るい砂色
    render('D3', t, DARK, DL, ink=np.array([214,150,86],np.float32))         # 暗い灰の紙（NYT の夜の地図に近い）
    render('D4', t, WHITE, LG, gamma=1.2, lim=1.5)               # 白紙・淡い所を抑える