"""Bake the raw OpenStreetMap extracts (raw/*.json.gz, see fetch.py) into one compact binary: ../data/prague.bin.

Data (c) OpenStreetMap contributors, ODbL 1.0 (see ../ATTRIBUTION.md). This file is a derived database.

Local frame: metres, x east, z south (three.js: y up, so north is -z), origin at the centre of the core box.
Everything is quantised (Q_NET = 0.5 m for the network, lamps, water and landmarks; Q_BLK = 1 m for building blocks),
simplified (Douglas-Peucker per tier) and written as zigzag LEB128 varints, delta-coded along each polyline.

Container:  'NLC2' | u32 version | u32 nSections | nSections x ( 4-char tag | u32 byteLength | payload )
Every payload except META is a multi-stream: varint nStreams | varint byteLength x nStreams | the streams, each a run of
varints (the layouts below list the fields; each field lives in its own stream, see the encoders).
Sections (all varints unless noted):
  META  utf-8 JSON: origin, quantisation, tier boxes, class table, counts, named places (bridges, squares, landmarks)
  NODE  n | n x (zz dx, zz dz)                       graph nodes in Morton order, Q_NET units
  EDGE  m | m x (da, zz db, code, k, k x (zz dx, zz dz))
        a = prevA + da (edges sorted by a), b = a + db; code = cls | flags << 4; k intermediate points, deltas from
        the previous point starting at node a. flags: 1 bridge, 2 lit=yes, 4 lit=no, 8 link, 16 oneway, 32 area,
        64/128 tier (0 core, 1 mid, 2 outer), 256 one-way runs b -> a
  FAR   n | n x (code, k, k x (zz dx, zz dz))        mid/outer streets as plain polylines, Q_FAR = 2 m, deltas run on
        across lines; code = cls | bridge << 4 | tier << 5
  WATR  nPoly | per poly: nRing | per ring: k | k x (zz dx, zz dz)     (first ring outer; deltas run on across rings)
  BLCK  nBlk  | per block: height m | nRing | per ring: k | k x (zz dx, zz dz)   Q_BLK units
  LAMP  n | n x (zz dx, zz dz, type)                 mapped highway=street_lamp nodes; type 0 other 1 gaslight 2 led
  RIVR  nLine | per line: k | deltas                 the Vltava centre line (waterway=river)
  LMRK  n | per landmark: kind, height m, nRing | rings  (floodlit buildings: 1 church 2 castle 3 civic 4 tower 5 theatre/museum)
"""
import gzip
import json
import math
import os
import struct
import sys
import time
from collections import Counter, defaultdict

from shapely.geometry import LineString, MultiPolygon, Polygon, box
from shapely.ops import linemerge, polygonize, unary_union
from shapely.strtree import STRtree

HERE = os.path.dirname(os.path.abspath(__file__))
RAW = os.path.join(HERE, 'raw')
OUT = os.path.join(HERE, '..', 'data', 'prague.bin')

CORE = (50.0700, 14.3900, 50.1150, 14.4650)
MID = (50.0280, 14.3350, 50.1530, 14.5280)
OUTER = (49.9700, 14.2300, 50.1950, 14.6500)
LAT0 = (CORE[0] + CORE[2]) / 2
LON0 = (CORE[1] + CORE[3]) / 2
KX = math.cos(math.radians(LAT0)) * 111320.0
KZ = 110574.0
Q_NET = 0.5
Q_BLK = 1.0
SIMPLIFY = {0: 0.6, 1: 1.5, 2: 3.0}  # per tier, metres

# class table (index = code)
CLASSES = ['motorway', 'trunk', 'primary', 'secondary', 'tertiary', 'residential', 'living_street', 'pedestrian',
           'service', 'busway', 'footway', 'rail']
HW = {
    'motorway': 0, 'motorway_link': 0, 'trunk': 1, 'trunk_link': 1, 'primary': 2, 'primary_link': 2,
    'secondary': 3, 'secondary_link': 3, 'tertiary': 4, 'tertiary_link': 4, 'residential': 5, 'unclassified': 5,
    'living_street': 6, 'pedestrian': 7, 'service': 8, 'busway': 9, 'footway': 10, 'path': 10, 'cycleway': 10,
    'steps': 10,
}
TIER_CLASSES = {0: set(range(12)), 1: {0, 1, 2, 3, 4, 5, 6, 7}, 2: {0, 1, 2, 3, 4}}


def proj(lat, lon):
    return ((lon - LON0) * KX, -(lat - LAT0) * KZ)


def in_box(lat, lon, b):
    return b[0] <= lat <= b[2] and b[1] <= lon <= b[3]


def box_m(b):
    x0, z1 = proj(b[0], b[1])
    x1, z0 = proj(b[2], b[3])
    return (x0, z0, x1, z1)


# ---------------------------------------------------------------------------------------------------------------- io
class W:
    def __init__(self):
        self.b = bytearray()

    def u(self, v):
        v = int(v)
        assert v >= 0, v
        while True:
            c = v & 0x7F
            v >>= 7
            if v:
                self.b.append(c | 0x80)
            else:
                self.b.append(c)
                return

    def z(self, v):
        v = int(v)
        self.u((v << 1) if v >= 0 else ((-v << 1) - 1))


class MS:
    """several varint streams in one section (gzip likes like with like): varint nStreams | varint len x n | data"""
    def __init__(self, n):
        self.s = [W() for _ in range(n)]

    def u(self, i, v):
        self.s[i].u(v)

    def z(self, i, v):
        self.s[i].z(v)

    def bytes(self):
        h = W()
        h.u(len(self.s))
        for w in self.s:
            h.u(len(w.b))
        return bytes(h.b) + b''.join(bytes(w.b) for w in self.s)


def q(v, s):
    return int(round(v / s))


def dp(pts, tol):
    """Douglas-Peucker keeping the endpoints (pts: list of (x, z))."""
    if len(pts) <= 2:
        return pts
    keep = [False] * len(pts)
    keep[0] = keep[-1] = True
    stack = [(0, len(pts) - 1)]
    while stack:
        i, j = stack.pop()
        ax, az = pts[i]
        bx, bz = pts[j]
        dx, dz = bx - ax, bz - az
        L2 = dx * dx + dz * dz
        best, bi = -1.0, -1
        for k in range(i + 1, j):
            px, pz = pts[k]
            if L2 < 1e-12:
                d = math.hypot(px - ax, pz - az)
            else:
                t = max(0.0, min(1.0, ((px - ax) * dx + (pz - az) * dz) / L2))
                d = math.hypot(px - ax - t * dx, pz - az - t * dz)
            if d > best:
                best, bi = d, k
        if best > tol:
            keep[bi] = True
            stack.append((i, bi))
            stack.append((bi, j))
    return [p for p, k in zip(pts, keep) if k]


def morton(x, z):
    def spread(v):
        v &= 0xFFFF
        v = (v | (v << 8)) & 0x00FF00FF
        v = (v | (v << 4)) & 0x0F0F0F0F
        v = (v | (v << 2)) & 0x33333333
        v = (v | (v << 1)) & 0x55555555
        return v
    return spread(x + 32768) | (spread(z + 32768) << 1)


def load(name):
    with gzip.open(os.path.join(RAW, f'{name}.json.gz')) as f:
        return json.load(f)['elements']


# ------------------------------------------------------------------------------------------------------------ streets
def way_class(t):
    h = t.get('highway')
    if h in HW:
        if h == 'footway' and t.get('footway') in ('sidewalk', 'crossing', 'traffic_island'):
            return None
        # paths only where they are lit (park and embankment walks); the rest are dark or duplicate a street
        if h in ('footway', 'path', 'cycleway', 'steps') and t.get('lit') not in ('yes', '24/7', 'automatic'):
            return None
        if h == 'service' and t.get('service') in ('driveway', 'drive-through', 'emergency_access', 'parking_aisle')                 and t.get('lit') != 'yes':
            return None
        return HW[h]
    r = t.get('railway')
    if r == 'rail' and t.get('service') not in ('siding', 'spur', 'yard') and t.get('usage') != 'industrial':
        return 11
    return None


def way_flags(t, tier, closed):
    f = 0
    if t.get('bridge') in ('yes', 'viaduct', 'aqueduct', 'boardwalk', 'cantilever', 'movable') or t.get('man_made') == 'bridge':
        f |= 1
    lit = t.get('lit')
    if lit in ('yes', '24/7', 'automatic', 'sunset-sunrise', 'dusk-dawn', 'interval', 'limited'):
        f |= 2
    elif lit == 'no':
        f |= 4
    if (t.get('highway') or '').endswith('_link'):
        f |= 8
    if t.get('oneway') in ('yes', '1', 'true') or t.get('junction') == 'roundabout':
        f |= 16
    if closed and t.get('area') == 'yes':
        f |= 32
    return f | (tier << 6)


def build_network(tiers):
    ways = {}
    for tier, elems in tiers:
        for e in elems:
            if e['type'] != 'way' or e['id'] in ways:
                continue
            t = e.get('tags', {})
            layer = t.get('layer', '0')
            if t.get('tunnel') in ('yes', 'culvert') or t.get('location') == 'underground':
                continue
            if layer.lstrip('-').isdigit() and int(layer) < 0 and t.get('tunnel') != 'building_passage':
                continue
            c = way_class(t)
            if c is None or c not in TIER_CLASSES[tier]:
                continue
            g = e.get('geometry') or []
            ids = e.get('nodes') or []
            if len(g) < 2 or len(ids) != len(g):
                continue
            ways[e['id']] = (tier, c, t, ids, g)
    use = Counter()
    for tier, c, t, ids, g in ways.values():
        for k, nid in enumerate(ids):
            use[nid] += 2 if k in (0, len(ids) - 1) else 1
    coord = {}
    edges = []  # [a_id, b_id, cls, flags, pts]
    for wid, (tier, c, t, ids, g) in ways.items():
        closed = ids[0] == ids[-1]
        flags = way_flags(t, tier, closed)
        seg_ids = [ids[0]]
        seg_pts = []
        for k, (nid, p) in enumerate(zip(ids, g)):
            if p is None:
                seg_ids, seg_pts = [], []
                continue
            xy = proj(p['lat'], p['lon'])
            coord[nid] = xy
            if not seg_pts:
                seg_ids = [nid]
                seg_pts = [xy]
                continue
            seg_pts.append(xy)
            seg_ids.append(nid)
            if use[nid] >= 2 or k == len(ids) - 1:
                if len(seg_pts) >= 2 and seg_ids[0] != nid or len(seg_pts) > 2:
                    edges.append([seg_ids[0], nid, c, flags, seg_pts])
                seg_ids = [nid]
                seg_pts = [xy]
    # merge chains through degree-2 nodes with identical attributes
    inc = defaultdict(list)
    for i, e in enumerate(edges):
        inc[e[0]].append(i)
        inc[e[1]].append(i)
    alive = [True] * len(edges)
    merged = 0
    for nid in list(inc.keys()):
        L = [i for i in inc[nid] if alive[i]]
        if len(L) != 2 or L[0] == L[1]:
            continue
        i, j = L
        ei, ej = edges[i], edges[j]
        if ei[2] != ej[2] or ei[3] != ej[3] or ei[0] == ei[1] or ej[0] == ej[1]:
            continue
        # orient ei to end at nid and ej to start at nid
        pi = ei[4] if ei[1] == nid else ei[4][::-1]
        ai = ei[0] if ei[1] == nid else ei[1]
        pj = ej[4] if ej[0] == nid else ej[4][::-1]
        bj = ej[1] if ej[0] == nid else ej[0]
        if ai == bj:
            continue  # would make a loop edge; keep the node
        oneway = ei[3] & 16
        if oneway and not ((ei[1] == nid) == (ej[0] == nid)):
            continue  # opposite one-way directions
        new = [ai, bj, ei[2], ei[3], pi + pj[1:]]
        if oneway and ei[1] != nid:
            new = [bj, ai, ei[2], ei[3], (pi + pj[1:])[::-1]]  # both were reversed: keep the traffic direction
        edges[i] = new
        alive[j] = False
        inc[bj] = [i if x == j else x for x in inc[bj]]
        inc[nid] = []
        merged += 1
    edges = [e for e, a in zip(edges, alive) if a]
    return edges, coord, merged


def encode_network(edges):
    # nodes present in edges
    nodes = {}
    for e in edges:
        for nid, p in ((e[0], e[4][0]), (e[1], e[4][-1])):
            nodes[nid] = p
    qn = {nid: (q(p[0], Q_NET), q(p[1], Q_NET)) for nid, p in nodes.items()}
    order = sorted(qn.keys(), key=lambda n: morton(*qn[n]))
    index = {nid: i for i, nid in enumerate(order)}
    # NODE streams: 0 dx, 1 dz
    ms = MS(2)
    ms.u(0, len(order))
    px = pz = 0
    for nid in order:
        x, z = qn[nid]
        ms.z(0, x - px)
        ms.z(1, z - pz)
        px, pz = x, z
    node_bytes = ms.bytes()
    # edges: simplify, orient a <= b, sort
    rows = []
    npts = 0
    for a, b, c, f, pts in edges:
        tier = (f >> 6) & 3
        sp = dp(pts, SIMPLIFY[tier])
        ia, ib = index[a], index[b]
        if ia > ib:
            ia, ib, sp = ib, ia, sp[::-1]
            if f & 16:
                f |= 256  # one-way against the stored a -> b direction
        rows.append((ia, ib, c, f, sp))
    rows.sort(key=lambda r: (r[0], r[1]))
    # EDGE streams: 0 da, 1 db, 2 code, 3 k, 4 dx, 5 dz
    ms = MS(6)
    ms.u(0, len(rows))
    pa = 0
    for ia, ib, c, f, sp in rows:
        ms.u(0, ia - pa)
        pa = ia
        ms.z(1, ib - ia)
        ms.u(2, c | (f << 4))
        mid = sp[1:-1]
        ms.u(3, len(mid))
        x0, z0 = qn[order[ia]]
        for p in mid:
            x, z = q(p[0], Q_NET), q(p[1], Q_NET)
            ms.z(4, x - x0)
            ms.z(5, z - z0)
            x0, z0 = x, z
        npts += len(mid)
    return node_bytes, ms.bytes(), len(order), len(rows), npts


# ---------------------------------------------------------------------------------------------------------- far field
Q_FAR = 2.0


def encode_far(edges):
    """mid/outer tiers: plain polylines per class (not routable), chains re-merged ignoring lit/oneway flags."""
    by = defaultdict(list)
    for a, b, c, f, pts in edges:
        tier = (f >> 6) & 3
        by[(c, f & 1, tier)].append(LineString(pts))
    lines = []
    for (c, br, tier), ls in by.items():
        m = linemerge(ls)
        parts = list(m.geoms) if hasattr(m, 'geoms') else [m]
        for p in parts:
            if p.length < 25:
                continue
            pts = dp(list(p.coords), 4.0 if tier == 1 else 6.0)
            lines.append((c | (br << 4) | (tier << 5), pts))
    # order the lines so each starts near where the previous one ended (short first deltas): greedy nearest start
    # (either end; a line may be reversed) through a 200 m grid
    G = 200.0
    grid = defaultdict(list)
    for i, (code, pts) in enumerate(lines):
        for end in (0, 1):
            p = pts[0] if end == 0 else pts[-1]
            grid[(int(p[0] // G), int(p[1] // G))].append((i, end))
    used = [False] * len(lines)
    ordered = []
    cur = (0.0, 0.0)
    remaining = len(lines)
    while remaining:
        cx, cz = int(cur[0] // G), int(cur[1] // G)
        best = None
        r = 0
        while best is None and r < 400:
            cand = []
            for gx in range(cx - r, cx + r + 1):
                for gz in (cz - r, cz + r) if r else (cz,):
                    cand += grid.get((gx, gz), [])
                if r:
                    pass
            if r:
                for gz in range(cz - r + 1, cz + r):
                    cand += grid.get((cx - r, gz), []) + grid.get((cx + r, gz), [])
            for i, end in cand:
                if used[i]:
                    continue
                p = lines[i][1][0] if end == 0 else lines[i][1][-1]
                d = (p[0] - cur[0]) ** 2 + (p[1] - cur[1]) ** 2
                if best is None or d < best[0]:
                    best = (d, i, end)
            r += 1
        if best is None:
            break
        _, i, end = best
        used[i] = True
        remaining -= 1
        code, pts = lines[i]
        if end == 1:
            pts = pts[::-1]
        ordered.append((code, pts))
        cur = pts[-1]
    lines = ordered
    # FAR streams: 0 code, 1 k, 2 dx, 3 dz
    ms = MS(4)
    ms.u(0, len(lines))
    px = pz = 0
    nv = 0
    for code, pts in lines:
        ms.u(0, code)
        ms.u(1, len(pts))
        for x, z in pts:
            X, Z = q(x, Q_FAR), q(z, Q_FAR)
            ms.z(2, X - px)
            ms.z(3, Z - pz)
            px, pz = X, Z
        nv += len(pts)
    return ms.bytes(), len(lines), nv


# -------------------------------------------------------------------------------------------------------------- water
def ring_pts(ring):
    return list(ring.coords)[:-1]


def assemble_relation(e):
    outers, inners = [], []
    for m in e.get('members', []):
        g = m.get('geometry')
        if not g or m.get('type') != 'way' or any(p is None for p in g):
            continue
        ls = LineString([proj(p['lat'], p['lon']) for p in g])
        (inners if m.get('role') == 'inner' else outers).append(ls)
    def polys(lines):
        if not lines:
            return []
        return list(polygonize(linemerge(lines) if len(lines) > 1 else lines))
    try:
        po = [p if p.is_valid else p.buffer(0) for p in polys(outers)]
        pi = [p if p.is_valid else p.buffer(0) for p in polys(inners)]
        if not po:
            return None
        shape = unary_union(po)
        if pi:
            shape = shape.difference(unary_union(pi))
        return shape if shape.is_valid else shape.buffer(0)
    except Exception:  # noqa: BLE001 - a broken multipolygon is skipped
        return None


def water_shape(elems_list, clip):
    parts = []
    seen = set()
    for elems in elems_list:
        for e in elems:
            if e['id'] in seen and e['type'] == 'relation':
                continue
            t = e.get('tags', {})
            is_water = t.get('natural') == 'water' or t.get('waterway') in ('riverbank', 'dock')
            if not is_water or t.get('water') in ('fountain', 'wastewater', 'basin') or t.get('covered') == 'yes':
                continue
            if t.get('tunnel'):
                continue
            if e['type'] == 'relation':
                seen.add(e['id'])
                s = assemble_relation(e)
                if s is not None and not s.is_empty:
                    parts.append(s)
            elif e['type'] == 'way':
                g = e.get('geometry') or []
                if len(g) >= 4 and e['nodes'][0] == e['nodes'][-1] and all(p is not None for p in g):
                    p = Polygon([proj(p['lat'], p['lon']) for p in g])
                    if not p.is_valid:
                        p = p.buffer(0)
                    parts.append(p)
    u = unary_union(parts).intersection(clip)
    u = u.buffer(0)
    geoms = list(u.geoms) if isinstance(u, MultiPolygon) else [u]
    return [g for g in geoms if g.area > 400]


def flatten(polys, extra=None):
    out, ex = [], []
    for i, pg in enumerate(polys):
        parts = list(pg.geoms) if hasattr(pg, 'geoms') else [pg]
        for p in parts:
            if isinstance(p, Polygon) and not p.is_empty:
                out.append(p)
                ex.append(extra[i] if extra is not None else None)
    return out, (ex if extra is not None else None)


def encode_polys(polys, qs, tol, extra=None):
    """streams: 0 structure (n | per poly: extra..., nRing | per ring: k), 1 dx, 2 dz (deltas run on across rings)"""
    polys, extra = flatten(polys, extra)
    ms = MS(3)
    out = []
    for i, pg in enumerate(polys):
        pg = pg.simplify(tol, preserve_topology=True)
        if pg.is_empty or not isinstance(pg, Polygon):
            continue
        ext = pg.exterior
        if len(ext.coords) < 4:
            continue
        rings = [ext] + [r for r in pg.interiors if len(r.coords) >= 4]
        out.append((extra[i] if extra is not None else [], rings))
    ms.u(0, len(out))
    px = pz = 0
    nv = 0
    for ex, rings in out:
        for v in ex:
            ms.u(0, v)
        ms.u(0, len(rings))
        for r in rings:
            pts = ring_pts(r)
            ms.u(0, len(pts))
            for x, z in pts:
                X, Z = q(x, qs), q(z, qs)
                ms.z(1, X - px)
                ms.z(2, Z - pz)
                px, pz = X, Z
            nv += len(pts)
    return ms.bytes(), nv


# ---------------------------------------------------------------------------------------------------------- buildings
def building_height(t):
    h = t.get('height') or t.get('building:height')
    try:
        if h:
            return max(3.0, min(120.0, float(str(h).split()[0].replace(',', '.'))))
    except ValueError:
        pass
    lv = t.get('building:levels')
    try:
        if lv:
            return max(3.0, min(120.0, float(str(lv).split(';')[0].replace(',', '.')) * 3.2 + 1.5))
    except ValueError:
        pass
    return None


LANDMARK_KIND = {'church': 1, 'cathedral': 1, 'chapel': 1, 'synagogue': 1, 'castle': 2, 'palace': 2}
# offices, town halls, ministries and the old wastewater plant are not the floodlit monuments of a night skyline
NOT_LANDMARK = ('úřad', 'magistrát', 'ministerstv', 'čistírn', 'finanční', 'banka', 'hotel')


def landmark_kind(t):
    b = t.get('building')
    hist = t.get('historic')
    name = (t.get('name') or '').lower()
    if any(w in name for w in NOT_LANDMARK) or b in ('office', 'government', 'commercial', 'industrial'):
        return 0
    if b in LANDMARK_KIND:
        return LANDMARK_KIND[b]
    if t.get('amenity') == 'place_of_worship':
        return 1
    if hist in ('castle', 'palace'):
        return 2
    if hist in ('tower', 'city_gate') or t.get('man_made') == 'tower':
        return 4
    if t.get('amenity') in ('theatre', 'concert_hall') or t.get('tourism') == 'museum':
        return 5
    if t.get('wikidata') and b in ('civic', 'government', 'public', 'train_station') :
        return 3
    return 0


def buildings(elems, clip):
    polys, heights, marks = [], [], []
    for e in elems:
        t = e.get('tags', {})
        if 'building' not in t or t.get('building') in ('no', 'roof', 'construction') or t.get('location') == 'underground':
            continue
        shape = None
        if e['type'] == 'way':
            g = e.get('geometry') or []
            if len(g) < 4 or any(p is None for p in g):
                continue
            shape = Polygon([proj(p['lat'], p['lon']) for p in g])
        elif e['type'] == 'relation':
            shape = assemble_relation(e)
        if shape is None or shape.is_empty:
            continue
        if not shape.is_valid:
            shape = shape.buffer(0)
        if shape.area < 12:
            continue
        polys.append(shape)
        heights.append(building_height(t))
        k = landmark_kind(t)
        if k and k != 3 and shape.area > 150:
            marks.append((k, building_height(t) or 22.0, shape, t.get('name')))
    t0 = time.time()
    u = unary_union([p.buffer(0.6, join_style=2) for p in polys]).buffer(-0.6, join_style=2)
    u = u.intersection(clip)
    blocks = list(u.geoms) if isinstance(u, MultiPolygon) else [u]
    blocks = [b for b in blocks if b.area > 40]
    tree = STRtree(polys)
    bh = []
    for b in blocks:
        idx = tree.query(b)
        acc = wsum = 0.0
        for i in idx:
            h = heights[i]
            if h is None:
                continue
            a = polys[i].area
            acc += h * a
            wsum += a
        bh.append(int(round(acc / wsum)) if wsum > 0 else 0)
    print(f'  blocks: {len(blocks)} from {len(polys)} buildings in {time.time() - t0:.1f} s; landmarks {len(marks)}')
    return blocks, bh, marks


# ------------------------------------------------------------------------------------------------------------- places
def places(elems):
    out = {'bridges': [], 'squares': []}
    by_name = defaultdict(list)
    # squares: every way named "... náměstí" (Prague's squares are mapped as the streets round them), plus Ungelt
    sq = defaultdict(list)
    for e in elems:
        t = e.get('tags', {})
        name = t.get('name') or ''
        if e['type'] != 'way' or not e.get('geometry') or not t.get('highway'):
            continue
        if ('náměstí' in name.lower() or name == 'Ungelt') and not name.lower().startswith('vestibul'):
            sq[name] += [proj(p['lat'], p['lon']) for p in e['geometry'] if p]
    for e in elems:
        t = e.get('tags', {})
        if e['type'] != 'way' or not e.get('geometry'):
            continue
        name = t.get('name') or ''
        if t.get('bridge') and t.get('highway') and 'most' in name.lower() and t.get('highway') not in ('footway', 'steps', 'path', 'cycleway', 'service'):
            by_name[name] += [proj(p['lat'], p['lon']) for p in e['geometry'] if p]
        if t.get('place') == 'square' and name and name not in sq:
            sq[name] += [proj(p['lat'], p['lon']) for p in e['geometry'] if p]
    for name, pts in by_name.items():
        # principal axis
        cx = sum(p[0] for p in pts) / len(pts)
        cz = sum(p[1] for p in pts) / len(pts)
        sxx = sum((p[0] - cx) ** 2 for p in pts)
        szz = sum((p[1] - cz) ** 2 for p in pts)
        sxz = sum((p[0] - cx) * (p[1] - cz) for p in pts)
        ang = 0.5 * math.atan2(2 * sxz, sxx - szz)
        d = [math.cos(ang), math.sin(ang)]
        proj1 = [(p[0] - cx) * d[0] + (p[1] - cz) * d[1] for p in pts]
        a, b = min(proj1), max(proj1)
        out['bridges'].append({'name': name, 'c': [round(cx, 1), round(cz, 1)],
                               'a': [round(cx + d[0] * a, 1), round(cz + d[1] * a, 1)],
                               'b': [round(cx + d[0] * b, 1), round(cz + d[1] * b, 1)]})
    out['bridges'].sort(key=lambda b: b['c'][1])
    for name, pts in sorted(sq.items()):
        if len(pts) < 2:
            continue
        cx = sum(p[0] for p in pts) / len(pts)
        cz = sum(p[1] for p in pts) / len(pts)
        out['squares'].append({'name': name, 'c': [round(cx, 1), round(cz, 1)]})
    return out


# --------------------------------------------------------------------------------------------------------------- main
def main():
    t0 = time.time()
    core = load('core')
    mid = load('mid')
    midwater = load('midwater')
    outer = load('outer')
    print(f'loaded in {time.time() - t0:.1f} s')

    all_edges, coord, merged = build_network([(0, core), (1, mid), (2, outer)])
    edges = [e for e in all_edges if (e[3] >> 6) & 3 == 0]
    node_b, edge_b, nN, nE, nP = encode_network(edges)
    print(f'network (core, routable): {nN} nodes, {nE} edges, {nP} inner points (merged {merged}); {len(node_b) + len(edge_b)} bytes')
    far_b, nFl, nFv = encode_far([e for e in all_edges if (e[3] >> 6) & 3 != 0])
    print(f'far field: {nFl} polylines, {nFv} vertices, {len(far_b)} bytes')

    cls_count = Counter(CLASSES[e[2]] + ('/bridge' if e[3] & 1 else '') for e in edges)

    core_m = box_m(CORE)
    mid_m = box_m(MID)
    outer_m = box_m(OUTER)
    water = water_shape([core, midwater], box(*mid_m))
    water_b, nWv = encode_polys(water, Q_NET, 1.0)
    print(f'water: {len(water)} polygons, {nWv} vertices, {len(water_b)} bytes')

    blocks, bh, marks = buildings(core, box(*core_m))
    blk_b, nBv = encode_polys(blocks, Q_BLK, 3.0, extra=[[h] for h in bh])
    print(f'blocks: {len(blocks)}, {nBv} vertices, {len(blk_b)} bytes')

    marks.sort(key=lambda m: -m[2].area)
    marks = marks[:160]
    mk_polys = [m[2] for m in marks]
    mk_b, nMv = encode_polys(mk_polys, Q_NET, 1.0, extra=[[m[0], int(round(m[1]))] for m in marks])
    print(f'landmarks: {len(marks)}, {nMv} vertices, {len(mk_b)} bytes')

    # lamps
    lamps = []
    for e in core:
        t = e.get('tags', {})
        if e['type'] == 'node' and t.get('highway') == 'street_lamp':
            x, z = proj(e['lat'], e['lon'])
            lt = t.get('lamp_type') or t.get('light:method') or ''
            typ = 1 if lt in ('gaslight', 'gas') else 2 if lt in ('led', 'LED') else 0
            lamps.append((q(x, Q_NET), q(z, Q_NET), typ))
    lamps.sort(key=lambda l: morton(l[0], l[1]))
    ms = MS(3)  # 0 dx (and n), 1 dz, 2 type
    ms.u(0, len(lamps))
    px = pz = 0
    for x, z, typ in lamps:
        ms.z(0, x - px)
        ms.z(1, z - pz)
        ms.u(2, typ)
        px, pz = x, z
    lamp_b = ms.bytes()
    print(f'lamps: {len(lamps)} ({Counter(l[2] for l in lamps)}), {len(lamp_b)} bytes')

    # river centre line
    rl = []
    for e in core + midwater + mid:
        t = e.get('tags', {})
        if e['type'] == 'way' and t.get('waterway') == 'river' and e.get('geometry'):
            rl.append(LineString([proj(p['lat'], p['lon']) for p in e['geometry'] if p]))
    rivers = []
    if rl:
        m = linemerge(rl)
        m = m.intersection(box(*mid_m))
        rivers = list(m.geoms) if hasattr(m, 'geoms') else [m]
        rivers = [r for r in rivers if r.length > 200]
    ms = MS(3)  # 0 structure, 1 dx, 2 dz
    ms.u(0, len(rivers))
    px = pz = 0
    for r in rivers:
        pts = dp(list(r.coords), 3.0)
        ms.u(0, len(pts))
        for x, z in pts:
            X, Z = q(x, Q_NET), q(z, Q_NET)
            ms.z(1, X - px)
            ms.z(2, Z - pz)
            px, pz = X, Z
    riv_b = ms.bytes()
    print(f'river lines: {len(rivers)}, {len(riv_b)} bytes')

    pl = places(core)
    pl['landmarks'] = [{'name': m[3], 'kind': m[0], 'c': [round(m[2].centroid.x, 1), round(m[2].centroid.y, 1)]} for m in marks if m[3]][:80]
    meta = {
        'name': 'Praha - Vltava bend (Staré Město, Josefov, Letná, Holešovice, Karlín, Nové Město)',
        'attribution': '© OpenStreetMap contributors (ODbL 1.0), https://www.openstreetmap.org/copyright',
        'fetched': time.strftime('%Y-%m-%d'),
        'origin': {'lat': LAT0, 'lon': LON0, 'kx': KX, 'kz': KZ},
        'frame': 'metres; x east, z south (north = -z); y up',
        'quant': {'net': Q_NET, 'blk': Q_BLK, 'far': Q_FAR},
        'tiers': {'core': [round(v, 1) for v in core_m], 'mid': [round(v, 1) for v in mid_m], 'outer': [round(v, 1) for v in outer_m]},
        'tierBoxesLatLon': {'core': CORE, 'mid': MID, 'outer': OUTER},
        'classes': CLASSES,
        'counts': {'nodes': nN, 'edges': nE, 'edgePoints': nP, 'water': len(water), 'blocks': len(blocks), 'lamps': len(lamps), 'landmarks': len(marks), 'byClass': dict(cls_count)},
        'places': pl,
    }
    meta_b = json.dumps(meta, ensure_ascii=False, separators=(',', ':')).encode('utf-8')
    sections = [(b'META', meta_b), (b'NODE', node_b), (b'EDGE', edge_b), (b'FAR ', far_b), (b'WATR', water_b), (b'BLCK', blk_b),
                (b'LAMP', lamp_b), (b'RIVR', riv_b), (b'LMRK', mk_b)]
    out = bytearray(b'NLC2' + struct.pack('<II', 2, len(sections)))
    for tag, b in sections:
        out += tag + struct.pack('<I', len(b)) + b
    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    with open(OUT, 'wb') as f:
        f.write(out)
    gz = gzip.compress(bytes(out), 9)
    print(f'wrote {OUT}: {len(out)} bytes, gzip {len(gz)} bytes')
    for tag, b in sections:
        print(f'  {tag.decode()}: {len(b):>8} raw  {len(gzip.compress(b, 9)):>8} gzip')
    print(f'total {time.time() - t0:.1f} s')


if __name__ == '__main__':
    main()