# -*- coding: utf-8 -*- """森林被覆（JAXA PALSAR-2 FNF）を SORAH の二色語法で描く。Z2 の規範を世界地図へ移植したもの。 規範の移植: - 紙＝ネイビー L12 / 陸＝L18 / 現象（森）＝白 L97 まで 8 段の帯 - 帯は「画素の中の森林被覆率」。面積は元格子（0.025°）で数え、表示だけ縮める - 線・文字は 4 倍スーパーサンプリング、下地からの相対シフトで置く - 階層 森 ＞ 陸・地名 ＞ 回帰線・赤道 """ import json, math, os, sys import numpy as np sys.path.insert(0, '/home/claude/loop/lib'); sys.path.insert(0, '/home/claude/dv/cards') import render3 as r, g1 from PIL import Image, ImageDraw from pyproj import Transformer DATA = '/home/claude/forest/data' W, H = 1080, 1350; SS = 4 NAVY = np.array(r.col(12, 270, 16), np.float32); WHITE = np.array([247, 247, 250], np.float32) ON_LIGHT = np.array(r.col(58, 265, 12), np.float32); ON_DARK = np.array(r.col(40, 265, 14), np.float32) FIELD = dict(L=12, hue=270, C=16) # 紙＝海 LAND = dict(L=18, hue=268, C=14) # 陸 FOREST_HI = 97 # 最上位の帯（森が画素を満たす） tr = Transformer.from_crs('EPSG:4326', 'EPSG:8857', always_xy=True) inv = Transformer.from_crs('EPSG:8857', 'EPSG:4326', always_xy=True) # ---------- 場（森林被覆率） ---------- _cache = {} def frac(year, ds=8, dense_only=False): """0.025° の分類から、ds×ds ブロックごとの森林被覆率と陸の割合を返す。面積は元格子で数える。""" k = (year, ds, dense_only) if k in _cache: return _cache[k] a = np.load(f'{DATA}/fnf_{year}_L1.npy') h, w = a.shape; h2, w2 = h // ds, w // ds a = a[:h2*ds, :w2*ds].reshape(h2, ds, w2, ds) forest = ((a == 1) if dense_only else ((a == 1) | (a == 2))).sum((1, 3)).astype(np.float32) landc = ((a >= 1) & (a <= 3)).sum((1, 3)).astype(np.float32) f = np.divide(forest, np.maximum(landc, 0.25*ds*ds), dtype=np.float32) # 海に接する画素で少数の陸が全体を代表しないように下限を置く _cache[k] = (f, landc / (ds * ds)); return _cache[k] def _pc_fwd(lon, lat): lon = np.asarray(lon, float); lat = np.asarray(lat, float) return lon * 111320.0, lat * 111320.0 def _pc_inv(x, y): x = np.asarray(x, float); y = np.asarray(y, float) return x / 111320.0, y / 111320.0 class Scene: """投影は Equal Earth（等積）か正距円筒。extent は投影座標の m。""" def __init__(self, bbox_lonlat=(-180, -60, 180, 84), size=(W, H), lon0=0.0, pad=0.02, proj='ee'): self.proj = proj self.fwd = (lambda a, b: tr.transform(a, b)) if proj == 'ee' else _pc_fwd self.invf = (lambda a, b: inv.transform(a, b)) if proj == 'ee' else _pc_inv self.lon0 = lon0; self.W, self.H = size xs, ys = [], [] for lon in np.linspace(bbox_lonlat[0], bbox_lonlat[2], 181): for lat in np.linspace(bbox_lonlat[1], bbox_lonlat[3], 91): x, y = self.fwd(self._wrap(lon), lat); xs.append(float(x)); ys.append(float(y)) x0, x1, y0, y1 = min(xs), max(xs), min(ys), max(ys) cx, cy = (x0 + x1) / 2, (y0 + y1) / 2 sx = (x1 - x0) * (1 + pad) / self.W; sy = (y1 - y0) * (1 + pad) / self.H self.s = max(sx, sy); self.cx, self.cy = cx, cy def _wrap(self, lon): return ((lon - self.lon0 + 180) % 360) - 180 def to_px(self, lon, lat): x, y = self.fwd(self._wrap(lon), lat); x = float(x); y = float(y) return ((x - self.cx) / self.s + self.W / 2, self.H / 2 - (y - self.cy) / self.s) def lonlat_grid(self): ys, xs = np.mgrid[0:self.H, 0:self.W] X = self.cx + (xs - self.W / 2 + 0.5) * self.s; Y = self.cy - (ys - self.H / 2 + 0.5) * self.s lon, lat = self.invf(X.ravel(), Y.ravel()) lon = np.array(lon).reshape(self.H, self.W); lat = np.array(lat).reshape(self.H, self.W) ok = np.isfinite(lon) & np.isfinite(lat) & (np.abs(lat) <= 90) # 投影の外側（丸い縁の外）は逆投影が値を返してしまうので、往復させて一致しない画素を落とす x2, y2 = self.fwd(np.where(ok, lon, 0).ravel(), np.where(ok, lat, 0).ravel()) x2 = np.array(x2).reshape(self.H, self.W); y2 = np.array(y2).reshape(self.H, self.W) ok &= (np.abs(x2 - X) < self.s) & (np.abs(y2 - Y) < self.s) lon = np.where(ok, ((lon + self.lon0 + 180) % 360) - 180, 0); lat = np.where(ok, lat, 0) return lon, lat, ok def sample(grid, lon, lat, ok): h, w = grid.shape c = np.clip(((lon + 180) / 360 * w).astype(int), 0, w - 1) rr = np.clip(((90 - lat) / 180 * h).astype(int), 0, h - 1) out = grid[rr, c]; return np.where(ok, out, 0) def sample_bilinear(grid, lon, lat, ok): """場を双線形で読む（最近傍だと帯の境目が元格子の形に折れる）。""" h, w = grid.shape fx = np.clip((lon + 180) / 360 * w - 0.5, 0, w - 1.001); fy = np.clip((90 - lat) / 180 * h - 0.5, 0, h - 1.001) x0 = fx.astype(int); y0 = fy.astype(int); tx = fx - x0; ty = fy - y0 g = np.nan_to_num(grid.astype(np.float32)) out = (g[y0, x0]*(1-tx)*(1-ty) + g[y0, x0+1]*tx*(1-ty) + g[y0+1, x0]*(1-tx)*ty + g[y0+1, x0+1]*tx*ty) return np.where(ok, out, 0) # ---------- 合成（Z2 と同じ規範） ---------- def lstar(rgb): x = rgb.astype(np.float32) / 255; y = 0.2126*x[..., 0] + 0.7152*x[..., 1] + 0.0722*x[..., 2] return 116 * np.cbrt(y) - 16 def aa_mask(size, fn): big = Image.new('L', (size[0]*SS, size[1]*SS), 0); d = ImageDraw.Draw(big); fn(d, SS) return np.array(big.resize(size, Image.LANCZOS)).astype(np.float32) / 255 def comp_rel(im, alpha, kd, kl): base = np.array(im).astype(np.float32); L = lstar(base) toward = np.where((L > 55)[..., None], base*(1-kl) + NAVY*kl, base*(1-kd) + WHITE*kd) a = alpha[..., None]; return Image.fromarray(np.clip(base*(1-a) + toward*a, 0, 255).astype(np.uint8)) def comp_label(im, alpha, tone_light=ON_LIGHT, tone_dark=ON_DARK): base = np.array(im).astype(np.float32); L = lstar(base) m = alpha > 0.02 if m.any(): from scipy import ndimage lab, n = ndimage.label(m) tone = np.zeros_like(base) for i in range(1, n + 1): sel = lab == i; mean = L[sel].mean() tone[sel] = tone_light if mean > 55 else tone_dark else: tone = np.zeros_like(base) a = alpha[..., None]; return Image.fromarray(np.clip(base*(1-a) + tone*a, 0, 255).astype(np.uint8)) def dashed(d, pts, width, s, dash): on, off = dash[0]*s, dash[1]*s; draw = True; acc = 0.0 for (x0, y0), (x1, y1) in zip(pts[:-1], pts[1:]): seg = math.hypot(x1-x0, y1-y0) if seg > 400*s or seg == 0: draw = True; acc = 0.0; continue t = 0.0 while t < seg: step = min(seg - t, (on if draw else off) - acc) if draw: d.line([(x0+(x1-x0)*t/seg, y0+(y1-y0)*t/seg), (x0+(x1-x0)*(t+step)/seg, y0+(y1-y0)*(t+step)/seg)], fill=255, width=int(width*s)) t += step; acc += step if acc >= (on if draw else off) - 1e-6: draw = not draw; acc = 0.0 # ---------- 地理レイヤー ---------- _land = None def land_polys(): global _land if _land is None: g = json.load(open('/home/claude/loop/data/raw/ne_50m_land.geojson')); 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 len(ring) >= 3: polys.append(ring) _land = polys return _land _lg = None def land_grid(ppd=40): """陸のマスクを等緯度経度の格子に一度だけ焼く（日付変更線の継ぎ目を作らないため）。南極は FNF の観測外なので除く。""" global _lg if _lg is None: p = f'{DATA}/land_{ppd}.npy' if os.path.exists(p): _lg = np.load(p) else: big = Image.new('L', (360*ppd, 180*ppd), 0); d = ImageDraw.Draw(big) for ring in land_polys(): if max(pt[1] for pt in ring) < -60: continue d.polygon([((lon+180)*ppd, (90-lat)*ppd) for lon, lat in ring], fill=255) _lg = np.array(big) > 127; np.save(p, _lg) return _lg def land_mask(sc): lon, lat, ok = sc.lonlat_grid() return sample(land_grid().astype(np.float32), lon, lat, ok) def grat(sc, lats=(-60, -30, 0, 30, 60), dash=(2.0, 5.0)): def fn(d, s): for lat in lats: dashed(d, [(x*s, y*s) for x, y in (sc.to_px(lon, lat) for lon in np.arange(-180, 180.1, 1))], 1.0, s, dash) for lon in range(-180, 180, 30): dashed(d, [(x*s, y*s) for x, y in (sc.to_px(lon, lat) for lat in np.arange(-88, 88.1, 0.5))], 1.0, s, dash) return fn def tropics(sc): def fn(d, s): for lat in (23.44, -23.44, 0): dash = (1.0, 2.5) if lat else (3.0, 3.0) dashed(d, [(x*s, y*s) for x, y in (sc.to_px(lon, lat) for lon in np.arange(-180, 180.1, 1))], 1.0, s, dash) return fn def text_fn(sc, item, size, track=0.16, dot=0, left=False, dy=0): def fn(d, s): f = g1.font('Regular', int(size*s)); name, lon, lat = item; x, y = sc.to_px(lon, lat); x *= s; y = y*s + dy*s if dot: d.ellipse([x-dot*s, y-dot*s, x+dot*s, y+dot*s], fill=255); y -= size*s*0.55 if left: wsum = sum(d.textlength(ch, font=f) + size*track*s for ch in name); x -= (dot+6)*s + wsum else: x += (dot+6)*s for ch in name: d.text((x, y), ch, font=f, fill=255); x += d.textlength(ch, font=f) + size*track*s return fn REGIONS = [('AMAZON', -62, -6), ('CONGO', 22, -4), ('BORNEO', 113, -1.5), ('SIBERIA', 98, 60), ('CANADA', -100, 54), ('CONGO BASIN', 22, -4)][:5]