# -*- coding: utf-8 -*-
"""地震（USGS ComCat, 1973–2025, M≥5.0）。v4 の表から導いた最初の一枚。
点×観測の時間 → 場＝画素あたりの回数（log1p）、8 段の帯（Q1）。
比較のため、釜田さんの提案（大きさ＝マグニチュード、薄さ＝深さ）を印で描いた Q2。"""
import sys, os, json
import numpy as np
sys.path.insert(0, '/home/claude/forest')
import forest1 as F, sylvania as S
from nyt import solid_grat
from PIL import Image, ImageDraw
from scipy import ndimage
OUT = '/home/claude/quakes/out'; os.makedirs(OUT, exist_ok=True)
PAPER = np.array([255, 255, 255], np.float32); LANDC = np.array([238, 238, 235], np.float32)
INK = np.array([28, 26, 24], np.float32)          # 白から最も遠い色（題材固有の色相は未定）
BBOX = (-180, -90, 180, 90); SIZE = (1350, 675); SS = 3
d = np.load('/home/claude/quakes/quakes.npz'); MMIN = 5.0
sel = d['mag'] >= MMIN
LON, LAT, DEP, MAG, T = d['lon'][sel], d['lat'][sel], d['dep'][sel], d['mag'][sel], d['t'][sel]
print('events', len(LON))

def land_grid_all(ppd=20):
    p = f'/home/claude/quakes/land_all_{ppd}.npy'
    if os.path.exists(p): return np.load(p)
    big = Image.new('L', (360*ppd, 180*ppd), 0); dr = ImageDraw.Draw(big)
    for ring in F.land_polys(): dr.polygon([((lo+180)*ppd, (90-la)*ppd) for lo, la in ring], fill=255)
    g = np.array(big) > 127; np.save(p, g); return g

def scene():
    sc = F.Scene(BBOX, (SIZE[0]*SS, SIZE[1]*SS), proj='pc', pad=0.0)
    lon, lat, ok = sc.lonlat_grid()
    land = F.sample(land_grid_all().astype(np.float32), lon, lat, ok) > 0.5
    return sc, lon, lat, ok, land

def field(cell=0.5, ref_q=0.995):
    """画素あたりの回数。格子 cell°、平滑化は格子 1 個分（上限）。log1p で尺度を固定。"""
    w, h = int(360/cell), int(180/cell)
    c = np.zeros((h, w), np.float32)
    ix = np.clip(((LON+180)/cell).astype(int), 0, w-1); iy = np.clip(((90-LAT)/cell).astype(int), 0, h-1)
    np.add.at(c, (iy, ix), 1)
    cs = ndimage.gaussian_filter(c, 1.0, mode=('nearest', 'wrap'))
    ref = float(np.quantile(cs[cs > 0], ref_q))
    v = np.log1p(cs) / np.log1p(ref)
    meta = dict(cell_deg=cell, sigma_cells=1.0, ref=ref, max_count=float(c.max()), events=int(len(LON)), mmin=MMIN)
    return np.clip(v, 0, 1).astype(np.float32), meta

def bands_img(bands=8):
    sc, lon, lat, ok, land = scene()
    v, meta = field()
    fr = F.sample_bilinear(v, lon, lat, ok)
    img = np.zeros((sc.H, sc.W, 3), np.uint8); img[:] = PAPER.astype(np.uint8); img[land] = LANDC.astype(np.uint8)
    a, b = S.to_lab(PAPER), S.to_lab(INK)
    cols = [S.from_lab(a + (b - a) * (k + 0.5) / bands).astype(np.uint8) for k in range(bands)]
    on = fr > 1.0/bands/2
    for k in range(bands):
        m = on & (fr > k/bands) & (fr <= (k+1)/bands + (1 if k == bands-1 else 0)); img[m] = cols[k]
    im = Image.fromarray(img)
    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}/Q1.png'); json.dump(meta, open(f'{OUT}/Q1.json', 'w'), indent=1); print('Q1', meta, flush=True)

def marks_img(name='Q2', alpha=0.30, r0=0.9, base=1.8, dep_full=300.0):
    """印。半径＝マグニチュード（r0·base^(M−5)）、薄さ＝深さ（浅い＝インク、深い＝紙へ）。"""
    sc, lon, lat, ok, land = scene()
    W, H = sc.W, sc.H
    img = np.zeros((H, W, 3), np.float32); img[:] = PAPER; img[land] = LANDC
    # 線形空間で減法的に重ねる（濃い所は自然に暗くなる）
    acc = np.zeros((H, W), np.float32)   # インクの被覆率
    dep_w = np.zeros((H, W), np.float32) # 深さの重み付き平均のための累積
    order = np.argsort(MAG)              # 小さいものから描く（大きいものを上に）
    x = (LON*111320.0 - sc.cx)/sc.s + W/2; y = H/2 - (LAT*111320.0 - sc.cy)/sc.s
    r = r0 * base ** (MAG - MMIN) * SS
    canvas = Image.new('F', (W, H), 0.0)
    # 深さの階級ごとに別レイヤー（薄さは深さで決まる）
    layers = []
    edges = [0, 35, 70, 150, 300, 800]
    for k in range(len(edges)-1):
        m = (DEP >= edges[k]) & (DEP < edges[k+1])
        cv = Image.new('L', (W, H), 0); dr = ImageDraw.Draw(cv)
        for i in np.nonzero(m)[0]:
            dr.ellipse([x[i]-r[i], y[i]-r[i], x[i]+r[i], y[i]+r[i]], fill=255)
        layers.append((np.array(cv).astype(np.float32)/255, k))
    a_ink, b_pap = S.to_lab(INK), S.to_lab(PAPER)
    depth_L = [0.0, 0.15, 0.35, 0.55, 0.72]   # 深いほど紙へ寄せる割合
    out = img.copy()
    for cov, k in layers[::-1]:   # 深いものから先に、浅いものを上に
        col = S.from_lab(a_ink + (b_pap - a_ink) * depth_L[k]).astype(np.float32)
        a = (cov * alpha)[..., None]
        out = out * (1 - a) + col * a
    # 密度は重なりで濃くなる。ただし alpha 合成なので飽和はインク色で止まる
    im = Image.fromarray(np.clip(out, 0, 255).astype(np.uint8))
    F.NAVY = INK; F.WHITE = PAPER
    im = F.comp_rel(im, F.aa_mask((W, H), solid_grat(sc)), 0.14, 0.10)
    im = im.resize(SIZE, Image.LANCZOS); im.save(f'{OUT}/{name}.png'); print(name, flush=True)

if __name__ == '__main__':
    bands_img(); marks_img()