# -*- coding: utf-8 -*- """SORAH Visual Loop v2 — Generator の描画部（プロジェクト版 gen-14 時点と同一）。""" import json, math, os import numpy as np, tifffile from PIL import Image, ImageDraw, ImageFilter from pyproj import Transformer from scipy import ndimage RAW = '/home/claude/loop/data/raw' W, H = 1080, 1350 GRID = dict(x0=-3850000.0, y0=5850000.0, px=25000.0, ncol=304, nrow=448) def lab_to_rgb(L, a, b): fy = (L + 16) / 116; fx = fy + a / 500; fz = fy - b / 200 def finv(t): return t**3 if t**3 > 0.008856 else (t - 16/116) / 7.787 X, Y, Z = 0.95047 * finv(fx), 1.0 * finv(fy), 1.08883 * finv(fz) r = 3.2406*X - 1.5372*Y - 0.4986*Z; g = -0.9689*X + 1.8758*Y + 0.0415*Z; bb = 0.0557*X - 0.2040*Y + 1.0570*Z def gam(u): u = max(0, min(1, u)); return 12.92*u if u <= 0.0031308 else 1.055*u**(1/2.4) - 0.055 return tuple(int(round(gam(c) * 255)) for c in (r, g, bb)) def col(L, hue=0, chroma=0): a = chroma * math.cos(math.radians(hue)); b = chroma * math.sin(math.radians(hue)) return lab_to_rgb(L, a, b) _cache = {} def conc(year): if year not in _cache: _cache[year] = tifffile.imread(f'{RAW}/N_{year}09_concentration_v4.0.tif').astype(np.int32) return _cache[year] _land = None def land_polys(): global _land if _land is None: g = json.load(open(f'{RAW}/ne_50m_land.geojson')) tr = Transformer.from_crs('EPSG:4326', 'EPSG:3411', always_xy=True) polys = [] for f in g['features']: geom = f['geometry']; rings = geom['coordinates'] if geom['type'] == 'Polygon' else [r for p in geom['coordinates'] for r in p] for ring in rings: if max(pt[1] for pt in ring) < 30: continue pts = [pt for pt in ring if pt[1] > 20] if len(pts) < 3: continue xs, ys = tr.transform([p[0] for p in pts], [p[1] for p in pts]) polys.append(list(zip(xs, ys))) _land = polys return _land def smooth_mask(mask, sigma_px): """表示用の平滑化。元格子は変えない。ガウス平滑→0.5で再二値化。統計を返す。""" if sigma_px <= 0: return mask, dict(sigma_px=0, max_disp_px=0.0, area_ratio=1.0, components_before=int(ndimage.label(mask)[1]), components_after=int(ndimage.label(mask)[1])) f = ndimage.gaussian_filter(mask.astype(np.float32), sigma_px) sm = f >= 0.5 # 最大変位: 対称差の画素から相手側の境界までの距離の最大 bnd_o = mask & ~ndimage.binary_erosion(mask); bnd_s = sm & ~ndimage.binary_erosion(sm) d_to_orig = ndimage.distance_transform_edt(~bnd_o) # 各画素から元境界までの距離 disp = d_to_orig[bnd_s] # 平滑化後の境界の各点が元境界からどれだけ離れたか md = float(disp.max()) if disp.size else 0.0; p95 = float(np.percentile(disp, 95)) if disp.size else 0.0 lab_o, n_o = ndimage.label(mask); lab_s, n_s = ndimage.label(sm) lost = 0; lost_area = 0 for i in range(1, n_o + 1): comp = lab_o == i if not (comp & sm).any(): lost += 1; lost_area += int(comp.sum()) return sm, dict(sigma_px=float(sigma_px), max_disp_px=md, p95_disp_px=p95, area_ratio=float(sm.sum() / max(1, mask.sum())), components_before=int(n_o), components_after=int(n_s), lost_components=lost, lost_area_px=lost_area, lost_area_ratio=float(lost_area / max(1, mask.sum()))) class Scene: def __init__(self, crop_km=5000, center=(0, 0)): self.wkm = crop_km; self.hkm = crop_km * H / W self.cx, self.cy = center def to_px(self, x, y): return ((x - self.cx) / (self.wkm*1000) + 0.5) * W, (0.5 - (y - self.cy) / (self.hkm*1000)) * H def raster(self, year): a = conc(year) ys, xs = np.mgrid[0:H, 0:W] X = self.cx + (xs / W - 0.5) * self.wkm * 1000 Y = self.cy - (ys / H - 0.5) * self.hkm * 1000 c = ((X - GRID['x0']) / GRID['px']).astype(int); r = ((GRID['y0'] - Y) / GRID['px']).astype(int) ok = (c >= 0) & (c < GRID['ncol']) & (r >= 0) & (r < GRID['nrow']) out = np.full((H, W), -1, np.int32); out[ok] = a[r[ok], c[ok]] return out def render(params, outdir): p = dict(DEFAULT); p.update(params) os.makedirs(outdir, exist_ok=True) ck = p.get('center_km') sc = Scene(p['crop_km'], (ck[0] * 1000.0, ck[1] * 1000.0)) if ck is not None else Scene(p['crop_km'], tuple(p['center'])) grid = sc.raster(p['year']) ice = (grid >= p['ice_threshold'] * 10) & (grid <= 1000) pole = grid == 2510 if p['fill_pole']: ice |= pole land_im = Image.new('L', (W, H), 0); d = ImageDraw.Draw(land_im) for poly in land_polys(): d.polygon([sc.to_px(x, y) for x, y in poly], fill=255) land = np.array(land_im) > 127 ice &= ~land stats = {} spx = p['smooth_km'] * 1000.0 / (p['crop_km'] * 1000.0 / W) if p['smooth_km'] else 0 ice_raw = ice.copy() if spx: ice, stats['ice'] = smooth_mask(ice, spx); ice &= ~land g79 = sc.raster(1979); ice79 = ((g79 >= p['ice_threshold'] * 10) & (g79 <= 1000)); if p['fill_pole']: ice79 |= (g79 == 2510) ice79 &= ~land if spx: ice79, stats['ice79'] = smooth_mask(ice79, spx); ice79 &= ~land edge79 = ice79 & ~ndimage.binary_erosion(ice79, iterations=max(1, p['line79_width'])) if p['line79_coast_gap'] > 0: near_land = ndimage.binary_dilation(land, iterations=p['line79_coast_gap']) edge79 &= ~near_land edge_med = None if p['median_line']: y0, y1 = p['median_years']; cnt = np.zeros((H, W), np.int32); n = 0 for y in range(y0, y1 + 1): gy = sc.raster(y); iy = (gy >= p['ice_threshold'] * 10) & (gy <= 1000) if p['fill_pole']: iy |= (gy == 2510) cnt += iy; n += 1 med = (cnt * 2 >= n) & ~land if spx: med, stats['median'] = smooth_mask(med, spx); med &= ~land edge_med = med & ~ndimage.binary_erosion(med, iterations=max(1, p['median_width'])) if p['line79_coast_gap'] > 0: edge_med &= ~ndimage.binary_dilation(land, iterations=p['line79_coast_gap']) img = np.zeros((H, W, 3), np.uint8) img[:] = col(p['field_L'], p['field_hue'], p['field_C']) if p['ref_mode'] in ('fill', 'both'): fh = p['field_hue'] if p['ref_fill_hue'] is None else p['ref_fill_hue'] fc = p['field_C'] if p['ref_fill_C'] is None else p['ref_fill_C'] img[ice79] = col(p['ref_fill_L'], fh, fc) img[land] = col(p['land_L'], p['land_hue'], p['land_C']) if p['land_mode'] == 'coast' and p['coast_width'] > 0: coast = land & ~ndimage.binary_erosion(land, iterations=p['coast_width'], border_value=1) img[coast] = col(p['coast_L'], p['coast_hue'], p['coast_C']) mode = p['subject_mode'] if mode == 'mass': img[ice] = col(p['ice_L'], p['ice_hue'], p['ice_C']) elif mode == 'hole': img[ice] = col(p['ice_L'], p['ice_hue'], p['ice_C']) elif mode == 'concentration': cval = np.clip(grid, 0, 1000) / 1000.0 lo, hi = p['ice_L_lo'], p['ice_L'] for k in range(8): band = ice & (cval >= k/8) & (cval < (k+1)/8 + (1 if k == 7 else 0)) img[band] = col(lo + (hi - lo) * (k + 0.5) / 8, p['ice_hue'], p['ice_C']) elif mode == 'boundary': pass elif mode == 'stack': # 47年重ね: 各9月で氷だった割合を field→stack_L_hi の明度（彩度は field_C→stack_C_hi）で塗る。 cnt = np.zeros((H, W), np.int32); n = 0 for y in range(p['stack_from'], (p['stack_to'] or p['year']) + 1): gy = sc.raster(y); iy = (gy >= p['ice_threshold'] * 10) & (gy <= 1000) if p['fill_pole']: iy |= (gy == 2510) cnt += iy; n += 1 freq = cnt / n if spx: # 表示用の平滑化は割合の場に一度だけ（各年の格子は変えない） freq = ndimage.gaussian_filter(freq.astype(np.float32), spx) stats['stack'] = dict(sigma_px=float(spx), note='freq field smoothed; per-year grids unchanged') ever = (freq > 0.01) & ~land lo, hi = p['field_L'], p['stack_L_hi']; c_lo, c_hi = p['field_C'], (p['field_C'] if p['stack_C_hi'] is None else p['stack_C_hi']) nb = p['stack_bands'] if nb and nb > 0: for k in range(nb): band = ever & (freq > k/nb) & (freq <= (k+1)/nb) t = (k + 0.5) / nb img[band] = col(lo + (hi - lo) * t, p['field_hue'], c_lo + (c_hi - c_lo) * t) else: # 連続 ys_, xs_ = np.nonzero(ever) # 明度のルックアップ（256段）で連続グラデーション lut = np.array([col(lo + (hi - lo) * t, p['field_hue'], c_lo + (c_hi - c_lo) * t) for t in np.linspace(0, 1, 256)], np.uint8) idx = np.clip((freq[ys_, xs_] * 255).astype(int), 0, 255) img[ys_, xs_] = lut[idx] if p['stack_paint_year']: img[ice] = col(p['ice_L'], p['ice_hue'], p['ice_C']) if p['ice_edge_width'] > 0 and mode != 'boundary': edge = ice & ~ndimage.binary_erosion(ice, iterations=p['ice_edge_width']) img[edge] = col(p['ice_edge_L'], p['ice_hue'], p['ice_C']) if mode == 'boundary': e = ice & ~ndimage.binary_erosion(ice, iterations=p['boundary_width']) img[e] = col(p['ice_L'], p['ice_hue'], p['ice_C']) if p['line79'] and p['ref_mode'] in ('line', 'both'): img[edge79] = col(p['line79_L'], p['line79_hue'], p['line79_C']) if edge_med is not None: img[edge_med] = col(p['median_L'], p['median_hue'], p['median_C']) Image.fromarray(img).save(f'{outdir}/render.png') Image.fromarray(img).resize((320, int(320*H/W)), Image.LANCZOS).save(f'{outdir}/render-320.png') Image.fromarray((ice * 255).astype(np.uint8)).save(f'{outdir}/mask-subject.png') Image.fromarray((land * 255).astype(np.uint8)).save(f'{outdir}/mask-land.png') field = ~ice & ~land Image.fromarray((field * 255).astype(np.uint8)).save(f'{outdir}/mask-field.png') json.dump(p, open(f'{outdir}/variables.json', 'w'), ensure_ascii=False, indent=1) stats['px_km'] = p['crop_km'] / W for k in ('ice','ice79','median'): if k in stats: stats[k]['max_disp_km'] = stats[k]['max_disp_px'] * stats['px_km'] json.dump(stats, open(f'{outdir}/smoothing.json', 'w'), ensure_ascii=False, indent=1) Image.fromarray((ice_raw * 255).astype(np.uint8)).save(f'{outdir}/mask-subject-raw.png') return outdir DEFAULT = dict( year=2025, ice_threshold=15, fill_pole=True, crop_km=5000, center=[0, 0], center_km=None, field_L=91, field_hue=250, field_C=8, land_L=91, land_hue=90, land_C=4, subject_mode='mass', ice_L=99.5, ice_L_lo=60, ice_hue=250, ice_C=0, ice_edge_width=0, ice_edge_L=70, boundary_width=4, stack_from=1979, stack_L_hi=62, line79=True, line79_L=45, line79_hue=55, line79_C=35, line79_width=3, line79_coast_gap=6, ref_mode='line', ref_fill_L=40, ref_fill_hue=None, ref_fill_C=None, median_line=False, median_years=[1981, 2010], median_L=68, median_hue=150, median_C=30, median_width=3, land_mode='fill', coast_L=45, coast_hue=250, coast_C=0, coast_width=1, smooth_km=0, stack_paint_year=True, stack_C_hi=None, stack_bands=8, stack_to=None, # 表示用の平滑化（ガウス sigma, km）。0 で従来と同一。面積・統計は元格子のまま )