# -*- coding: utf-8 -*-
"""エアロゾル光学的厚さ（AVHRR Deep Blue V001, NOAA-18, 1° 月次, 2006–2011）。面×観測の時間（月＝周期）。
v4 の表から: 8 段の帯、欠測は補間して脱色＋枠、色は茶とオレンジの中間。"""
import sys, os, json
import numpy as np
sys.path.insert(0, '/home/claude/forest'); sys.path.insert(0, '/home/claude/quakes'); sys.path.insert(0, '/home/claude/stripes')
import forest1 as F, sylvania as S, q1 as Q
from PIL import Image
from scipy import ndimage
import temp2 as T
OUT = '/home/claude/aerosol/out'
PAPER = np.array([255,255,255], np.float32); LANDC = np.array([238,238,235], np.float32)
INK = np.array([168, 96, 34], np.float32)      # 茶とオレンジの中間
BBOX = (-180, -90, 180, 90); SIZE = (1350, 675); SS = 3; LIM = 0.6   # 帯の飽和 AOD
D = np.load('/home/claude/aerosol/data/aod.npz'); AOD = D['aod'][:, ::-1, :]; CNT = D['cnt'][:, ::-1, :]; MONTHS = list(D['months'])   # 北が上

def field(idx=None, min_cnt=5):
    """月の集合の平均。count が少ない月は除く。"""
    a = AOD if idx is None else AOD[idx]; c = CNT if idx is None else CNT[idx]
    ok = np.isfinite(a) & (c >= min_cnt)
    v = np.where(ok, a, 0).sum(0) / np.maximum(ok.sum(0), 1)
    return np.ma.masked_where(ok.sum(0) == 0, v.astype(np.float32))

def render(name, g, ss=SS, grat=True, fill=True, sigma_deg=1.0, ring=True):
    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
    v, conf = T.fill_field(g, up=10, sigma_deg=sigma_deg*5)     # fill_field は 5° 基準なので ×5
    fr = F.sample_bilinear(np.nan_to_num(v), lon, lat, ok); cf = F.sample_bilinear(conf, lon, lat, ok)
    img = np.zeros((sc.H, sc.W, 3), np.float32); img[:] = PAPER; img[land] = LANDC
    a, b = S.to_lab(PAPER), S.to_lab(INK)
    cols = np.array([S.from_lab(a + (b-a)*(k+0.5)/8) for k in range(8)], np.float32)
    grey = np.array([T.grey_out(c) for c in cols], np.float32)
    k = np.clip((fr / LIM * 8).astype(int), 0, 7); on = (fr > LIM/16) & (cf > 0.12)   # 観測から遠すぎる所（極夜など）は紙のまま
    col = cols[k] * cf[..., None] + grey[k] * (1 - cf[..., None])
    img = np.where(on[..., None], col, img)
    # 海岸線: 場が陸も海も覆うので、地理は極細の相対シフト線で置く（気温 V4 と同じ）
    coast = (land ^ ndimage.binary_erosion(land, iterations=max(1, ss))).astype(np.float32)
    if ring:
        interp = cf < 0.5
        rg = interp & ~ndimage.binary_erosion(interp, iterations=max(1, ss)) & on
        img[rg] = img[rg]*0.68 + np.array([120,120,124], np.float32)*0.32
    im = Image.fromarray(np.clip(img,0,255).astype(np.uint8))
    F.NAVY = INK; F.WHITE = PAPER; im = F.comp_rel(im, coast, 0.22, 0.16)
    if grat:
        from nyt import solid_grat
        F.NAVY = INK; F.WHITE = PAPER; im = F.comp_rel(im, F.aa_mask((sc.W, sc.H), solid_grat(sc)), 0.14, 0.10)
    im = im.resize(SIZE, Image.LANCZOS); im.save(f'{OUT}/{name}.png'); print(name, flush=True); return im

if __name__ == '__main__':
    render('A1', field())                                        # 6 年の平均
    jan = [i for i, m in enumerate(MONTHS) if m.endswith('01')]; jul = [i for i, m in enumerate(MONTHS) if m.endswith('07')]
    render('A2', field(jan)); render('A3', field(jul))           # 1 月と 7 月の平均（周期）