# -*- coding: utf-8 -*-
"""陰影起伏の地（釜田さんの参照画像の系統）。ETOPO 2022 から計算。海は平坦な灰、陸は標高で明るく、北西からの光で陰影。"""
import sys; sys.path.insert(0,'/home/claude/quakes'); sys.path.insert(0,'/home/claude/forest')
import numpy as np, forest1 as F, sylvania as S, q1 as Q
from scipy import ndimage
from PIL import Image
DEM = np.load('/home/claude/quakes/dem/etopo_0p05.npy')

def hillshade(z, cell_deg=0.05, az=315, alt=45, vex=3.0):
    lat = np.linspace(90-cell_deg/2, -90+cell_deg/2, z.shape[0])
    dx = 111320.0*cell_deg*np.cos(np.radians(lat))[:,None]; dy = 111320.0*cell_deg
    gy, gx = np.gradient(z*vex, dy, axis=0), np.gradient(z*vex, axis=1)/dx
    slope = np.arctan(np.hypot(gx, gy)); aspect = np.arctan2(-gx, gy)
    azr, altr = np.radians(az), np.radians(alt)
    hs = np.sin(altr)*np.cos(slope) + np.cos(altr)*np.sin(slope)*np.cos(azr - aspect)
    return np.clip(hs, 0, 1).astype(np.float32)

def base(sc, lon, lat, ok, land, ocean_L=48, land_lo=68, land_hi=92, shade=0.55, vex=3.0, ocean_shade=0.0, a=0.0, b=0.0):
    """L の合成: 海＝ocean_L（必要なら海底の陰影を薄く）、陸＝標高で land_lo→land_hi、陰影を乗算的に。"""
    z = np.where(DEM > 0, DEM, 0)
    hs = hillshade(z, vex=vex); hs = ndimage.gaussian_filter(hs, 0.8)
    hs_all = hillshade(DEM, vex=vex) if ocean_shade > 0 else None
    zz = F.sample_bilinear(DEM, lon, lat, ok); h = F.sample_bilinear(hs, lon, lat, ok)
    zl = np.clip(zz, 0, 6000) / 6000.0
    L = land_lo + (land_hi - land_lo) * np.sqrt(zl)
    L = L * (1 - shade) + L * shade * (0.35 + 0.65 * h * 1.6).clip(0, 1.6)   # 陰影は明度を上下両方に
    L = np.where(land, L, ocean_L)
    if ocean_shade > 0:
        ho = F.sample_bilinear(hs_all, lon, lat, ok); L = np.where(land, L, ocean_L * (1 - ocean_shade) + ocean_L * ocean_shade * (0.4 + 1.2 * ho)).astype(np.float32)
    L = np.clip(L, 0, 100)
    Ls = np.linspace(0, 100, 1001); tab = np.array([S.from_lab(np.array([l, a, b])) for l in Ls])
    return tab[np.clip((L*10).astype(int), 0, 1000)].astype(np.float32)

# ---- F5 系（反転: 海が明るい紙、陸が暗い。光は北西）----
_DEM33 = None
def dem33():
    global _DEM33
    if _DEM33 is None: _DEM33 = np.load('/home/claude/quakes/dem/etopo_0p033.npy')
    return _DEM33
def f5_base(sc, lon, lat, ok, land, ocean_L=90.0, ocean_amp=0.0, ocean_depth=0.0, vex_land=10.0, vex_sea=10.0):
    """陸: L=(26+14√z)·(0.25+1.5·h)、海: L=ocean_L·(1 − ocean_amp·(1−h_sea)) − ocean_depth·(深さ/6000)。"""
    D = dem33()
    zz = F.sample_bilinear(D, lon, lat, ok); zl = np.clip(zz, 0, 6000)/6000.0
    hl = F.sample_bilinear(hillshade(np.where(D > 0, D, 0), cell_deg=1/30, vex=vex_land), lon, lat, ok)/0.72
    Lland = np.clip((26 + 14*np.sqrt(zl)) * (0.25 + 1.5*hl), 6, 70)
    if ocean_amp > 0 or ocean_depth > 0:
        hsea = F.sample_bilinear(hillshade(np.where(D < 0, D, 0), cell_deg=1/30, vex=vex_sea), lon, lat, ok)/0.72
        dep = np.clip(-zz, 0, 6000)/6000.0
        Lsea = ocean_L * (1 - ocean_amp*(1 - hsea)) - ocean_depth*dep
    else: Lsea = np.full_like(zz, ocean_L)
    L = np.where(land, Lland, np.clip(Lsea, 0, 100))
    Ls = np.linspace(0, 100, 1001); tab = np.array([S.from_lab(np.array([l, 0, 0])) for l in Ls])
    return tab[np.clip((L*10).astype(int), 0, 1000)].astype(np.float32)