# -*- coding: utf-8 -*-
"""航空路線（OpenFlights routes.dat, ODbL）。線×時間なし。大圏で結び、太さ＝その路線を飛ぶ航空会社数の階級。"""
import sys, os
import numpy as np, pandas as pd
sys.path.insert(0, '/home/claude/forest'); sys.path.insert(0, '/home/claude/quakes')
import forest1 as F, sylvania as S, q1 as Q
from PIL import Image, ImageDraw
OUT = '/home/claude/transport/out'; os.makedirs(OUT, exist_ok=True)
ap = pd.read_csv('/home/claude/transport/data/airports.dat', header=None, usecols=[0,4,6,7], names=['id','iata','lat','lon'])
rt = pd.read_csv('/home/claude/transport/data/routes.dat', header=None, usecols=[0,2,4], names=['airline','src','dst'])
pos = ap.dropna(subset=['iata']).set_index('iata')[['lat','lon']]
rt = rt[rt.src.isin(pos.index) & rt.dst.isin(pos.index)]
# 無向の路線にまとめ、航空会社数を数える
key = rt.apply(lambda r: tuple(sorted((r.src, r.dst))), axis=1)
cnt = key.value_counts()
pairs = pd.DataFrame(list(cnt.index), columns=['a','b']); pairs['n'] = cnt.values
print('routes rows', len(rt), 'undirected', len(pairs), 'n quantiles', pairs.n.quantile([.5,.9,.99]).tolist())
ORANGE = np.array(F.r.col(66, 58, 62), np.float32); NAVY = np.array(F.r.col(12, 270, 16), np.float32)
WHITE = np.array([255,255,255], np.float32); LANDC = np.array([238,238,235], np.float32)
BBOX = (-180, -58, 180, 83); SIZE = (1350, 1080); SS = 3
CLS = [(1, 2, 0.35, 0.35), (2, 4, 0.55, 0.55), (4, 8, 0.8, 0.75), (8, 999, 1.2, 0.9)]

def gc(lat1, lon1, lat2, lon2, n=48):
    p1, l1, p2, l2 = map(np.radians, (lat1, lon1, lat2, lon2))
    d = 2*np.arcsin(np.sqrt(np.sin((p2-p1)/2)**2 + np.cos(p1)*np.cos(p2)*np.sin((l2-l1)/2)**2))
    if d < 1e-9: return None
    f = np.linspace(0, 1, n)
    A = np.sin((1-f)*d)/np.sin(d); B = np.sin(f*d)/np.sin(d)
    x = A*np.cos(p1)*np.cos(l1) + B*np.cos(p2)*np.cos(l2); y = A*np.cos(p1)*np.sin(l1) + B*np.cos(p2)*np.sin(l2); z = A*np.sin(p1) + B*np.sin(p2)
    return np.degrees(np.arctan2(z, np.hypot(x, y))), np.degrees(np.arctan2(y, x))

def build(name, paper, landc, ink, lon0=0.0, cls=CLS, grat=True):
    sc = F.Scene(BBOX, (SIZE[0]*SS, SIZE[1]*SS), proj='pc', pad=0.02, lon0=lon0); W, H = sc.W, sc.H
    lon, lat, ok = sc.lonlat_grid()
    land = F.sample(Q.land_grid_all().astype(np.float32), lon, lat, ok) > 0.5
    out = np.zeros((H, W, 3), np.float32); out[:] = paper; out[land] = landc
    def px(la, lo):
        lo = ((lo - lon0 + 180) % 360) - 180
        return (lo*111320.0 - sc.cx)/sc.s + W/2, H/2 - (la*111320.0 - sc.cy)/sc.s
    for lo_, hi_, wpx, al in cls:
        sel = pairs[(pairs.n >= lo_) & (pairs.n < hi_)]
        cv = Image.new('L', (W, H), 0); d = ImageDraw.Draw(cv)
        for a, b in zip(sel.a, sel.b):
            g = gc(pos.at[a,'lat'], pos.at[a,'lon'], pos.at[b,'lat'], pos.at[b,'lon'])
            if g is None: continue
            x, y = px(g[0], g[1]); pts = list(zip(x, y))
            run = [pts[0]]
            for p0, p1 in zip(pts[:-1], pts[1:]):
                if abs(p1[0]-p0[0]) > W/2:
                    if len(run) > 1: d.line(run, fill=255, width=max(1, int(round(wpx*SS))))
                    run = [p1]
                else: run.append(p1)
            if len(run) > 1: d.line(run, fill=255, width=max(1, int(round(wpx*SS))))
        a = (np.array(cv).astype(np.float32)/255*al)[..., None]; out = out*(1-a) + ink*a
    im = Image.fromarray(np.clip(out,0,255).astype(np.uint8))
    if grat:
        from nyt import solid_grat
        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); return im

if __name__ == '__main__':
    build('A1', WHITE, LANDC, NAVY)                    # 白紙・ネイビー（紙から最も遠い色）
    build('A2', NAVY, np.array(F.r.col(18,268,14),np.float32), ORANGE)   # ネイビー紙・オレンジ（人の活動）
    build('A3', WHITE, LANDC, NAVY, cls=[(2,4,0.45,0.5),(4,8,0.8,0.75),(8,999,1.2,0.9)])   # 1 社だけの路線を間引く